Systematic parameter space study for the UHECR originfrom GRBs in models with multiple internal shocks

We scrutinize the paradigm that conventional long-duration Gamma-Ray Bursts (GRBs) are the dominant source of the ultra-high energy cosmic rays (UHECRs) within the internal shock scenario by describing UHECR spectrum and composition and by studying the predicted (source and cosmogenic) neutrino fluxes. Since it has been demonstrated that the stacking searches for astrophysical GRB neutrinos strongly constrain the parameter space in single-zone models, we focus on the dynamics of multiple collisions for which different messengers are expected to come from different regions of the same object. We propose a model which can describe both stochastic and deterministic engines, which we study in a systematic way. We find that GRBs can indeed describe the UHECRs for a wide range of different model assumptions with comparable quality; the required heavy mass fraction at injection is found to be larger than 70% (95% CL). We demonstrate that the post-dicted (from UHECR data) neutrino fluxes from sources and UHECR propagation are indeed below the current sensitivities but will be reached by the next generation of experiments. We finally critically review the required source energetics with the specific examples found in this study.


INTRODUCTION
Gamma-Ray Bursts (GRBs) have been proposed to be powerful enough to describe Ultra-High Energy Cosmic Rays (UHECRs). The interactions of cosmic rays with photons in the source can lead to substantial neutrino production during the prompt emission phase of gamma rays, see e.g. Waxman and Bahcall (1997); Hümmer et al. (2012). In this work, we focus on conventional long-duration GRBs with isotropic luminosities around 10 52.5 erg s −1 (Gruber et al. 2014). The searches for GRB neutrinos are highly sensitive because the directional, timing and energy information can be used to suppress background. Therefore, the neutrino flux limits from the IceCube Observatory for GRBs are the best among all potential source classes (Abbasi et al. 2012;Aartsen et al. 2017). Besides that, the observation of a temporally and spatially coincident neutrino emission is one of the few ways to locate transient cosmic ray accelerators such as GRBs. So far, no source neutrinos from GRBs have been observed. The interaction of cosmic rays with the cos-E-mail: jonas.heinze@desy.de mic backgrounds during propagation to Earth can also lead to a "cosmogenic" neutrino flux (Berezinsky and Zatsepin 1969;Aab et al. 2017; but no cosmogenic neutrinos have been found either (Aartsen et al. 2018;Aab et al. 2019). A model-dependent correlation analysis of the UHECR arrival directions with source catalogs, performed by the Pierre Auger Observatory for Active Galactic Nuclei and Starburst Galaxies, indicates an intermediate-scale anisotropy related to the distribution of nearby UHECR sources (Aab et al. 2018). According to this result, only a limited fraction of the observed UHECR flux can be attributed to known starburst galaxies or active galactic nuclei, leaving the main contribution to the UHECR flux unconstrained. GRBs may be therefore a conceivable option for the remaining part.
In the present study, we follow the hypothesis that GRBs are the dominant sources of the UHECRs and describe the UHECR spectrum and composition within the internal shock scenario (Kobayashi et al. 1997;Daigne and Mochkovitch 1998). It has been shown that simple one-zoneemission models are in tension with the neutrino flux limits for most of the parameter space, both for the case of a pure proton composition  and for a heavier primary composition (Biehl et al. 2018). This argument follows from the over-production of source neutrinos due to high radiation densities in the source (Biehl et al. 2018) for luminosities and collision radii typically assumed for GRBs -which points towards alternative sources, such as low-luminosity GRBs (Zhang et al. 2018;. Multi-collision models, such as the one discussed in this work, exhibit a lower neutrino flux (Bustamante et al. 2014;Globus et al. 2015a), which is typically dominated by cosmicray interactions close to the photosphere where the radiation densities are high. If only a small fraction of the collisions occurs there, the overall neutrino flux is lower than in the single zone model for the same proton injection luminosity. For a discussion of the production regions of cosmic messengers, the impact of the engine properties on the different messengers and light curves in different energy bands, and the dependence on the collision model for the shells, see Bustamante et al. (2014Bustamante et al. ( , 2017; Rudolph et al. (2019).
In contrast to one-zone models, multi-collision models add the possibility of studying different jet structures and corresponding light curves. A multi-peaked light curve can be generated by a smooth Lorentz factor variation in the outflow, short-time variability by additional random variations (stochasticity) in the Lorentz factor profile. A discussion on different outflow structures and corresponding light curves can be found in e.g. Daigne and Mochkovitch (1998); Bosnjak et al. (2009)). Also Bustamante et al. (2017) investigate different Lorentz factor profiles in connection with neutrino and gamma-ray light curves for proton-loaded jets, however without fitting the source parameters to data. Most GRB multi-collision models do not jointly describe the UHECR spectrum and composition -except Globus et al. (2015a,b), who performed a complete investigation of the multi-collision model in the context of UHECR including nuclei for a smooth, continuous outflow that corresponds to a single-peaked light curve without a short-time variability or an intermittent engine. They draw a self-consistent picture for one set of parameters and one collision model, with a neutrino flux prediction close to the current stacking limit. No substantial progress on the subject has been made since then, owing to the computational complexity of the problem.
In Section 2 of this paper, we present our multi-collision GRB model for different outflow patterns accounting for a smooth temporal profile (a ramp-up of the Lorentz factor towards later times) in combination with a stochastic engine behavior, mainly based on an extension of the techniques developed in our previous works. Section 3 describes the systematic scan over the engine parameters and how these are estimated from UHECR data. For the first time, we obtain contours on the allowed parameter space from UHECR data and discuss the result in detail based on four representative cases. Section 4 is dedicated to the multi-messenger signatures, i.e. the light-curves and the neutrino flux predictions. We aim to conclude whether the absence of associations of IceCube neutrinos with GRBs Abbasi et al. (2012) really disfavors, or even excludes, those as a dominant source of UHECRs when a more sophisticated model for the multimessenger production is used.

Multi-collision dynamics
We follow the model first presented in Kobayashi et al. (1997) using the implementation described on detail in Bustamante et al. (2017); Rudolph et al. (2019) in order to describe the relativistic outflow of the GRB jet. We briefly review the model here and discuss a few modifications and improvements to the model implemented since then. We also discuss the key differences to the collision model from Globus et al. (2015a), which is based on Daigne and Mochkovitch (1998).
The relativistic outflow with non-uniform density and velocity profiles is approximated by a series of plasma shells with distinct masses mi, Lorentz factors Γi and widths li, separated by distances di. Due to the gradient in velocity, the more rapid shells will eventually catch up to slower shells. A collision between a rapid shell (index r) and a slow shell (index s) forms a new, merged shell (index m). Due to energy conservation, the energy dissipated during the collision is given by the difference of kinetic energy before and after the collision: (1) Momentum and mass conservation implies that the Lorentz factor of the merged shell Γm is given by: Γm mrΓr + msΓs mr/Γr + ms/Γs .
In Daigne and Mochkovitch (1998), most of the energy is dissipated during the initial phase where the less massive shell sweeps up a mass equal to its own in the other shell, such that Γ diss = √ Γr · Γs. We follow this assumption here, diverting from Bustamante et al. (2017). Individual and merged shells continue to propagate in the fireball until they collide with other shells or reach the circumburst medium, where they are taken out of the simulation. The circumburst medium is assumed to dominate the outflow dynamics at a distance R circumburst 5.5 ·10 11 km from the central engine.
In each collision particles are accelerated and emitted. We assume that the timescale of emission δtem is given by the time the reverse shock takes to cross the rapid shell. Since both shells get compressed during the collision process, the width of the merged shell is reduced during the collision (lm,C < ls + lr). The a full derivation for the formulas for δtem and lm,C is contained in Kobayashi et al. (1997); Bustamante et al. (2017). In a different scenario, the internal energy remaining in the shell after the collision will result in an expansion of the merged shell, which is discussed in more detail in Rudolph et al. (2019). Here we assume that the merged shell after the collision recovers the width of a single shell before the collision (l m,final = ls = lr). At all times, the shells contained in the fireball thus share the same width. In Globus et al. (2015a), shell widths do not play a role since the energy densities are computed assuming a constant kinetic wind luminosity. This approach is strictly valid only for the first round of collisions and roughly corresponds to the assumption of a constant shell width.
The final emitted particle spectra are computed by summing over all single collision spectra. A collision that occurs at time tC at a distance RC from the engine will start to be observed at time t obs = tC − RC/c. We calculate the time-dependent flux from each collision assuming a 'Fast Rise and Exponential Decay' (FRED) shape normalized to the total gamma-ray luminosity of the collision (see Kobayashi et al. (1997); Bustamante et al. (2017) for details).
Initial setup The simulations start with 1000 shells with equal energies E kin = 10 52 erg, 1 shell widths l = c · 0.002 s and shell separations d = c · 0.002 s. The innermost shell is assumed to be located at radius Rmin = 10 3 km. After the fireball simulation completes, the kinetic energy E kin is re-normalized to ensure a gamma-ray output Eγ 10 53 erg over all collisions.
The initial shell distribution follows a log-normal distribution around a deterministic (temporal) Lorentz factor profile Γ 0,k that assigns higher velocities to shells emitted at later stages of the fireball evolution: On top of that, the Lorentz factor of the k-th shell is assumed to be stochastically distributed around that profile by where x is sampled from a Gaussian P (x)dx = exp(−x 2 )/ √ 2πdx. The spread AΓ describes the strength of the fluctuations, and the values of Γmin and Γmax control the mean and dynamic range of the entire distribution. Note that the bulk Lorentz factor, which (given equal energy) is calculated as depends non-linearly on the parameters Γmin, Γmax and AΓ.
Benchmark engine profiles In order to simplify the discussion and the visualization of our results, we define four distinct initial benchmark shell configurations. With the formulae from the last section those correspond to four choices of Γmin, Γmax and AΓ that define the strength of the rampup and the stochasticity of the engine.
The stochastic fluctuations of the Lorentz factors are considered to be related to the short time variability on top of a pulsed light curve and have been studied in the fireball framework in the past, see e.g. Daigne and Mochkovitch (1998). Globus et al. (2015a) examine the case of a disciplined engine, adding no stochasticity to the ramp-up structure. In addition to the former works, we study the impact of stochastic engine behavior and the strength of the ramp-up in the multi-messenger context including a fit to the observed UHECRs.
The four cases are shown in Fig. 1 and, sorted by increasing stochasticity, correspond to: (1) Strong ramp-up, no stochasticity (SR-0S): Strong ramp-up of the mean Lorentz factor towards later times (thus smaller radii), no stochasticity.
(2) Strong ramp-up, low stochasticity (SR-LS): Strong ramp-up of the mean Lorentz factor, some stochasticity. The shell distribution is mainly dominated by the ramp-up.
(3) Weak ramp-up, medium stochasticity (WR-MS): Moderate ramp-up of the mean Lorentz factor, medium stochasticity. The profile and the stochastic features are both pronounced.
(4) Weak ramp-up, high stochasticity (WR-HS): Weak ramp-up of the mean Lorentz factor toward later times, high stochasticity. The shell distribution is dominated by the stochasticity, the structure of the profile Γ 0,k is barely visible.
The initial Lorentz factors are illustrated in Fig. 1 where we also display the corresponding values of Γmin, Γmax and AΓ.

Radiation model
The radiation produced by each collision is computed independently using the time-dependent radiation code Neu-CosmA (Biehl et al. 2018). The total spectrum is then obtained as the sum over all collisions. The energy dissipated in a collision is distributed among cosmic rays ( CR), electrons ( e) and magnetic field ( B ) such that CR + B + e = 1 and injected into the radiation model. We define the baryonic loading as f b = 1/fe = CR/ e and calculate the magnetic field following Bustamante et al. (2017) assuming equipartition between photons and the magnetic field ( B ∼ e).
Assuming acceleration of all electrons as well as fast cooling, each collision deposits Eγ,C = eEC in gamma rays. We do not explicitly model the photon fields but instead assume a fixed shape resembling observations: A broken power law peaking at = 1 keV (primed indices refer to the comoving frame of the emitting material) with spectral indices α = −1 and β = −2. We do not explicitly model the target photon spectrum here, but instead postulate that observations are described by the underlying radiation model. More explicit radiative modeling of the photon fields, such as the one in Globus et al. (2015a), would lead to a dependence of the (synchrotron) peak energy on the collision parameters, namely the magnetic field and dissipated energy per mass. For an engine as SR-0S, we expect a simple evolution from higher to lower peak energies as collisions move outside. Stochasticity in the Lorentz factor distribution is expected to add a (stochastic) spread in the distribution of peak energies. The more detailed treatment of these effects goes beyond the scope of this study and would additionally require a radiation model for the electromagnetic processes in the presence of hadronic interactions. We ensure that the maximal photon energy is high enough in order not to impact the multi-messenger production unless the optical thickness to pair-production exceeds unity, when we impose a cutoff there. For detailed modeling of GRB spectral energy distribution in the internal shock model see e.g. Bosnjak et al. (2009); Daigne et al. (2011); Bošnjak and Daigne (2014).
We simulate the nuclear system with the timedependent radiation code NeuCosmA that iteratively solves the transport equations in order to calculate the cosmicray spectra for each collision, see Boncioli et al. (2017);   Biehl et al. (2018) by balancing acceleration with losses due to photo-hadronic interactions, photo-disintegration, photopair production, synchrotron emission and adiabatic expansion of the emitting material. For the acceleration, the maximal efficiency of η = 1 is assumed (Bohm limit). Instead of the interaction timescale tint in Biehl et al. (2018), we assume the effective energy loss timescale t loss = A · t int limits the maximal energy for photo-disintegration, which is a rough estimate assuming that a single nucleon is ejected per interaction. This assumption is justified because we only consider different mass groups for the injection, namely hydrogen, helium, nitrogen, silicon and iron, which means nearby isotopes are not distinguished. Protons are limited by the effective cooling timescale for photo-meson production. The integral fractions of the injection elements are defined as They are free parameters of the simulation and are determined by the fit to UHECR data later. We define the heavy mass fraction (HMF) as ( Particle injection is assumed to persist throughout the dynamical timescale t dyn = c · l m,C . The injection luminosity is accordingly normalized to L inj = E diss /t dyn and the radiation densities are computed by distributing that luminosity/energy over an isotropic volume V = 4πR 2 C l m,C . The system is evolved over the dynamical timescale t dyn , at the end of which the spectra are extracted. Note that this assumption is slightly different from Daigne and Mochkovitch (1998); Globus et al. (2015a), where the corresponding time scale is the expansion timescale t ex .
In order to compute the emitted particle spectra, additional assumptions on the escape mechanisms have to be made. In Baerwald et al. (2013); Biehl et al. (2018) neutral particles free-stream while charged particles escape if the edge of the emitting region is within their Larmor radius. This yields an effective escape rate which is comparable to a Bohm-like diffusion process. However, the simulations presented in Globus et al. (2015a) suggest a different behavior more similar to a high-pass filter, leading harder, bell-shape escape spectra. Similarly hard escape spectra were analytically derived in Ohira et al. (2010) as ∝ exp(− ln 2 (E/Emax)), where the escape is most efficient at the maximal energy Emax. For this study, we employ the analytical form derived in Ohira et al. (2010) that yields similar results to Globus et al. (2015a) -supported by the argument that it well describes UHECR data.
A sizable fraction of collisions may occur below the photosphere defined by τ Th = 1 (where τ Th is the optical depth to Thomson scattering) for each individual shell collision, see App. A.4 in Bustamante et al. (2017) for details. Sub-photospheric collisions, in principle, may lead to high neutrino fluxes if the observed photon spectrum is simply extrapolated to below the photosphere. However, achieving efficient particle acceleration below the photo-sphere requires a fine-tuned combination of conditions (Sironi and Spitkovsky 2011;Beloborodov 2017), and the target photon spectrum almost certainly has a different shape. A simple extrapolation of the radiation model that we apply for the prompt phase therefore tends to overestimate the expected neutrino emission when assuming similar acceleration efficiency. The absence of neutrino observations from GRBs in IceCube may indicate that if GRBs are indeed sources of UHECR, the contribution of sub-photospheric collisions must not be very large in order not to violate the current limits. Therefore, there is no particle emission in our model from sub-photospheric collisions. To avoid an unphysical bias in the parameter scan due to this penalty, we exclude models with less than 40% of energy dissipated above the photosphere from the interpretation of the parameter scan.

Multi-messenger emission regions
Here, we discuss the two extreme examples SR-0S and WR-MS to illustrate the basic behavior of the source model and the effect of the source parameters. For this purpose the injection composition is fixed to identical values for the four examples, whereas in later sections the IA become free parameters of the fit. The particle emission regions in the left panel of Fig. 2   3. Since collisions occur with larger variation of Lorentz factors along different radii, the photospheric condition can be encountered within the light-shaded area. The difference between the subphotospheric extrapolation (dashed) and regular emission curves reflects that the number such collisions within the light-shaded band is small. The second bump, visible at R C ∼ 7 · 10 9 km is related to the bulk collision and controlled by the gradient of the initial Lorentz factor distribution as described in the text. While the radial dependence at the bulk collision resembles the behavior of the SR-0S case in Fig. 2, the fireball develops efficient particle emission at much larger radii.
quence of the initial Lorentz factor distribution. Collisions occur first where the gradient in speed is maximal since ∆β = βr − βs ≈ 1/2 Γ 2 s − 1/2 Γ 2 r . The ramp-up in the Lorentz factor distribution results in the bulk of the particle emission around happening at radii above RC ∼ 7 · 10 9 km. Due to absence of stochastic engine behavior in SR-0S, collisions with subsequent particle emission occur up to a radius of 5 · 10 10 km, which corresponds to the point where the latest emitted shell runs into the bulk of slow merged shells. Neutrino emission peaks close to the photosphere where col-lisions are optically thick. The radial dependence of the nuclear composition, shown in the right panel of Fig. 2, is determined by the interplay between the dominant cooling process and maximal energy and is different for each type of nucleus.
Emission of heavier nuclei starting from the CNO group is suppressed at lower radii due to photo-nuclear interactions, while efficient acceleration retains a high emission up to the maximal radius. H and He emission is highest at inter-mediate radii and drops toward outer radii since acceleration becomes less efficient.
To some extent, the SR-0S scenario is orchestrated since the initial Lorentz factor distribution is chosen such that the bulk of the contributing collisions occurs near the bulk collision at RC ∼ 7 · 10 9 km. By adding a spread AΓ > 0 to the Lorentz factor distribution a two-bump shape develops (Fig. 3) where the second bump comes from the bulk collision and the collisions at lower radii from collisions between neighboring shells due to random fluctuations in their Lorentz factors. The height of the second peak, and thus efficient UHECR emission, depends on the choice of AΓ and the amount of energy left over from the random collisions at lower radii or earlier times. If the stochasticity is excessively high, the energy remaining in the fireball for the bulk collision at larger radii is too low to produce a significant output of UHECR. At the same time, the bulk Lorentz factor decreases and harsher assumptions have to be made for Γmin and Γmax.

FITTING UHECR DATA
In this section, we include the effects of UHECR propagation and study the parameter space which describes UHECR data. We also discuss the requirements for the source energetics and the heavy mass fraction at injection from these results.

Source Population and UHECR transport models
We adopt the cosmological distribution of GRBs derived by Wanderman and Piran (2010) based on the analysis of Swift GRBs: Due to computational constraints all GRBs are computed at an injected luminosity of 10 53 erg/s instead of using the full luminosity distribution. The suppression of radiation from subphotospheric collisions yields an emitted luminosity close to 10 52.5 erg/s, which is close to the observed break luminosity in Liang et al. (2007).
For the extragalactic propagation we employ identical tools as in , i.e. the numerical code PriNCe and the distributed fitting framework. All unstable nuclear isotopes lighter than iron including neutrons decay at injection or immediately at production, since the decay length is typically much shorter than the interaction length. The chain of decay products is followed down to protons, stable nuclei and neutrinos, which are produced in charged pion and muon decays. The model for the extragalactic background light (EBL) is Gilmore et al. (2012) and the photonuclear disintegration model is Talys (Koning et al. 2007). It is noteworthy that in this study, the nuclear disintegration within the source and during propagation have been simulated with the same nuclear disintegration model.

UHECR interactions and fit method
Each GRB model (defined by a triple Γmin, Γmax and AΓ) produces an individual set of nuclear and neutrino spectra at Earth for each of the five injection masses. The superposition of the nuclear spectra yields ln A and σ(ln A) that are converted into the corresponding Xmax and σ(Xmax) following Abreu et al. (2013).
The purpose of our study is to perform a systematic parameter space scan over the model parameters and composition. Unfortunately the number of free parameters in the model (initial distribution of Lorentz factors Γ k , kinetic energies E k , width l k , and separation d k of the shells) exceeds the constraints from the (integrated and angular-averaged) spectrum and composition measurements of UHECRs. We make an initial guess for the parameters described in Section 2.1, in particular for the initial shell distribution and injections compositions, that appear suitable to fit the spectrum. The assumption of equal energy per emitted shell imposes constrains on multiple shell parameters such as separation and width. Another criterion for a suitable model is a sufficient maximal rigidity similar to that observed in the UHECR spectrum. It results in constraints on the average emission radius and magnetic fields, which are mainly controlled by the initial distribution of Lorentz factors that translate into shell speeds.
The Minuit 2 fitter (James and Roos 1975) is used for the minimization of the goodness of fit estimator The total χ 2 includes individual contribution from the allparticle spectrum and the Xmax . The data are taken from the Pierre Auger Observatory (Fenu, F. et al. 2017;Bellido, J. et al. 2017) and include a systematic energy-scale uncertainty δE = ±14%. Note that we do not include the second moment of the Xmax distribution in the total χ 2 . This is not a technical limitation and the reasoning behind this choice is explained later. In order to optimize the computation, individual simulations are performed for each injection isotope separately, to be superimposed later, i.e. in total we compute NΓ min × NΓ max × NA Γ × 5 injection masses = 12 × 10 × 11 × 5 = 6600 fireball evolutions. We minimize the five injection fractions IA and the energy shift δE of the Pierre Auger Observatory data. Note that while in the propagation computation the relation between the injected and ejected mass fraction is linear, this is only approximately true for the source model since a higher baryonic fraction may shift some collisions into the photosphere and change the emission spectra, see discussion in Bustamante et al. (2017). We verified in separate simulations that the impact of this non-linearity is small within our current parameter ranges.
The confidence contours are drawn using ∆χ 2 = χ 2 − χ 2 min projected onto planes of two parameters by minimizing over all other fit parameters. The best fit point is found by minimizing over all points in the (Γmin, Γmax, AΓ) cube.

Systematic parameter space search
As outlined in Section 2.3 the engine parameters Γmin, Γmax and AΓ affect the distribution of collisions along the jet and their properties, and thus impact the particle interactions and shift the emission regions of messengers. We systematically scan the engine parameter space in Γmin, Γmax and AΓ assuming that the model represents the emission of all GRB that power the UHECR flux. We continuously adjust the fraction that goes into non-thermal baryons f b and the integral injection fractions for the nuclear species such that the observed UHECR spectrum and composition is fitted for each engine configuration. As we will see, the fit contours enclose a relatively wide range of GRB realizations, implying that it is possible to obtain similar fits but with a superposition of multiple GRB models. We do not attempt this kind of generalization since it simply increases the number of free parameters.
The fit result is shown in Fig. 4 where the third variable is integrated out or summed over for each two-dimensional distribution. For each point in Γmin, Γmax and AΓ, the five baryonic loadings, the energy scale uncertainty of the data has been obtained by (continuous) minimization. The contours for the 3σ confidence interval enclose a large parameter space, demonstrating that the model (fit) is relatively robust against parameter changes. Engines without stochastic behavior AΓ = 0 are disfavored as well as those with Γmax < 200. The 1-and 2σ contours favor engines that produce an average Γ bulk ∼ 200 − 400 indicated by the isocontours in the upper left panel. That is per se interesting, as a bulk Lorentz factor around 300 has been reconstructed from observations (see, for example, Ghirlanda et al. (2018)), which we recover by completely different means. The best fit belongs to the WR-MS engine type with a weak ramp up and medium stochasticity. Our model is not applicable within the white shaded region at lower Γmin in which more than 60% of energy is dissipated below the photosphere. Except a small overlap with the 1-σ contour and a narrower secondary minimum with a very high stochasticity, the UHECR fit prefers GRB realizations for which most energy is dissipated above the photosphere.
For the discussions that follow further below, we choose four examples indicated by the markers in Fig. 4. The rough criteria to choose these examples we motivated by a similar Γ bulk ∼ 320 as the best fit (WR-MS), very different stochasticity (SR-0S vs WR-HS) and one example with a low AΓ (SR-LS). All examples lie within (or close to) the three sigma contours.
The details on the fit result are summarized in Tab. 1, which we will also come back to later. The baryonic loadings are comparable to the simplistic expectation if equal number rates of electrons and protons are picked up at similar (low energies) by the acceleration process, for which one obtains f b (mp/me) 1/2 44, see e.g. Pohl (1993). The dissipation efficiency is around 20%, which is within the typical range expected for GRB internal shock scenarios, see Rudolph et al. (2019) for a more detailed discussion.

Source spectra
To illustrate how different engine behaviors affect the ejected cosmic-ray spectra, we show in Fig. 5 the ejected spectrum from the source for each of the four example cases. The weights for the five intregral fractions of the injection masses are obtained from the fit, the initial composition of the jet is thus different for each of the example cases.
By adding stochasticity the emission of nuclei is spread out among a larger range of radii compared to the case of no stochasticity (SR-0S). With stochasticity, collisions can occur at smaller radii, where optical depths for photonuclear interactions are high and more secondary nucleons (red curves) are produced through disintegration. This also leads to higher neutrino emission.
The smaller cosmic-ray production region close to the engine of SR-0S is also reflected in smaller maximal energies of cosmic-ray nuclei, as the upper left panel of Fig. 5 clearly shows. Higher maximal energies are reached by stochastic engines through collisions that occur at outer radii where the acceleration of heavier nuclei is still efficient but the source is transparent enough for nuclei to escape. These collisions will negligibly contribute to the total neutrino output. The distribution of collisions among many radii in stochastic models can, however, yield neutrino bright and UHECR dim, and, UHECR bright and neutrino dim collisions within the same astrophysical object in contrast to one-zone models as in Biehl et al. (2018) or multi-collision models dominated by one bulk collision (SR-0S).
Concerning the spectral shape, the contributions of collisions that are spread out along all radii in models with higher stochastic component lead to a softening/broadening of the emission spectra for individual mass groups. This is demonstrated in Fig. 5 where the spectra from each collision are illustrated by thin and the total output by thick curves. For models dominated by the bulk collision (SR-0S and SR-LS) the spectra are clearly narrower.

UHECR spectra at Earth and general aspects of the fit
Following the ejected spectra we discuss the corresponding observed spectra and general properties of the fit for the four example cases. Comparing observed to ejected spectra, one has to bear in mind general effects of extragalactic propagation. Interactions with background photons during propagation generally lead to a flux depletion of nuclei with energies above the threshold for photo-disintegration. Secondary nuclei that emerge from photo-nuclear interactions populate the spectra at lower energies, softening the spectrum at Earth compared to that at the source. Because the boost is approximately conserved during disintegration, the energy Ep of the interacting nucleus of mass Ap and a secondary nucleus of mass As is fixed by the relation Es = As/Ap · Ep. The lighter secondary nuclei from disintegration, therefore, show up at lower energy.
The spectra at Earth are shown in Fig. 6. It is remarkable how well the spectrum and the Xmax are reproduced with such different source spectra and engine properties.
To fit the observed UHECR spectrum for different engine realisations, we adjust the integral fractions of the injection elements at the source (Eq. 6) for the four benchmark cases. The integral fractions resulting from the fit are listed in Table 1. In all cases data require more than 50% of intermediate nuclei (N and Si) injected at the source. The . Parameter space in Γmax, Γmax and A Γ for the fit to the UHECR spectrum and composition observable Xmax . The orange, green and blue contours are given for ∆χ 2 = 1, 4 and 9, respectively, corresponding to 1, 2 and 3σ for 1 d.o.f. for a Gaussian likelihood; see also color scale. For each two dimensional panel, the ∆χ 2 is minimized over the third (unshown) parameter and the injected composition.
In the gray shade area, less than 40% of the energy is dissipated in super-photospheric collisions. The black contours in the upper left panel correspond to the bulk Lorentz factor of the initial shell setup. Our benchmark scenarios are marked as in the figure legend. While the actual best-fit is in the left orange contour in the upper left panel, we use WR-MS as similar scenario which has a ∆χ 2 very close-by (and a slightly higher bulk Lorentz factor which is comparable to the other benchmarks).  WR-MS and WR-HS cases require an even higher (> 70%) integral fraction of silicon and iron at injection to compensate for their depletion due to efficient photo-disintegration in the source. An unavoidable counterpart to efficient photodisintegration is the high abundance of secondary neutrons in the source (thin red curves in Fig. 5). This is reflected in the integral fractions that do not require any additional proton injection in the source when stochasticity is present.
In the absence of stochasticity the SR-0S case requires a 22% proton integral fraction at injection to compensate for the low disintegration in the source. Since the high maximal energies required by data are barely reached at the source, the fit pulls the energy scale (δE < 0 means data move toward lower energy) to minimize the tension. Because maximal energies are smaller, the missing abundance of light elements is not well compensated through disintegration during propagation. The high maximal energies at the source for WR-MS and WR-HS (visible in the lower panels of Fig. 5) affect the spectrum at and beyond the cutoff energies in the lower panels of Fig. 6 that show an onset of recovery above 2 · 10 11 GeV.
Interestingly, the SR-0S case does not require any (primary) helium to be injected into the source. Therefore, all helium is secondary, i.e. a product of photo-disintegration in the source or the propagation. However, we notice that this might also be an effect of the disintegration model (Talys) used in the computation. The production of helium, and its subsequent disintegration, is strongly affected by the uncertainties in the disintegration cascade, due to absence of data and models (Alves Batista et al. 2015;Boncioli et al. 2017). For the stochastic cases the origin of helium in UHECR at Earth is both primary and secondary.
In simpler models, the succession of increasingly heavier mass spectra toward higher energies is often assumed to be caused by a maximal rigidity reached by the accelerator or one acceleration zone. As we show with these stochastic multi-collision models, this assumption is not essential to describe the spectrum and the Xmax . In SR-LS and WR-MS the heaviest mass group at the cutoff is mostly silicon and not iron. For WR-MS and WR-HS the cutoff for the proton spectrum at Earth reaches or exceeds that of the helium or the nitrogen group, confirming that the data do not require that the maximal energy at ejection follows the Peters cycle, as already found for example in Biehl et al. (2018). This is already visible in the source spectra and is enhanced at Earth since protons are abundantly produced during propagation of heavier elements with energies close to, or above the observed cutoff.
A real discrimination of the models and a crucial piece of the origin of UHECR puzzle may come from the width of the Xmax distribution, σ(Xmax). The data requires that while the composition has to become heavier toward higher energies, this transition has to happen through a smooth gradual succession between neighboring mass groups (see e.g. Aab H He N Fe Figure 6. Observed UHECR spectrum (large panels) for two moments of Xmax (small panels) for the four example cases for the best-fit composition (see Tab. 1). Only the spectrum and mean Xmax are included in the fitting procedure. The gray shaded area indicates the range below 6 · 10 10 GeV which is excluded from the fit.
et al. (2017)) ,which minimizes the overlap of different mass spectra. For example, the Xmax could in principle be described by a mixture of just protons and iron. The overlap of protons and heavier spectra would, however, lead to large fluctuations in Xmax producing σ(Xmax) distribution that is flat or increasing in energy, contrary to what is observed. Clearly, our model does not describe σ(Xmax) sufficiently well. If we include σ(Xmax) in the χ 2 definition in Eq. (8), the result only marginally improves at the cost of narrow contours. The preferred models tend to have fewer optically thin collisions contributing to the emission and mostly lie at the edge of the range in which our model is applicable. The reasons for this behavior are twofold and can be understood by comparing Figures 5 and 6: (i) Since collisions from different radii contribute to the total output, the maximal energies of individual collisions cover a wider range (see the colored "bands" in the SR-0S panel of Fig. 5 that are comprised of successive thin curves). When those are summed (thick, colored curves for each mass group) the total spectra are much wider than the individual bell shapes, equivalent to a softening of the mass group spectra. For stochastic engines this effect is stronger, since there is a wider spread in the maximal energy of individual positions. In combination with a positive source evolution, softer spectra are known to produce worse fits.  showed that this behavior persists for positive redshift evolutions across various model combinations used for the simulation of air showers and extragalactic propagation.
(ii) The second reason is the strong overlap of light and heavy mass spectra in stochastic models, in particular the abundance of protons at the highest energies in the lower panels of Fig. 6. The required smooth succession of heavier masses is violated and hence σ(Xmax) can not decrease as in the data. This is related to the level of disintegration in the source, the maximal energy required to fit the observed spectrum, and the level of disintegration of heavy elements on the extragalactic photon fields during propagation. A significant suppression of in-source disintegration is not easily possible if electrons and nuclei are accelerated by the same mechanism and radiate within the same volume.
We note that similar effects will occur when giving up the assumption of identical GRBs. If multiple GRBs or other generic accelerators with different maximal rigidities or intrinsic luminosities contribute to the observed average UHECR flux, the average of source spectra will be softer/broader than that of individual objects. The same problem in describing σ(Xmax) will therefore occur for any model that uses a population of distant sources with different maximal rigidities. On the other hand, a mechanism that generally suppresses photo-disintegration may reduce the tension of the model with σ(Xmax) observations, but would likely also reduce emission of secondary messengers such as photons or neutrinos, rendering the attempt to discover UHECR sources with multimessenger techniques a stiffer challenge than it already is.

Source energetics
A well known problem for the GRB origin of UHECRs is the large required isotropic-equivalent energy emitted in baryons around 2 − 3 · 10 53 erg per GRB in the UHECR range (Baerwald et al. 2015) (here computed for the Wanderman-Piran GRB evolution). This estimate is a consequence of the required local emissivity ∼ 10 44 erg Mpc −3 yr −1 to sustain the flux of UHECRs with a local GRB rate around 1 Gpc −3 yr −1 . This implies that the required kinetic energy must be much larger, even if most of it is efficiently transferred into escaping non-thermal baryons with the highest energies. The dissipation efficiency (transfer from kinetic energy to nonthermal radiation) is only one part of the problem, the other is how to accelerate baryons to the highest energies without retaining much of it below the UHECR energy range.
The four examples have been chosen to describe GRBs with gamma-ray energies in the ball park of a few times 10 52 erg (row Eγ in Tab. 1). In fact, the overall normalization is Eγ ≡ 10 53 erg but sub-photospheric collisions do not contribute to this budget and hence the resulting energy output is lower. The fraction of energy emitted in super-photospheric collisions (see respective row) is fsup ≡ Eγ/10 53 erg ranging between 40% and 80% for the majority of parameter combinations.
The initial kinetic energy E kin,init is mostly converted into baryons due to relatively high baryonic loading. The energy ejected as UHECRs can be written as 3 where diss is the fraction of kinetic energy dissipated into non-thermal radiation according to the collision model, 3 There is also a small correction factor taking into account that not all injected energy (counted in f b ) is available for E src CR because of radiation losses, which we neglect in these considerations. f bol ≡ E src UHECR /E src CR is a bolometric correction describing how much energy is deposited into the UHECR range, and fesc ≡ E esc UHECR /E src UHECR is the fraction of energy escaping the source. From the table, we find diss 0.13 − 0.28, fesc 0.5, and f bol 0.07 − 0.08 for the chosen four examples. Apart from the known dissipation efficiency problem a large correction comes from f bol that describes the fraction of non-thermal baryons in the UHECR energy range. Although our model implies that the baryons are picked up at low energies, such as in a thermal bath, this factor would be an order of magnitude larger if all non-thermal particles were accelerated up to UHECR energies or if the acceleration spectra were harder than E −2 .
We find that the (isotropic-equivalent) kinetic energy of the outflow is ∼ 10 55 − 10 56 erg per GRB (see E kin,init row in Tab. 1) if the long-duration GRBs are the sole sources of UHECRs. Recent afterglow observations indicate that a substantial fraction of the bulk kinetic energy can escape observation (Acciari et al. 2019;Abdalla et al. 2019), supporting earlier arguments (Fan and Piran 2006;Beniamini et al. 2015Beniamini et al. , 2016) about a systematic underestimation of kinetic energy. But even in view of the current results our findings may sound excessive. Some ingredients of this estimate are certainly model-dependent. For example, an order of magnitude lower values ∼ 10 54 − 10 55 erg are found for hard (significantly harder than E −2 ) acceleration spectra, or when assuming an acceleration mechanism with an extremely high transfer efficiency to the highest energies. Nevertheless, from the energetics point of view a paradigm shift is required if conventional GRBs are indeed the sources of UHECRs. Neutrino observations are an opportunity to independently test this scenario.
It is interesting to compare our result to Gottlieb et al. (2020), who find that "all intermittent jets are subject to heavy baryon contamination that inhibits the emission at and above the photosphere". As we shown this baryon contamination is needed to describe UHECRs. Our radiative (not dissipation) efficiency (energy in photons versus kinetic energy) is between 5 · 10 −4 and 0.003 (see Tab. 1), which is comparable to their result. So perhaps their negative conclusion for the gamma-ray signal actually indicates that GRBs with intermittent engines are indeed efficient cosmic ray accelerators -for which a large baryon contamination is required. In contrast to the choice of generic intermittance patters in Gottlieb et al. (2020), our empirical Lorentz factor profile and its parameters tend to suppress early subphotospheric collisions and avoid the rapid slow down of the wind. Within our framework, we can not answer if such patterns can be realized in "ab initio" hydrodynamical simulations.

Heavy mass fraction at injection
The heavy mass fraction (HMF), defined as the fraction of integrated energy released into isotopes heavier than He, is reported in the last row of Tab.1, corresponding to the example cases and the best fit. The increase in the stochasticity clearly enhances the efficiency of the interactions in the engine, imposing a higher injection of species heavier than He in order to fit the measured composition. The WR-HS, being extreme in terms of stochasticity, is less efficient in providing enough energy at larger radii, and the requested HMF are smaller. The HMF is above 0.7 within the 95% CL (see Fig.7, left panel), and above 0.4 at the 3σ CL. This is mostly due to the increase in the proton fraction (see Fig.7, right panel) which is allowed within the three sigma contour, that includes a smaller level of stochasticity with respect to the one sigma contour. The large value of the HMF is an outcome of the fit, compared to Globus et al. (2015a), where it is ten times larger than found in Galactic cosmic rays.
The separation of heavy and light masses at injection becomes milder when the σ(Xmax) is included in the fit procedure. For this reason, the allowed HMF is larger than 0.6 at 3σ if the fit includes the σ(Xmax), and the spread of the mass fractions is much less pronounced.
In Zhang et al. (2018) several pre-supernova models have been selected from Woosley and Heger (2006), to test the GRBs as sources of UHECRs. In particular, models which provide a heavy distribution of nuclear mass fractions at the onset of the core collapse are chosen in this paper, and such compositions are used as output of low-luminosity GRBs with internal shock model, that are supposed to power the UHECR flux at Earth. These have HMF greater than 0.9, which is found to be compatible with our result independent of the confidence level.

LIGHT CURVES AND MULTI-MESSENGER IMPLICATIONS
Since the UHECRs alone cannot discriminate among different model assumptions, we seek here for additional observables and multi-messenger signals to study the underlying model, in particular, light curves and the neutrino flux.

Predicted light curves
In Fig. 8, we show the light curves, both in gamma rays and neutrinos, for an individual GRB for the four different model assumptions Fig. 1. Note that these examples are sorted by increasing stochasticity and decreasing engine ramp-up. The case SR-0S corresponds to a very disciplined engine in Bustamante et al. (2017), leading to a single-pulsed light curve in gamma rays. The increasing engine stochasticity adds time variability to that light curve, see SR-LS and WR-MS. In the extreme case WR-HS, the ramp-up is sub-dominant, the stochasticity of the engine leads to a very spiky light curve, and the underlying (longer) pulse structure is gone. Thus these gamma-ray light curves are clearly qualitatively different, and may be used as model discriminators for the individual GRB. Note that the height of the individual peaks is determined by the luminosity of each collision, which increases with the difference in the Lorentz factors of the colliding shells.
If the light curves are to be used as model discriminator for the UHECR model, these results need to be compared to the whole population of GRBs. That, however, is not trivial, as GRBs come in a large variety of light curves and there are limitations to resolve the time variability especially for GRBs detected with low statistics; thus it is likely that there are selection effects. For example, while most detected GRBs appear to be single-pulsed with very simple structure, they may be detected at too low statistics to resolve any features or time structure. It is therefore conceivable that the whole GRB population consists of a mix of different light curve types, similar to the ones in Fig. 8. That is a subject beyond the scope of this work which requires further study if light curves are to be used as model-discriminator.
The neutrino light curves are correlated with the gamma-rays, although there is no one-to-one correspondence. The reason for that is that the pion production efficiency scales ∝ R −2 C and that the innermost collisions lead to higher neutrino production. In stochastic models (e.g. WR-HS) the collision radii are randomly distributed and, in particular, not correlated with observation time. The neutrino peaks are randomly enhanced compared to the gamma ray peaks depending on where that collision occurred. In deterministic models (e.g. SR-0S), collision radius and observation time are correlated such that the first collisions come from the innermost radii and that the neutrino production efficiency quickly drops with time. In such GRBs, the veryhigh energetic (TeV) gamma rays can be surpressed early-on because of gamma-gamma pair production, see Bustamante et al. (2017) for a more detailed discussion.  Fig. 1  (see figure titles). The light curves are obtained by assuming that each collision emits a fast rise and exponential decay peak ('FRED'), which is normalized to its total energy ouput. Each light curve is shown for a single GRB assuming a redshift of two.
One could use these observations for the optimization of the GRB stacking searches: For single-pulsed light curves, the neutrinos are expected to arrive early after the gammaray trigger (in the first few seconds for the s), whereas for highly-variable light curves, the neutrinos could arrive at any given time during T90.

Post-dicted neutrino fluxes from GRBs
Neutrinos have been proposed in the past as a model discriminator to potentially rule out the UHECR origin from GRBs. We therefore derive the "post-dicted" prompt and cosmogenic neutrino fluxes from the 3σ contour of the UHECR fit in Fig. 4. They are shown for the source/prompt neutrinos and cosmogenic neutrinos as shaded regions in the left and right panels of Fig. 9, respectively. The SR-0S example is outside of the shaded range, since it is not within the 3σ contours. The deterministic engine (SR-0S) produces the lowest neutrino flux since the typical collision radius is large. The neutrino flux increases with stochasticity, as there are more collisions at low collision radii (but above the photosphere), see Fig. 2 and Fig. 3.
The post-dicted prompt neutrino fluxes (left panel of Fig. 9) are well below the GRB stacking limits from the Ice-Cube Observatory (Aartsen et al. 2017). One way to interpret this result is that even with generous variations among individual GRBs, the absence of neutrino associations in the present detectors does not exclude GRBs as the origin of UHECRs. Our model does not impose exotic fireball parameters, a high baryonic loading or excessive Lorentz factors, and that this result is strictly derived in a post-dictive unbiased way from the UHECR fit. This means that there is no a priori bias from the neutrinos included, and it is somewhat surprising that no model within the 3σ UHECR fit produces a prompt neutrino flux detectable in IceCube. The fluxes are compatible with earlier estimates in Bustamante et al. (2014Bustamante et al. ( , 2017; Rudolph et al. (2019) for a fixed baryonic loading. The baryonic loadings are derived from the UHECR fit (see Tab. 1) by integrating over the entire source population. This result demonstrates that the initial interpretation of IceCube's non-observation (Abbasi et al. 2012) can be regarded as too strong, and that the UHECR origin from GRBs cannot be ruled out based on neutrino observations, yet.
The next generation detectors, such as KM3Net-ARCA in the Mediterranean Sea (Aiello et al. 2019), demonstrate promising full sky sensitivity estimates for point-source detection and hence may detect some of the GRBs. The planned IceCube-Gen2 detector at the South Pole will also have an enhanced sensitivity to GRBs because of the significantly larger effective area and because the stacking search should be still statistics-limited. It is therefore conceivable that the next-generation of experiments can finally constrain the UHECR paradigm with neutrinos by at least about a factor of ten larger exposure -which is the curve we show 10 4 10 5 10 6 10 7 10 8 10 9 10 10 E [GeV] in Fig. 9 for comparison. It is expected that these future detectors can at least exclude the cases with high source stochasticity. Concerning cosmogenic neutrinos, Heinze et al. (2019) and Alves  demonstrated that for alike, homogeneously distributed UHECR sources that accelerate nuclei up to a maximal rigidity, the detection of cosmogenic neutrinos from UHECR nuclei is out of reach for the next generation detectors. This statement is valid considering model systematics of the propagation and the air-shower model ). The present model captures some of the variety of observed light curves by scanning over engine properties. The neutrino post-dictions in right panel of Fig. 9, therefore, include some non-trivial scenarios with (for example) a high-energy proton component at energies higher than in simple rigidity-dependent sources models. This sub-leading proton contribution increases with the level of stochasticity in our model, see Fig. 6. Since these protons reach the threshold for CMB interactions, the cosmogenic neutrino flux is significantly enhanced, see van Vliet et al.
(2019) for a detailed discussion. This is prominently visible in the right panel of Fig. 9 that indicates the possibility to observe a diffuse component from GRB with the next generation radio detectors.
We notice here that the use of the σ(Xmax) in the fit would reduce those cases in which protons/neutrons reach very high energies in the escape spectra, with the consequence of suppressing the production of neutrinos in the extragalactic space. For this reason, the prediction for the cosmogenic neutrino within 3σ of the UHECR fit would be below the GRAND sensitivity. Since the production of prompt neutrinos is mostly dependent on the efficiency of the in-source interactions and not on the maximum energy of the cosmic rays, the result on the source neutrino flux is qualitatively not affected by the use of the σ(Xmax) in the fit.

SUMMARY AND CONCLUSIONS
We have performed a systematic parameter space study of the engine properties of GRBs in the internal shock scenario in the context of multi-messenger observations. Since neutrino observations are known to constrain the UHECR origin from GRBs in one-zone models, our model includes multiple internal shocks, and can describe both stochastic engines and deterministic ramp-ups over a wide parameter space. The main target of our study has been the question if long-duration GRBs can describe UHECR spectrum and composition data in spite of existing limits from the nonobservation of GRB (prompt and cosmogenic) neutrinos.
We have demonstrated that UHECR data (spectrum and Xmax ) can indeed be described in a wide range of engine parameters, and that UHECR data alone cannot discriminate among different options. The observed σ(Xmax), however, indicates a rather pure composition at the highest energies, which is intrinsically difficult to obtain in multiregion or population models in which a diversity of production regions or sources contribute. Here the superposition of UHECR emissions from different collision radii results in a larget overlap of light and heavy masses than observed in data for the majority of the tested engine parameters. We speculate that either σ(Xmax) is more strongly affected by hadronic uncertainties than anticipated, or that significant fine-tuning is needed to describe σ(Xmax) if GRBs are to describe the diffuse UHECR flux -otherwise GRBs may not be the dominant source of the UHECR flux. A similar argument can be made for other source populations that are dominated by higher redshifts, unless there is a mechanism to fine-tune the maximal rigidity of each source. Therefore the trend in σ(Xmax) might be better described by a small population of local sources.
We have "post-dicted" the neutrino flux from the 3σ allowed region of the UHECR fit, and found that the expected neutrino flux is well below the current stacking bound in consistency with earlier estimates from multi-collision models. This means that the UHECR paradigm cannot yet be ruled out with neutrino data alone; we however expect that the next generation of neutrino telescopes could observe GRB neutrinos if GRBs powered the UHECR flux. If no neutrinos are observed, especially stochastic engine models will be effectively constrained with these detectors.
We have also shown that the fraction of nuclei heavier than helium to be injected at the source has to be larger than 70% (at the 95% CL) in order to describe the UHECR data. This reflects the distribution of isotopes at the end of the core-collapse, which can be attributed to several different characteristics of the pre-supernova models.
As possible model discriminators, which are sensitive to the stochasticity of the engine, the GRB light curves have been identified. While these can be easily used for individual GRBs, it may be difficult to associate them with a whole population of GRBs.
Furthermore, we have demonstrated that the (isotropicequivalent) kinetic energy of the outflow has to be ∼ 10 55 erg − 10 56 erg per GRB implying a low transfer efficiency into electromagnetic radiation, if the population of long-duration GRBs is powering the UHECRs. The dissipation efficiency problem (transfer from kinetic energy to nonthermal radiation including baryons) is in that case only one part of the problem, the other is how to accelerate baryons to the highest energies without leaving too much energy below the UHECR energy. The energy requirement can be somewhat relaxed for hard (significantly harder than E −2 ) acceleration spectra, or an acceleration mechanism with an extremely high transfer efficiency to the highest energies. Nevertheless, this extreme energy requirement is difficult to circumvent, and may be tested by afterglow observations. We conclude that the GRB paradigm for UHECRs cannot be uniquely excluded at this point. On the one hand, the description of UHECR data has revealed that spectrum and Xmax can in fact described for a wide range of engine parameters with an inferred bulk Lorentz factor 200 − 400 within expectations, and the post-dicted prompt and cosmogenic neutrino fluxes are below current limits without any prior or bias. On the other hand, major obstacles have been the description of the observed σ(Xmax) and the required energetics. Only independent observations, such as limits or detections from the next generation of neutrino detectors, will be able to robustly exclude the GRB-UHECR connection.