MultI-messenger probe of Cosmic Ray Origins: MICRO project

. We present the milestones achieved by the MICRO (MultI-messenger probe of Cosmic Ray Origins) project. This multi-institute project consists of a study of bursting astrophysical sources as candidate sources for ultra-high-energy cosmic rays (UHECRs). We aim at identifying source classes that correlate best with existing observational data (arrival direction distribution, energy spectrum and primary mass composition) and to constrain the origin of UHECRs on the basis of the largest datasets acquired to this day. The study of the secondary ﬂuxes of neutrinos and gamma rays will provide a powerful test of the most suitable astrophysical scenarios. For these achievements, a public software enabling a joint ﬁt of both composition and ﬂux data as a function of energy and direction is being developed. In addition, a modelling of transient sources is foreseen, including the treatment of hadronic interactions within the sources, which represents a novelty for UHECR simulation codes. This contribution will highlight the ﬁrst results of the project and detail upcoming software releases, which will beneﬁt the astroparticle community at large.


Introduction
Identifying the sources of Ultra-High Energy Cosmic Rays (UHECR) is one of the most pressing questions in highenergy astrophysics. Despite the advent of high-statistics and high-quality data, most prominently obtained by the Pierre Auger Observatory, providing an astrophysical explanation to the observed features at Earth is not trivial. To shed light on the origin of UHECRs, many studies, such as for instance refs [1][2][3] and references therein, have considered a simple cosmological model, assuming all the sources identical, pointlike and steady in time. MICRO (MultI-messenger probe of Cosmic Ray Origins) is a project which aims to go beyond the conventional modelling of steady sources of UHECRs. A novelty of this approach is a complete study of bursting sources starting from the modelling of selected source classes, including hadronic interactions within the source and over the propagation down to Earth, to predict the UHECR sky as a function of energy and mass composition. The objectives specifically addressed in the MICRO project will provide important answers to the leading questions of identifying the sources of UHECRs: (i) how do burst-like signatures such as Gamma Rays Bursts (GRBs) or Active Galactic Nuclei (AGN) flares fit the cosmic-ray data, (ii) how can we constrain the 3D distribution of sources from available UHECR observables, and (iii) could astrophysical high-energy neutrinos, some high-energy gamma rays, and UHECRs come from the same bursting sources. Two different subjects are addressed within the MICRO project: first the development of a public software, enabling a joint fit of both UHECRs composition and flux data as functions of energy and direction; second the modelling the source environment, studying how this can affect UHECR ejection.

Towards a public software for UHECR observables
One of the milestones of the project is to release a public software capable of providing acceptable solutions to the astrophysical origin of UHECRs. The idea is to define it an astrophysical scenario, assuming bursting sources, following a certain distribution in redshift. The injected flux for each mass at the source can be written as: where E 0 is a reference energy, ε tot j · k governs the normalization for the mass j, which represents the injected energy, e.g. per stellar mass per rate , R cut is a cutoff rigidity, γ is the spectral index at escape from the sources and f cut (E i , Z · R cut ) is the cutoff function. In order to connect properties at the sources with observables at Earth, the extra-galactic propagation has to be performed. This is done considering Monte Carlo codes, such as SimProp [4] or CRPropa [5]. All the events generated are stored in a normalized 5D tensor T E det , A det , z E inj , A inj . In this way, the probability that a particle with energy E inj and mass A inj at a redshift z would be detected with an energy E det and a mass A det is performed. After properly taking into account the possible interactions of UHECRs with the cosmic microwave background (CMB) and the extra-galactic background light (EBL) during the extra-galactic propagation, the propagated flux at Earth can be calculated and then compared to the experimental data. In Fig. 1 the expected and the experimental spectra are shown, as well as the first two moments of the depth of the shower maxima X max distributions.
As in similar works [1,2], it was found that it is pos- Figure 1. Top: Fluxes of nuclei with primary injected particles as labelled. Curves with different colours show the sum of the flux of primaries with a given mass number and all secondaries produced by the same nuclear species. The fit is performed only at energies above the ankle (10 18.7 eV). Experimental data are the flux measurements by the Pierre Auger Collaboration [6]. Bottom: black points are the average and standard deviation of the X max distribution at each energy, as measured by the Pierre Auger collaboration ; coloured lines are predicted by Sibyll2.3d [7] for various pure compositions (thin lines) and for the composition predicted by our model at the best fit (thick brown line).
sible to have a good description of both UHECR energy spectrum and mass composition in terms of a hard spectral index γ < 2, and a low maximum energy of protons 2 · 10 18 eV. An advantage of using a Monte Carlo simulations for computing the extra-galactic propagation is that the secondary productions can be easily computed. Neutrinos propagate through the inter-galactic space undergoing only adiabatic energy losses due to the expansion of the Universe, hence their energy spectrum at Earth, unlike that of UHECRs, is expected to be also affected by the presence of sources at large redshifts and by their cosmological evolution. For this reason, the flux of predicted cosmogenic neutrinos associated with the best fit results of each chosen scenario can be compared with the measured fluxes (or, at higher energies, with the experimental upper bounds) and such comparison could possibly constrain the cosmological evolution of sources. This is done in Fig. 2, where the expected neutrino flux is compared to the IceCube flux [8] and different computed limits for cosmogenic neutrinos.
The expected neutrino flux could serve as multimessenger constraint for the sources of UHECRs. A neutrino flux exceeding the measured one improves the constraining capability with respect to UHECR data only and then it can improve the search for the best-fit Figure 2. The predicted flux of cosmogenic neutrinos for a combined fit above the ankle, compared to the IceCube neutrino flux [8], the GRAND sensitivity after three years of operation [9] and to the limits for the cosmogenic neutrinos by the Pierre Auger [10] and IceCube [11] collaborations. parameters at the source.

Organization of the code
The code has been divided in several different modules in Python, which separate the different steps of the fit.
• Tensor: This module handles the reading of the tensor, stored in a numpy's compress array format for each injected mass. This module allows us to take into account the extra-galactic propagation: UHECRs travel through the extra-galactic environment and interact with photon backgrounds changing their energy and, in the case of nuclei, their species too. This module computes the weights that will be used to transform the spectra at the source to the expected spectra at Earth.
• Spectrum: This module contains methods to simulate the emission of UHECR nuclei from sources. Nuclei of mass A and charge Z are injected in the extra-galactic space with a power-law energy spectrum with spectral index γ and rigidity cutoff R cut , as shown in eq. (1). Starting from the injected spectra at the source and considering a defined evolution of the source to a certain z max , methods within this class can be used to compute the expected energy spectrum at Earth up to the atmosphere. In addition, this module reads in the experimental spectrum, to provide both the model and the data to the module minimizer, which deals with the minimization of the deviance.
• Mass and Gumbel: These classes deal with the effects of the cosmic-ray cascade production in the atmosphere which, due to stochastic effects, introduce shower to shower fluctuations in particular for the slant depth of the shower maximum (X max ). The expected X max distributions are calculated making use of the generalized Gumbel probability distribution functions [12], and then corrected for the detector effects. A great advantage of using this parametric function is that it allows us to evaluate the model X max distribution for any mixture of nuclei without the need to simulate showers for each primary. In this analysis all the nuclei from A = 1 up to A = 56 are taken into account. The spectra are unaffected by atmosphere, so they are not managed in this class. Alternatively, it is possible to read and compute the mean mass ln A and its variance at Earth, or the first two moments of the X max distributions, using parametrizations obtained by CONEX showers [13]. As in the module spectrum, in mass the reading of the experimental mass composition indicator, as well as the computation of the expected one, is performed.
• Minimizer: This class implements the methods for calculating the deviances of the spectrum and mass composition. The deviances are calculated comparing experimental data with model expectations, where the latter depend on fit parameters. The X max deviance calculation is performed at each energy bin of the corresponding X max distribution.
Others modules used in the code are: Draw, which contains all the utilities to produce the plots, Xmax_tools, which is a library containing all the parametrizations for computing X max and σ(X max ). A module called Skymaps for taking into account the arrival directions in the combined fit is under development.

Including the hadronic interaction in UHECRs code
A key part of the MICRO project regards the inclusion of the hadronic interaction, i.e. interactions with of UHECRs with other hadrons, in the UHECR propagation codes. Hadronic interactions are negligible for the extra-galactic propagation, due to the low density of the gas in the extra-galactic medium; however, for propagation inside the source environment, they are a crucial quantity to be taken into account. The computation of hadronic process is a new entry in UHECR propagation codes. In fact, as shown in [14], where hadronic interactions have been implemented in SimProp, hadronic interactions generate more secondary particles than photo-hadronic process, consequently more photons, neutrinos and a larger number of light nuclear fragments. This different composition of propagated nuclei produces, in turn, a change in the evolution of the nuclear cascade inside the source environment. Fig. 3 illustrates the mass distributions of nuclei escaping from a prototype starburst galaxy. In particular, nuclear fragments are shown considering a single type of target: the photo-hadronic-only scenario (blue) is shown separately from the hadronic-only one (red). It can be seen that, while the hadronic scenario produces efficiently all lighter nuclei, the photo-hadronic scenario does not produce efficiently intermediate mass (C, N, O) nuclei.
An effort to include hadronic interaction in CRPropa is ongoing. This is aimed to be done using impy [15]. A good advantage of calling impy consists in directly comparing different hadronic interaction models and their predictions, for example, the secondary production, as shown in Fig. 4; in this way it is possible to study the systematic uncertanties related to the choice of one model. An  . Percentage of particles produced in p-p interactions according to two different hadronic interaction model: PHOJET [16] and SIBYLL [7]. Credits: L. Morejon.
ongoing work inside the CRPropa team concerns the extension of this module for nuclei interactions. Finally, a comparison between the implementation of the hadronic interactions between SimProp and CRPropa is foreseen, which will benefit the astroparticle community at large.

Modelling the source
A goal of the project consists in modelling the source environment of potentially interesting class of sources. Within this spirit, many collaborators are involved in modelling of multimessenger flares from blazars [17]. The idea is to understand the light curves of gamma and radio fluxes and their temporal correlations to putative neutrino emission from the blazar PKS 1502+106 . This requires detailed nu-  To properly take into account the physics on sub-kpc scales, several ingredients are required: 1) The possibility to select custom photon fields for the chosen source, e.g. the thermal field emitted by the accretion disc or the synchrotron-radiation field emitted by relativistic electrons; 2) To have time-dependent interaction rates realizing the temporal evolution during propagation along the jet axis; 3) To include hadronic interactions; 4) To tag the secondary emission, in order to resolve origin of secondaries by processes; 5) the spatial resolution of (turbulent and ordered) magnetic fields in sources are important quantities to determine. It was chosen to use CRPropa, using the same framework of the model exploited in [18]; in this way it was possible to have time resolved multimessenger signatures of relativistic plasmoids, i.e. coherent structure of plasma and magnetic fields, as shown in Fig. 5.

Future perspectives
The MICRO project is reaching its prefixed milestones and aims to have a useful impact in the scientific community, especially releasing a public software for UHECR analysis and including hadronic interactions in UHECR propagation codes. The great novelties consists in developing a cosmological scenario for describing the UHECR observables using transient sources. Injecting bursts of cosmic rays will be important, since it will provide insights on the UHECR spectrum prior to interactions within the cosmic environment and thus it will constrain the sources. In addition, we aim at accounting for the in-source propagation scenario in the individual treatment of AGN/GRBs, with a broad scope of the possible emission scenarios. As for the emission of secondaries, the different photon fields, magnetic fields, and gas distributions obviously play an important role as well. As another important product of the MICRO project, we aim at predicting the flux of neutral secondaries produced within the sources: this will be made possible by implementing the in-source propagation scenario, by determining whether the flaring scenarios being explored are consistent with current limits from Fermi-LAT, Ice-Cube, and Auger for diffuse emissions [19], and whether neutrino and γ-ray point-like emission could arise from the propagation of nuclei within the source environment.