The cosmic-ray content of the Orion-Eridanus superbubble

The nearby Orion-Eridanus superbubble, which was blown by multiple supernovae several Myr ago, has likely produced cosmic rays. Its turbulent medium, still energised by massive stars, can impact cosmic-ray transport locally. The gamma rays produced in cosmic-ray interactions with interstellar gas were used to compare the GeV to TeV cosmic-ray spectrum in the superbubble and in other regions near the Sun. We used ten years of Fermi-LAT data in the 0.25-63 GeV energy range to study the closer (Eridanus) end of the superbubble. We modelled the spatial and spectral distributions of the gamma rays produced in the different gas phases of the clouds found in this direction. We found that the gamma-ray emissivity spectrum of the gas along the outer rim and in a shell inside the superbubble is consistent with the average spectrum measured in the solar neighbourhood. This result calls for a detailed assessment of the recent supernova rate and census of massive stellar winds in the superbubble in order to estimate the epoch and rate of cosmic-ray production and to constrain the transport conditions that can lead to such homogeneity and little re-acceleration. We also found significant evidence that a diffuse cloud lying outside the superbubble, at a height of 200-250 pc below the Galactic plane, is pervaded by a 34\% lower cosmic-ray flux, but with the same particle energy distribution as the local one. Super-GeV cosmic rays should freely cross such a diffuse atomic cloud without significant loss or spectral distorsion. We tentatively propose that the cosmic-ray loss relates to the orientation of the magnetic field lines threading the cirrus, which point towards the halo according to the dust polarisation data. We gathered past and present emissivity measurements near the Sun to show how the local cosmic-ray flux decreases with Galactic height and to compare this trend with model predictions.


Introduction
Superbubbles form around starburst regions under the combined and sustained influence of the ionising radiation and energetic winds of massive stars plus a series of supernovae (SNe).They have long been considered as potential sites for cosmic-ray (CR) production because of the collective action of the multiple shock waves they contain (Montmerle 1979;Cesarsky & Montmerle 1983;Bykov & Toptygin 2001;Bykov 2001;Parizot et al. 2004;Ferrand & Marcowith 2010).The abundances of CR nuclei also require a mix of massive-star outflows, supernova ejecta, and interstellar matter that is easily found in superbubbles (Lingenfelter 2018, and references therein).The relative fractions in the mix depend on the uncertain stellar yields for heavy nuclei (Binns et al. 2019).Superbubbles can also alter the distribution of Galactic CR diffusing through them because of the broad spectrum of magnetic turbulence powered by their multi-scale flows.
Modelling the acceleration, transport, and losses of CR inside the turbulent and multiphasic medium of a superbubble is particularly complex (see the review by Bykov 2014).Intermittent acceleration by individual shocks, followed by interactions with large-scale compression and rarefaction waves, can produce hard CR spectra in the MeV-TeV range (Bykov 2001).The momentum distribution of the particles evolves over a few million years (Myr).It asymptotically tends toward a power law for p 2 f (p) ∝ p −γ , with a γ index ranging between two and three depending on the particle injection momentum: γ = 3 for the acceleration of suprathermal particles and γ = 2 for the reacceleration of relativistic CR.Such superbubbles can transfer 10% or more of their kinetic power to CR (Bykov 2001).
Superbubbles can also affect Galactic CR entering them.Particle interactions with magnetohydrodynamic (MHD) waves in the turbulent plasma can re-accelerate the CR despite pion losses along their path.For efficient superbubbles, with acceleration time scales much shorter than the diffusive escape time scale from the bubble, the resulting hadronic γ-ray emission can exhibit power-law differential spectra as hard as E −0.3 γ in the 1-100 GeV energy band (Tolksdorf et al. 2019).The re-accelerated CR should accumulate at the periphery of the superbubble as the enhanced turbulence prevents their penetration deep into the interior (Tolksdorf et al. 2019).Conversely, CR that are produced internally can remain confined in the superbubble medium for typically 300 kyr for 100 GeV particles in a 25-50 pc superbubble if the diffusion lengths are 100 times shorter than in the standard interstellar medium (ISM; Ackermann et al. 2011).Superbubbles have therefore enough energy and magnetic turbulence to substantially modify the CR properties, but whether A&A 635, A96 (2020) they have a positive (acceleration) or negative (confinement and losses) feedback is still an open question that needs clarifications from observations.Ackermann et al. (2011) reported an extended excess of hard γ-ray emission from the Cygnus X starburst region, which they interpreted as a cocoon of freshly accelerated CR inside the 60 pc-wide superbubble.The latter has been carved by the stellar winds and ionisation fronts from young stellar clusters including Cyg OB2 and NGC 6910.Cyg OB2 is one of the most massive clusters in the Galaxy at a distance of 1755 +23  −19 (stat) +373 −261 (syst) kpc (Berlanas et al. 2019) and the NGC 6910 cluster is part of the Cyg OB9 association at a distance of 1600 ± 200 pc (Kolaczkowski et al. 2004).The E −2.1 γ differential γ-ray spectrum extends to 100 GeV in the Fermi-LAT data and softens to E −2.6 γ in the 0.4-10 TeV range in the MILAGRO, ARGO-YBJ, and HAWC data (Bartoli et al. 2014;Hona et al. 2017).The apparent break energy around 1 TeV is uncertain because of the confusion with other extended TeV sources in these crowded directions.The γ-ray luminosity does not exceed a few percent of the wind power in the stellar clusters.Star formation has progressed across the superbubble from the low longitude end 10 Myr ago towards the younger Cyg OB2 cluster (Berlanas et al. 2018).The latter has formed stars more or less continuously between 1 and 7 Myr ago and its mass distribution suggests that its most massive members have already evolved to their supernova end state (Wright et al. 2015).Five or six of them may have exploded within the last 1-2 Myr (Lingenfelter 2018).NGC 6910 has a comparable age of 6 ± 2 Myr (Kolaczkowski et al. 2004), but is also part of the older and broader Cyg OB9 association.
A second possible cocoon of young CR has been reported in γ rays further away in the Galactic disc, around l = 25 • and b = 0 • , in a gas cavity found near a candidate massive OB association, G25.18+0.26,seen in X rays (Katsuta et al. 2017).The distance estimate of 7.7 kpc implies a three times larger bubble (210 × 170 pc) and a 14 times larger γ-ray luminosity than for the Cygnus X cocoon, but with a similarly hard emission spectrum extending to 100 GeV.The large distance and heavy extinction, however, hamper further optical analyses of the stellar cluster.Confusion with other γ-ray sources and with the bright Galactic emission further hampers γ-ray analyses of high-energy CR until better resolved observations with the Cherenkov telescope array (CTA) become available (Acharya et al. 2013).Both the Cygnus X and G25.0+0.0 cocoons unfortunately lie in confused directions respectively tangent to the Local spiral arm and Scutum-Centaurus arm.
The Orion-Eridanus superbubble is the nearest example of a superbubble.This over-pressured and elongated cavity, about 200 pc in width and 250 pc in length, slowly expands at a velocity of 20 km s −1 (Joubaud et al. 2019, hereafter Paper I).Its close end, in the Eridanus constellation, approaches the Local Bubble at a distance of 150-200 pc from the Sun.Its far end lies slightly beyond the 388 ± 5 pc distant, 1-2 Myr old stellar clusters of the Orion Nebula (Pon et al. 2016;Kounkel et al. 2017;Zari et al. 2017).It has been formed by the ionising radiation and energetic winds of tens of massive stars and by a series of 10-20 supernovae (Bally 2008) that have occurred at a rate of about 1 Myr −1 over the past 12 Myr (Voss et al. 2010).The composite structure of the superbubble has likely evolved in space and time from near to far, along the blue stream of massive stars identified by Pellizza et al. (2005) and Bouy & Alves (2015) in front of the Orion clouds.The stream extends over a length of 150 pc along the major axis of the superbubble.The stellar ages span from 20 Myr (near) down to 1 Myr (far).The near end of the cavity is filled with a mix of hot and warm plasmas with temperatures of order (3-9) and (0.3−1.2) MK, and gas densities of order 0.005 and 0.05 cm −3 , respectively, Paper I. It is bounded by shells of neutral gas that have been swept-up and compressed by the expanding outer shock wave.Several shells have been isolated in position, velocity, and distance (Paper I) and they were used here to probe the CR flux in various parts of the superbubble.They are sketched against the H α intensity map of the superbubble in Fig. 1 and their H I column density and CO intensity maps are presented in Fig. 2 of Joubaud et al. (2019).
We studied the near part of the superbubble using γ-ray data obtained with the Fermi-LAT (Atwood et al. 2009).This part lies at large (negative) Galactic latitudes, from −50 to −15 • , away from Galactic confusion.To do so, we modelled the γray emission coming from the hadronic interactions of CR with the different gas phases and clouds of the superbubble.We coupled γ-ray and dust analyses to trace the total gas, in particular the dark neutral medium that is not detected via H I and CO observations.Our analysis yielded γ -ray emissivities per gas nucleon in the atomic phase of the different clouds, where the H I line data allow a direct estimation of the gas mass pervaded by the CR.We could thus compare the CR flux in the superbubble with the average measured in the local ISM and in the Solar System (Casandjian 2015).Gamma-ray emissivity estimates in the local ISM have indeed shown that the CR flux is quite uniform within a few hundred parsecs from the Sun and that it is in agreement with the direct measurements obtained in the Solar System (see the review by Grenier et al. 2015).Other works, based on different interaction cross-sections (Strong & Fermi-LAT Collaboration 2015;Orlando 2018), found a ∼30% discrepancy in the emissivities between direct measurements in the Solar System and the local ISM average.Residual cloud to cloud variations of order 10% are commensurate with systematic uncertainties in the atomic gas mass that has served for the measurements.
The γ-ray emissivities per gas nucleon have also been shown to decrease with increasing altitude above the Galactic plane, by about 50% within 2 kpc (Tibaldo et al. 2015).Our analyses allowed the study of the CR flux in a diffuse cirrus cloud lying in the vicinity of the superbubble, below it with respect to the Galactic plane.
The paper is structured as follows : data are presented in Sect.2, we describe the γ-ray and dust models in Sect.3. The results are presented in Sect. 4 for the superbubble first and then the cirrus.The results are discussed in Sect. 5 in the same order.

Data
The present analysis is based on the same gas data and on the same atomic and molecular cloud separation as in the study of the dynamics and gas content of the superbubble presented in Paper I. We used the same analysis region towards the Eridanus part of the superbubble, extending in equatorial coordinates from 43 • to 78 • in right ascension and from −29 • to 21 • in declination, as shown in Fig. 1.We masked out two 5 • -wide areas on the western and eastern sides of the region to avoid complex gas distributions in the background.All maps were projected onto the same 0. • 25-spaced Cartesian grid, that is appropriate for the γ-ray counts observed at very high Galactic latitudes.

Gamma-ray data
We used ten years of Pass 8 photon data provided by the LAT between 0.25 and 63 GeV (Atwood et al. 2013;Bruel et al. 2018).This energy range was chosen to preserve an angular resolution better than approximately 2 • for the 68% containment angle at low energy and to have sufficient photon statistics at high energy.We used the associated instrument response functions P8R3_SOURCE_V2 (for the various PSF event types), and the corresponding isotropic spectrum for the extragalactic and residual instrumental backgrounds 1 .Tight selection criteria were used: SOURCE class selection, PSF 1, 2 and 3 event types, photon arrival directions within 100 • or 105 • of the Earth zenith depending on the photon energy and PSF type (100 • below 0.6 GeV for PSF 1-2 and below 0.4 GeV for PSF 3, 105 • otherwise).Such criteria reduce the contamination by residual cosmic rays and by Earth atmospheric γ rays in the photon data.The instrument functions, the exposure map, the γ -ray emissivity spectrum of the local interstellar gas, q LIS (Casandjian 2015), and the spectrum of the isotropic background were evaluated in 12 energy bins, 0.2 dex in width and centred from 102.5 to 10 4.7 MeV.We took the energy resolution of the LAT into account when modelling the data.To ensure photon statistics robust enough to follow details in the spatial distributions of the different interstellar components, we analysed the data in seven broad and independent energy bands, bounded by 10 2.4 , 10 2.6 , 10 2.8 , 10 3.0 , 10 3.2 , 10 3.6 , 10 4.0 , and 10 4.8 MeV.

H I and CO data, as well as cloud separation
In order to trace the atomic gas, we used the 16. 2 resolution HI4PI survey (Ben Bekhti et al. 2016), with a velocity resolution of 1.49 km s −1 in the local standard of rest (LSR).We selected velocities between −90 and +50 km s −1 to exclude the H I emission coming from the high-velocity clouds that lie in the hot Galactic corona, far behind the local medium we are interested in (Wakker et al. 2008).We checked that these high-velocity clouds were not detected in γ rays and we removed them from our model.
In order to trace the molecular gas, we used the 8. 5 resolution 12 CO (J = 1−0) observations at 115 GHz from the momentmasked CfA CO survey (Dame et al. 2001;Dame & Thaddeus 2004).We completed this dataset with the CO observations of the MBM 20 cloud that were obtained with the Swedish-ESO Submillimetre Telescope (SEST).They were kindly provided by D. Russeil (Russeil et al. 2003).
We decomposed the H I and CO velocity spectra into individual lines and we used this information to identify and separate eight nearby cloud complexes that are coherent in position, velocity, and distance.Details of the method are presented in Paper I. The resulting maps were projected onto a 0. • 25-spaced Cartesian grid in equatorial coordinates.The main entities associated with the superbubble are the North Rim , West Rim , and South Loop along the outer rim, and the East Shell in the interior.We have also identified the Eridu 2 atomic cirrus cloud that lies outside the superbubble, at a comparable distance, but further away from the Galactic plane.The compact MBM 20 molecular cloud lies just in front of the superbubble edge, between the Local Bubble and the Orion-Eridanus superbubble.The relative positions of all these clouds in the sky are displayed in Fig. 1.The analysis region partially intercepts other cloud complexes, the edge of the nearby Taurus cloud, and part of the Cetus and North Taurus complexes which are located behind the superbubble (Remy et al. 2017).Their contributions to the γ-ray emission were taken into account in the model, as well as the faint atomic background from the Galactic disc.

Dust
In order to trace the dust column density, we used the optical depth map obtained at 353 GHz, τ 353 , by Planck Collaboration Int.XVII (2014).It was obtained by modelling the dust thermal emission recorded by Planck and the Infrared Astronomical Satellite (IRAS) with a modified black body, I ν = τ 353 GHz B ν (T d )(ν/353 GHz) β , with T d the dust colour temperature and β the spectral index.We degraded the optical-depth map from its original resolution of 5 down to 16.2 to match the resolution of the H I data and we projected it onto a 0. • 25spaced Cartesian grid.The dust opacity, or optical depth per gas nucleon, σ ν , is defined according to τ ν = I ν B ν (T ) = σ ν N H .

Ionised gas
Ionised gas is visible in H α emission.It is displayed in Fig. 1 using the data of Finkbeiner (2003) traces only the recombining gas and it does not linearly scale with the column density of ionised gas, N H II .Its intensity is proportional to the emission measure: I H α ∝ n i n e ds, with n i the ion density, n e the electron density and the integration is along the line of sight.Free-free emission at mm wavelengths can also trace the ionised gas, but the most recent estimates inferred from the Planck data are still too heavily contaminated by dust emission in our analysis region.We therefore created a uniform template in order to keep the morphology of the main H α arcs.
We set the template value to one within the 14 R contour of H α emission and zero outside.The normalisation was left free in our fits (see Sect. 3).We note this template K H α .
The hot ionised gas that fills the superbubble unfortunately yields column densities below 5 × 10 19 cm −2 that are too low to be detected in γ rays with the LAT.We thus did not include this gas in the model.

Models and analyses
The intensity of the hadronic γ radiation corresponds to the integral along the lines of sight of the CR flux times the gas volume density.The LAT data we have used probe CR with energies above a few GeV that uniformly permeate all gas phases up to the molecular phase seen in CO at 115 GHz (Grenier et al. 2015;Remy et al. 2017, and references therein).The γ-ray intensity therefore scales with the total gas column density, N H , and can be modelled by a linear combination of the column densities present in the ionised, atomic, molecular, and dark neutral phases.The average γ-ray emissivity per gas nucleon in a cloud can only be inferred in the H I phase where we can directly measure the hydrogen column densities, provided corrections for self absorption of the H I lines.This is not possible in the other gas phases where the knowledge of the CO-to-H 2 conversion factors and of the dust emission opacities is required to infer the gas mass.These factors are, moreover, known to vary from cloud to cloud and with environmental conditions (Remy et al. 2017).To model the total γ-ray emission, additional ancillary components needed to be considered like the Galactic inverse Compton (IC) radiation, point sources, the isotropic emission accounting for extragalactic and instrumental backgrounds, and the emissions from the Sun and the Moon.

The dark neutral medium (DNM)
The dark neutral phase contains large amounts of gas at the H I -H 2 interface, but it is invisible in H I and CO line emission because of H I self-absorption and because of the photodissociation and weak excitation of CO molecules in diffuse H 2 .In order to trace the DNM gas, we iteratively coupled the γ-ray and dust analyses as dust column densities also trace the total N H in the local ISM.This is illustrated in Fig. 2 where we convolved the dust optical depth map with the LAT PSF on the one hand, and, on the other hand, we subtracted from the γ-ray data the non-gaseous emissions obtained from our best-fit model (Sect.3.2). Figure 2 shows strong similarities in the spatial distributions of the dust and γ-ray gas tracers, but it also reveals differences in their dynamical range in several places.It has been shown in particular that the dust optical depth does not linearly trace the total gas column density in the dense molecular phase (N H > 3 × 10 21 cm −2 ) because the dust grains evolve and their emission properties change (Planck Collaboration Int. XVII 2014;Planck Collaboration Int. XXVIII 2015;Remy et al. 2017).The DNM phase for which we use dust data in our analysis is, however, diffuse enough for the linear approximation to hold.We therefore extracted from the dust and interstellar γ-ray maps the significant signals found above the H I and CO expectations and we used the spatial correlation between these signals to infer the A96, page 4 of 19 T. Joubaud et al.: The cosmic-ray content of the Orion-Eridanus superbubble additional gas column densities that correspond to the DNM at the H I -H 2 transition.

Gamma-ray model
Earlier studies indicated that the bulk of the Galactic CR radiating at 0.1-100 GeV have diffusion lengths much larger than typical cloud dimensions and that they permeate all the H Ibright, DNM, and CO-bright gas phases.The observed interstellar γ-ray emission can therefore be modelled, to first order, by a linear combination of the contributions of the different clouds seen along the lines of sight, including their different gas phases (H α , H I, DNM, CO).In order to convolve the gas maps with the LAT response functions, we assumed that the γ-ray emissivity spectrum of each gas phase of each cloud follows the average one obtained in the local ISM (q LIS (E), Casandjian 2015).We left, however, a free normalisation factor for the emissivity of each component in each energy band to account for possible deviations in CR density and spectrum.
The model includes other radiation components such as the Galactic IC intensity, I IC (α, δ, E), the isotropic intensity mentioned above, I iso (E), intensities from the Sun and the Moon, I SM (α, δ, E), and point sources with individual flux spectra S j (E).The γ-ray intensity I(α, δ, E), expressed in cm −2 s −1 sr −1 MeV −1 , can thus be modelled in each energy band as: where N H I,i (α, δ) denotes the H I column density map of each component, W CO,i (α, δ) gives the CO line integral intensity map, and K H α (α, δ) is the H α template described in Sect.2.4.τ DNM 353 stands for the DNM dust map resulting from the joint dust and γ -ray analyses (further detailed in Sect.3.4).To account for the spill-over of emission produced outside the analysis region, but reconstructed inside it, we modelled point sources and interstellar contributions in a region 6 • wider than the analysis region.
All the components are illustrated in Fig. 3. Three clouds have atomic envelopes and molecular cores (the superbubble North Rim, the foreground MBM 20 cloud, and the background Cetus-North Taurus clouds).The other five clouds (the superbubble South Loop , West Rim, and East Shell , the Eridu cirrus, and the edge of the foreground Taurus cloud) as well as the Galactic disc background are only detected in H I .The DNM map is dominated by dark gas in the North Rim and a small contribution from MBM 20.
The input q LIS spectrum was based on four years of LAT data and on the correlation between the γ radiation and the N H I column densities derived from the LAB survey, for a spin temperature of 140 K, at latitudes between 7 • and 70 • (Casandjian 2015).The q H I,i (E) normalisation factors in the model can therefore compensate for cloud-to-cloud variations in CR flux or spectrum.For each cloud, the final average γ-ray emissivity spectrum per gas nucleon in the atomic phase is the product of the q LIS (E) spectrum and of the best-fit q H I,i (E) normalisation in each energy band.
The model includes 151 points sources from the Fermi-LAT 8-yr FL8Y source list3 inside the analysis region.Their flux spectra, S j (E), were computed with the spectral characteristics provided in the catalogue.Their q S j (E) normalisation in the model allows for possible changes due to the longer exposure used here and to the use of a different interstellar background for source detection in FL8Y.The sources in the analysis region were fitted individually.The sources present in the 6 •wide peripheral band around the analysis region were merged into a single map, S ext (α, δ, E), and its normalisation, q S ext (E), was left free in each energy band.We checked that its best-fit normalisation was compatible with one.
The isotropic emission spectrum was determined with the FL8Y interstellar background model4 .The γ -ray emission observed in the region also includes a contribution from the large-scale Galactic IC emission emanating from the interactions of CR electrons with the Galactic interstellar radiation field.The GALPROP5 parameter file 54-LRYusifovXCO4z6R30-Ts150-mag2 was used to generate an energy-dependent template of the Galactic IC emission across the analysis region (Ackermann et al. 2012).The Sun and the Moon paths also cross this region, so we took their γ-ray emissions into account.This was done following the method described by Johannesson & Orlando (2013), reproducing the interactions of CR electrons with the solar radiation field and the hadronic interactions of CR nuclei with the Sun and the Moon.The Galactic IC emission model, the isotropic emission spectrum, and the Sun and Moon emission models were derived from earlier LAT data.Therefore, we left their respective scalings free in each energy band to allow for possible changes due to the longer exposure used here and to the use of a different interstellar background.
In order to compare the modelled photon map predictions with the LAT photon data in the different energy bands, we multiplied each component map by the energy-dependent exposure and we convolved it with the energy-dependent PSF.We also convolved the spectra in each pixel with the LAT energy resolution.The convolution of the gas maps with the response functions was done for each PSF event type independently in small energy bands.We then summed the photon count maps in energy and for the different PSF types to get the modelled count maps in the seven energy bands of the analysis.Figure 3 gives the photon yields that were obtained in the entire energy band from the best fit.It shows that the emission originating from the gas dominates over other types of emissions and that the LAT angular resolution makes it possible to separate the various clouds, as well as the different gas phases within the clouds, based on their spatial distribution.The q model parameters were fitted to the LAT photon data in each energy band using a binned maximum-likelihood with Poisson statistics.

Dust model
In the case of a uniform dust-to-gas mass ratio and uniform mass emission coefficient of the grains, the dust optical depth linearly scales with the total N H .We therefore modelled τ 353 (α, δ) in each direction as a linear combination of the same gaseous contributions as in the γ-ray model, with free normalisations to be fitted to the data, as in Planck Collaboration Int.XXVIII (2015); Remy et al. (2017).We added a free isotropic term, y iso , to account for the residual noise and the uncertainty in the zero level of the dust data (Planck Collaboration XI 2014).The τ 353 (α, δ) model can be expressed as where N H DNM (α, δ) stands for the DNM gas column density map iteratively built from the coupled analyses of the γ-ray and dust data (further detailed in Sect.3.4).
The y model parameters were estimated using a χ 2 minimisation.We expect the model uncertainties to exceed the measurement errors in τ 353 (α, δ) because of potential variations in grain properties through the clouds and because of the limitations of the gas tracers (survey sensitivities, emission saturation, self-absorption, etc.).As we cannot precisely determine the model uncertainties, we set them to a fractional value of the model and we determined this fraction to be 18% by reaching a reduced χ 2 of unity.This fraction is larger than the 3-9% measurement uncertainties in the τ 353 values across this region (Planck Collaboration XI 2014).

Analysis iteration
In order to extract the DNM gas present in both the dust and γ-ray data in addition to the gas seen in H I and CO emission, we iteratively coupled the γ-ray and dust models.We built residual maps between the data (dust or γ rays) and the best-fit contributions from the N H II , N H I , W CO and other non-gaseous components.We kept only the positive residuals with high significance above the noise (a simple cut at zero would induce an offset bias in the next model).To do so, we de-noised the residual maps using the multi-resolution support method implemented in the MR filter software (Starck & Pierre 1998).We used a multi-resolution thresholding filter, with six scales, Gaussian noise and a 2σ threshold to filter the dust residuals.For the γ-ray ones, we used seven scales, Gaussian and Poisson noise, and a 3σ threshold.The DNM template estimated from the dust was provided to the γ-ray model (τ DNM 353 in Eq. ( 1)); conversely, the DNM column density map derived from the γ rays was provided to the dust model (N H DNM in Eq. ( 2)).We started the iterative process by fitting the dust optical depth with the H I , CO, H α and isotropic components to build the first DNM map for the γ-ray model.We then iterated between the γ-ray and dust models until reaching a saturation in the log-likelihood value of the fit to the γ-ray data (from the third to the fourth iteration).We checked that the resulting DNM maps were consistent in morphology and column density with the maps obtained by Remy et al. (2017) in the intersection region between the two analyses.
At each iteration, the estimates of the q and y model coefficients as well as the DNM maps changed, and the likelihood significantly improved since there was less and less need for the other components, in particular the H I and CO ones, to compensate for the missing gas.They still do at some level because the DNM template provided by the γ rays or dust optical depth have limitations (e.g., dust opacity variations, limited γ-ray sensitivity).

Results
Before presenting in Sects.4.7 and 4.8 the γ -ray emissivity spectra obtained in the different clouds from the best-fit models described in Eqs. ( 2) and (1), we discuss how the fits to the γ -ray and dust data improved with optical depth corrections to the atomic N H I column densities (Sect.4.1) and with the addition of gas templates other than H I and CO (Sects.4.2 and 4.3).We also assess the robustness of the best-fit models with jackknife tests in Sect.4.4 and we present residual maps in Sect.4.5.

H I optical depth correction
In the range of N H I column densities studied here, Nguyen et al. (2019) showed that a simple isothermal correction of the H I emission spectra, with a uniform spin temperature (T S ) across the cloud, provides better than 10% estimates of the more precise N H I values inferred from the combination of emission and absorption H I spectra.Other correction methods are less efficient.In order to account for the unknown level of H I opacity in the different N H I maps, we have repeated our dust and γ-ray analyses for ten spin temperatures (100, 125, 150, 200, 300, 400, 500, 600, 700, 800 K) and for the optically thin case which yields the minimum amount of gas.We applied the same spin temperature to all the clouds as we could not test all the possible combinations of individual spin temperatures.The best-fit residual maps show that a single spin temperature can reproduce the γ-ray data without any significantly deviant cloud.
The γ rays can help constrain the average level of H I opticaldepth correction applicable to the whole region by comparing the T S -dependent contrast of the N H I maps with the spatial structure of the γ-ray flux emerging from the H I gas.We found that the maximum log-likelihood value of the fits significantly preferred a spin temperature of 100 K in this region (see Appendix A).We could not test spin temperatures below 100 K in our isothermal approach as its value needs to be above the highest values of the brightness temperature, which are close to 100 K in our region.In the following, we present the results obtained for this temperature of 100 K, unless otherwise mentioned.

DNM detection
In order to assess the presence of additional gas that is not traced by H I and CO line intensities, we first performed the dust and γ-ray fits with only the H I, CO and H α templates as gaseous components.The residuals are displayed on the left-hand side of Fig. 4. The data significantly exceed the model in several extended regions, coherently in dust and in γ rays though the two fits are independent.A coincident variation in CR density and in dust-to-gas ratio is very unlikely and hadronic CR interactions with dust grains produce undetectable γ-ray emission (Grenier et al. 2005), so this correlated excess has a gaseous origin.The DNM gas distribution is built by correlating and filtering the excess dust and γ-ray signals in the iteration between the dust and γ-ray fits (see Sect. 3.4).Adding this template improves the fits quality (see the right-hand side of Fig. 4).An improvement with a log-likelihood ratio of 5418 when adding the DNM template to the γ-ray fit implies a very significant DNM detection.
The left panel of Fig. 4 shows that most of the DNM gas is within the North Rim.A small contribution comes from the MBM 20 cloud around (α,δ) = (68, −14).Keeping the two clouds in the same DNM map assumes a common dust opacity and γ-ray emissivity for both.We have not split the map in two since the DNM mass in MBM 20 is negligible compared to the North Rim one and it does not bias the CR results about the North Rim.
More details on the DNM column densities and mass will be presented in a forthcoming paper, together with a dedicated analysis of MBM 20 at a higher spatial resolution and with an independent DNM template.These results yield much improved dust residuals towards MBM 20 and they give support to keeping the DNM map as a whole for the superbubble analysis.We do not present any CR result for MBM 20 here as its small size requires the higher angular resolution analysis to separate the γ-ray contributions from the different gas phases.

Detection of ionised gas
The H α template was added to the models to take into account hadronic CR interactions in the arcs of warm ionised gas.In γ rays, this component is weakly detected at a 4.6σ confidence level over the whole energy band, which is low compared A96, page 7 of 19 A&A 635, A96 (2020) to the ∼74σ detection for the DNM.The H α template is, moreover, only detected below 2 GeV because of too low photon yields at higher energies from this faint component.Figure 3 shows that this component contributes less than 2 photons per pixel (0. • 25 × 0. • 25) over the entire 0.25-63 GeV energy band.This contribution can be converted into an N H column density assuming the same γ-ray emissivity per gas nucleon as in the North Rim .The result gives values around 7 × 10 19 cm −2 that are realistic, but at the limit of detection with the LAT.These values are much lower than the 2 to 30 × 10 20 cm −2 column densities that characterise the other atomic clouds in the field.We have kept the H α template in the models to avoid biasing the other components, but we do not further discuss the physical implications of this faint component.

Best γ-ray and dust fits and jackknife tests
The best-fit model to the γ -ray data includes all the components described in Eq. ( 1) as they have all been significantly detected.The faint N H I column densities from the Galactic-disc background have not been detected in the dust fit at low spin temperatures because of the overlap with the much more massive North Rim, but all the other components give significant contributions to the best-fit model.
The χ 2 minimisation and binned maximum-likelihood methods used respectively for the dust and γ-ray fits give access to statistical errors on the best-fit coefficients.They are inferred from the information matrix of the fit (Strong 1985) and they include the effect of correlations between parameters.Because of the large spatial extents of the cloud complexes and because the gas, dust, and γ-ray distributions are tightly correlated, the statistical errors are small: typically 1-6% and 1-7% for the gas parameters of the local clouds in dust and γ rays, respectively.Only the contributions of the Galactic background and H II template to the γ-ray data have larger uncertainties of 9 and 17%, respectively.
We also checked the magnitude of the systematic uncertainties in the linear approximations of the models, hence of spatial changes in the model parameters and/or in the mean level of H I and CO self absorption.We did so by repeating the last iteration of the dust and γ-ray fits a thousand times over random subsets of the analysis region, namely after masking out 20% of the pixels with a sum of 3. • 75-wide, randomly selected squares.In γ rays, the jackknife tests have been performed only for the total 0.25-63 GeV energy band.We have assumed the same relative deviations, σ q /q, for the individual energy bands as for the total one since the jackknife sampling tests potential non-uniformities in the models at larger scales than the angular resolution of the data.
Appendix B shows the distributions thus obtained for the best-fit coefficients.Most of them show Gaussian-like distributions, thereby indicating that the results presented in Appendix C are statistically stable and that the average coefficients that describe our models are not driven by subset regions in each cloud complex.A few distributions exhibit long, non-Gaussian tails when the corresponding clouds subtend small solid angles (e.g., the small CO clouds from the North rim, or the Sun and Moon emissions).The long tails reflect the indeterminacy of the parameter when a large fraction of the non-zero pixels in the corresponding maps is masked.The standard deviations found in the jackknife distributions amount to 1-9% in dust and to 1-5% in γ rays, except for the H I Galactic background (11%) and the H II template (23%) in γ rays.
We quadratically added the 1σ fitting errors and the corresponding standard deviations of the jackknife distributions to give the final statistical errors on a given parameter.The results on the best-fit values and final errors are given in Appendix C.

Residual maps
The residual maps obtained from the best dust and γ-ray fits for the final iteration are presented in Fig. 4 (right-hand side).When including all gaseous components, in particular the DNM one, we see that our linear model provides an excellent fit to the γ-ray data (bottom right).The residuals are consistent with noise except, marginally, towards a compact CO cloud near (α, δ) = (72,18) where the model slightly over-predicts the data.This may be due to a lower N H 2 -to-W CO conversion factor in this cloud compared to the larger values that characterise the other, more tenuous, clouds in the North Rim molecular map.Such variations have been observed in other molecular clouds (Remy et al. 2017).
The dust residual map (top right) presents structured deviations.The ring-like excess and central deficit towards MBM 20 is due to the 0. • 25 sampling resolution of the present analysis that is not sufficient to capture the multi-phasic structure of this very compact molecular cloud.These residuals vanish in the higher angular resolution analysis that will be presented in the forthcoming paper.The other excesses, towards the North Rim and West rim, are likely due to the non-linear increase in dust opacity in the dense parts of these clouds and will be discussed in the forthcoming paper.

Cosmic-ray penetration in the different gas phases
For a given cloud, spectral changes in the ratio between the γ -ray emissivities measured in the DNM and CO-bright gas and that in the diffuse atomic envelope can reveal changes in the underlying CR distribution as the particles penetrate into the molecular phases.One expects low-energy CR (below 10 GeV) to be depleted as they diffuse away on self-excited MHD turbulence in the cloud envelope.The turbulence is triggered by the CR density gradients induced by severe hadronic losses in the cloud interior (Skilling & Strong 1976;Ivlev et al. 2018).
Figure 5 shows that we found a stable DNM to H I emissivity ratio as a function of energy in the massive North Rim of the superbubble.It implies that the same CR spectrum permeates the atomic and DNM phases even though the gas is denser and less ionised at the H I -H 2 transition than in the warm H I envelope.The North Rim DNM mass is large (of order 1.2 × 10 4 M at 160 pc), but the DNM gas column densities of 2 to 20 × 10 20 cm −2 are one to two orders of magnitude below the critical values for significant low-energy CR depletion (Ivlev et al. 2018).
Likewise, for the North Rim and North-Taurus-Cetus clouds, we found no spectral evolution of their respective CO to H I emissivity ratios.The weak down-going trend seen for the North Rim ratio is not statistically significant.The data is consistent with the same CR spectrum flowing through the H I envelope and CO-bright gas, the latter having charateristic volume densities of a few 10 3 cm −3 and magnetic field strengths of 1-10 µG (Crutcher 2012).The molecular gas column densities in these clouds are indeed far below the critical values for low-energy CR depletion.These results confirm the same trends found in other nearby clouds (Grenier et al. 2015;Remy et al. 2017, and references therein).

Cosmic rays in the Orion-Eridanus superbubble
We first probed the average γ-ray emissivity, q SB , of the gas in the superbubble by gathering the four superbubble clouds into a single component.This emissivity was derived for a spin temperature of 150 K in order to compare with the average emissivity, q LIS , derived in the local ISM for a spin temperature of 140 K (Casandjian 2015).The result is displayed in Fig. 7.The superbubble emissivity presents the same spectral energy distribution (SED) as the local ISM average.The 6% difference is highly significant (12σ) over the entire energy band.Figure 7 shows, however, that the difference is consistent with the systematic uncertainties in the H I optical depth of the different clouds that have served to perform both measurements.The difference is also compatible with the ±9% dispersion found in other nearby clouds for the same reason (Remy et al. 2017).The amplitude and spectral shape of the mean q SB emissivity in the superbubble implies that the CR spectrum pervading this region is not modified, in particular at high energy, despite the perturbed conditions and past supernova history of the superbubble.
We compared the γ-ray emissivity spectrum in the superbbuble with that inferred from the CR spectrum measured in the Fig. 7. Average spectral energy distribution of the γ-ray emissivity per gas nucleon, q SB , found in the atomic gas of the superbubble for an H I spin temperature of 150 K (magenta diamonds, solid lines) and for the optically thin H I case (magenta circles, dashed lines).The blue squares show the average emissivity, q LIS , found in the local ISM for a spin temperature of 140 K (Casandjian 2015).The black line gives the γ-ray emissivity per gas nucleon inferred from the CR spectrum measured in the Solar System (see text).
Solar System.To do so, we used the parametrisations of the CR proton and helium spectra proposed by Corti et al. (2016) on the basis of the Voyager 1, AMS-02, PAMELA, and BESS data and the use of an updated, rigidity-dependent, force-field approximation to account for the solar modulation.We also used the CR electron spectrum inferred from the Voyager 1, AMS-02, PAMELA, and Ulysses data and the HELMOD code for CR propagation inside the heliosphere (Boschini et al. 2018).We used the hadronic-interaction cross sections from Kamae et al. (2006) and the multiplication factors to take elements heavier than H into account from Mori (2009, see also Kachelriess et al. 2014).The result for the γ-ray emissivity per gas nucleon in hadronic and bremsstrahlung interactions is presented in Fig. 7.It shows that the CR spectrum diffusing through the superbubble is compatible within the uncertainties with the particle distribution passing through the Solar System.
We probed the uniformity of the γ-ray emissivity across the superbubble.The emissivity SED found in individual H I clouds associated with the superbubble are shown in Fig. 8 for the spin temperature of 100 K favoured by the γ-ray fit for this region.They are compared with the average q SB emissivity in the superbubble in order to search for spatial variations.The emissivity SED in the North Rim and South Loop are consistent with q SB at all energies within the uncertainties.The slightly lower emissivity found in the East Shell at the lowest and highest energies is not statistically significant (0.6σ).The 8% lower emissivity in the West Rim differs from the q SB average at a 5.6σ confidence level, but the result for the optically thin case shows that the small difference may stem from cloud to cloud variations in the H I optical depth correction.

Cosmic rays in the Eridu cirrus
The Eridu cirrus cloud lies outside the superbubble, at an altitude of about 200-250 pc below the Galactic plane.Its distance was derived in Paper I using the dust reddening map of Green et al. (2018).The SED of the γ-ray emissivity we found in this cloud is displayed in Fig. 9. Compared with the local ISM q LIS average, we found a highly significant (14σ) drop in emissivity by 34%.This difference cannot be attributed to H I optical Fig. 8. Spectral energy distributions of the γ -ray emissivity per gas nucleon found in different atomic clouds associated with the superbubble, for an H I spin temperature of 100 K favoured by the fit (solid, diamonds) and for the optically thin H I case (dashed, circles).The grey squares show the average emissivity in the superbubble, q SB , for an H I spin temperature of 100 K.
depth corrections as the optically thin case, which sets the minimum amount of gas, gives an upper limit to the emissivity that is still 25% below q LIS .We do not detect any significant amount of DNM gas in the Eridu direction, but adding any missing gas would imply an even smaller γ-ray emissivity per H nucleon to produce the same γ-ray intensity.Confusion in the cloud separation in velocity is not possible as it is the only H I structure seen towards these directions.At Galactic latitudes near −50 • , it is indeed too far from the Galactic plane for significant background gas contamination.
Given the modest gas column densities in the cloud (1 − 5 × 10 20 cm −2 ), we checked the influence of large-scale γ-ray structures from the IC and isotropic emission components.Varying those by 3σ around their best-fit spectra yields a maximum variation of 6% for the Eridu emissivity spectrum.We also assessed the influence of the spatial distribution of the Galactic IC emission across the region by trying different GALPROP models.The radial distribution of CR sources in the Galaxy as well as the halo size directly influence the IC intensity gradient across the analysis region, but in a limited way.Relative variations in intensity do not exceed 10% if we compare the case of a steep Galactic CR gradient (with CR sources following the peaked distribution of OB stars and for a small halo height of 4 kpc) and the case of a shallow CR gradient (with CR sources following the distribution of supernova remnants and for a halo height of 10 kpc) (Ackermann et al. 2012).The Eridu emissivity spectrum changes by less than 3% if we replace the present IC model by one of those more extreme cases.These tests show that the spatial distribution of the gas in the Eridu cloud is significantly detected in the γ-ray fit, with little influence from the large-scale structures in non-gaseous components.
Fig. 9. Spectral energy distribution of the γ -ray emissivity per gas nucleon found in the atomic gas of the Eridu cirrus, for an H I spin temperature of 100 K favoured by the fit (solid, diamonds) and for the optically thin H I case (dashed, circles).The blue squares show the average emissivity, q LIS , found in the local ISM for a spin temperature of 140 K (Casandjian 2015).

H I optical depth correction
Among the set of spin temperatures we used to account for H I optical depth corrections, the maximum log-likelihood values favour a temperature of 100 K, as shown in Appendix A. This temperature is lower than in other nearby clouds (Planck Collaboration Int. XXVIII 2015;Remy et al. 2017).It is also low compared to the values closer to 300 K expected for N H I column densities in the range of (2−30) × 10 20 cm −2 sampled here.Such column densities only partially shield the gas from the interstellar UV and soft X-ray radiation field (0.1 A V 1.5) (Kanekar et al. 2011;Bohlin et al. 1978).Our model assumes a single spin temperature over the entire line of sight line (monophasic T S ) and in the different clouds.The latter implies that the spin temperature determination is biased by the most extended clouds because of their preponderant weight in the fit.Among the different clouds in the field, the main H I contribution comes from the North Rim which, being rich in DNM gas at the dense H I and diffuse H 2 transition, likely contains a large fraction of cold neutral medium structures.The latter have spin temperatures well below 200 K (Murray et al. 2015).
It should, however, be noted that, due to the relatively low column densities involved here, the optical depth correction has only a small impact of order 10% on the CR spectrum derivation.The relative differences in γ-ray emissivities between the T S = 100 K and the optically thin cases are 1, 10, 10, 17, and 14% for the East Shell , West Rim, South Loop, North Rim, and Eridu cirrus, respectively.

The cosmic-ray flux in the superbubble
The CR spectrum, probed here at energies between approximately 2 and 500 GeV, appears to be uniform across the solar neighbourhood and the older end of the Orion-Eridanus superbubble.It is consistent with the direct CR measurements in the solar system.This energy spectrum in E −2.7 indicates that CR propagation near the superbubble environment should be dominated by the same energy-dependent mode of diffusion as in the Galaxy at large.We note that the present γ-ray emissivity measurements mainly probe the CR flux in the outer shells of the superbubble and not directly in its hot interior plasma.The latter yields too low column densities to be detectable in γ rays with the LAT, and even the denser warm ionised gas is not firmly detected in our analysis.CR acceleration or re-acceleration may then take place in the superbubble interior or in the more active zone at the younger end of the superbubble (not studied here), but it leaves no detectable trace in the outer shells closer to us that we have studied.The emissivity measurement in the East Shell suggests that the CR flux inside the superbubble is consistent with that of the local ISM as this cloud is considered to be inside the hot gas.
Galactic CR may not penetrate deep into the superbubble if there is strong turbulence inside it (Tolksdorf et al. 2019).Conversely, for an inefficient accelerator as the Orion-Eridanus superbubble seems to be, CR can pass through the hot plasma without significant changes.Yet, the magnetic field lines appear to be wrapped around the superbubble, having been ordered and compressed along the outer rim by the shock expansion.Figure 10 shows the orientation of the plane-of-the-sky magnetic field (B sky ) in the West Rim , where strengths up to 15 µG have been inferred (Paper I).CR entering the shell will mainly diffuse along these field lines.
Transport perpendicular to the average field lines is possible with the field line random walk (Jokipii 1966).Particles can also drift across the field because of large-scale variations in magnetic direction or strength, or by scattering off smallscale MHD structures (Desiati & Zweibel 2014).Perpendicular transport is, however, always much slower than parallel transport.Even if advected with Alfvén waves, the bulk of CR would stream at the Alvén speed along the large-scale field lines.The ordered magnetic field in the outer shell could then limit the mixing between the internal and external CR populations.On the other hand, we expect the older part of the superbubble to be filled with numerous weak secondary shocks, reflected off the corrugated outer shell and by interior cloudlets, as observed in simulations (Kim et al. 2017).This turbulence, scaling as k −2 , may help CR to diffuse in an energy-independent way across the outer shell despite its large-scale magnetic configuration (Bykov & Toptygin 1987).
In this context, the comparison between the γ-ray emissivities of the East Shell and West Rim can provide insight into CR transport.The East Shell is considered to be inside the hot plasma phase because of the lack of X-ray absorption against the bright interior plasma (Paper I), but the distance uncertainties may as well place it as an outer-rim shell in the lower, rear part of the superbubble.Future absorption studies, in X rays and with dust, will be crucial to precisely locate this cloud and its γ-ray flux measurement with respect to the edge of the superbubble.
The lack of a CR enhancement or hardening in the Orion-Eridanus superbubble is in marked contrast with the case of the Cygnus X cocoon where a bright E −2.1 γ differential emissivity spectrum has been detected from 1 to 100 GeV (Ackermann et al. 2011).We discuss below possible origins of the apparent inefficiency of the Orion-Eridanus superbubble as a potential CR (re)-accelerator.

Intermittent shock acceleration
The first possible origin concerns the acceleration conditions inside the superbubble.The acceleration efficiency by intermittent shocks and magnetic turbulence depends on the spectral energy density of the magnetic fluctuations that results from the complex interplay between energy cascading, wave dissipation, and CR-driven instabilities.The outcome is not known and the spectrum is often modelled as a k −b T power law in wavenumber, scaling from the injection scale L * .For strong magnetic turbulence, with an amplitude η T = δB 2 /(B 2 + δB 2 ) ≈ 1 close to saturation, the particles with gyroradius R g and velocity v diffuse isotropically with a diffusion coefficient A&A 635, A96 (2020) that largely exceeds the Bohm value D B = R g v/3.The stellar density in the superbubble and the results of superbubble MHD simulations (Ntormousi et al. 2017) suggest a typical injection scale L * ≈ 10 pc and interior magnetic-field strengths of order 1-3 µG.In these conditions, the diffusion coefficient spans from for a Kraichnan-like spectrum b T = 3/2 (Perez et al. 2012), to a tenfold larger value These values are respectively 200 and 20 times smaller at 10 GeV than the average uniform coefficient inferred in the Milky Way using the GALPROP propagation code and elemental spectra from AMS-02 and ACE/CRIS (Jóhannesson et al. 2019).The size of the TeV haloes found around the nearby Geminga and PSR B0656+14 pulsars requires slow diffusion between 1.5 and 7.0 × 10 26 (E/10 GeV) cm 2 s −1 (Abeysekara et al. 2017).These values compare well with our superbubble estimates and they apply to the same distance range and anticentre direction.The TeV data is best described with a 3 µG ambient field and a coherence length 5 pc of the magnetic turbulence (López-Coto & Giacinti 2018).These conditions are close to the values used for the superbubble.Yet, even though Geminga is born in the superbubble (Pellizza et al. 2005), it now lies 100 pc outside it, as do PSR B0656+14, and the extent of the slow diffusion zone around the pulsars is still debated.It could be restricted to the relic pulsar wind nebula (Tang & Piran 2019) or to a region 100 pc controlled by streaming instabilities (Jóhannesson et al. 2019, and references therein).So, one cannot directly compare diffusion lengths inside the superbubble and around the two pulsars.
In the Kraichnan case, the confinement time, L 2 /[6D], in the L ≈ 200-pc-wide bubble is long enough (>3.9 Myr) to accelerate particles from supra-thermal energies up to a moderate maximum energy of 30-80 GeV set by D(E max ) ≈ u sh L * /3 (Bykov 2001) if we assume typical velocities u sh ≈ 500 km s −1 for the shocks and large-scale turbulent flows in the old part of the superbubble.The case of Kolmogorov turbulence is much less favourable.The maximum energy is limited to a few GeV as the particles diffuse out of the system in less than one Myr.Ferrand & Marcowith (2010) reached the same conclusion about a fast escape restraining CR acceleration in a 300 pc-long superbubble filled with 1 µG plasma around Orion OB1.
Acceleration is slow in the weakly magnetised interiors of superbubbles, but the spectral index of the turbulence at the small scales that govern how fast the particles diffuse away also plays an important role in the overall acceleration efficiency.Hence the Orion-Eridanus superbubble may be less efficient than the Cygnus X one because of a weaker magnetic field inside the bubble and/or a steeper turbulence spectrum.MHD superbubble simulations can inform us on the turbulence spectrum inside CR-inefficient bubbles like Orion-Eridanus, but one needs to include CR feedback processes to model the magnetic turbulence inside active systems like Cygnus X (Bykov 2014).The feedback involves streaming instabilities, resonant wave-CR interactions, and the energy transfer from the plasma to CR.

Time evolution
The non-linear models by Bykov (2001Bykov ( , 2014) ) include intermittent shocks, a description of the correlations between the shocks and the long-wavelength turbulent motions of the compressible plasma, and energy conservation between the background plasma and CR population.The results show that the Orion-Eridanus system has the potential to accelerate CR up to a few hundred GeV (see Fig. 2 of Bykov 2001).The CR spectrum, however, evolves considerably with time, so one should consider the late development stage of the Orion-Eridanus superbubble.Even though the supernova rate is roughly constant over tens of Myr for stellar populations following the local initial mass function (Ferrand & Marcowith 2010), the models predict a softhard-soft evolution of the output CR spectrum with time.The late softening is due to the damping of the shocks and of the large-scale compressible turbulence because of the energy transfer to CR.It is not due to adiabatic losses to the superbubble expansion as the models are quasi-static because of the slow expansion.This is a valid assumption given the present expansion velocity of 20 km s −1 of the Orion-Eridanus superbubble (Paper I).Radiative losses are also unimportant in the rarefied plasma inside the superbubble (n H = 0.005−0.05cm −3 , Paper I).Given its age, the Orion-Eridanus system is now in the late softening stage tending towards p 2 f (p) ∝ p −3 .The hardest period (∝p −2 ) occurred earlier, after 1 Myr of evolution in the nonlinear model of Bykov (2001) and after 10 Myr in the linear version of Ferrand & Marcowith (2010).
One can further test this spectral ageing by investigating the younger part of the superbubble, closer to Orion.The stellar history in the superbubble, as pictured by the stream of blue stars spanning from the closest wall (∼20-Myr old stars) down to Orion (∼1-Myr old stars), indicates that the energy injection progressed in space over a period of several Myr (Bouy & Alves 2015).The X-ray emission from the hot interior plasma, tracing its current energy content (Krause et al. 2014), presents a clear hardening towards Orion (Ochsendorf et al. 2015, Paper I).The energetic sources have shifted to larger distances with time.But the gradual propagation of massive-star formation may not readily explain the CR contrast between Orion-Eridanus and Cygnus X since star formation has also progressed across Cygnus X over a comparable period, between 10 and 1 Myr ago (Berlanas et al. 2018).

Stellar content
A notable difference between the inefficient Eridanus and efficient Cygnus X accelerators is the presence of a super-massive cluster in the latter : Cygnus OB2 still gathers 143 supernova progenitors (>8 M , Wright et al. 2015).A similar total of 62+48 is found across Orion OB1 and the blue stream (Bouy & Alves 2015), but Cygnus OB2 has produced ∼6 SN/Myr for the last Myr (Lingenfelter 2018), whereas 1 SN/Myr has exploded on average in Orion-Eridanus over the last 12 Myr (Voss et al. 2010;Bouy & Alves 2015).The higher space density and frequency of supernovae in the last Myr in Cygnus X may have favoured CR acceleration compared to Orion-Eridanus.Colliding shocks, in particular, are more frequent in a rich and compact cluster such as Cygnus OB2.The interaction of a supernova shock with one or several fast approaching and massive stellar winds can rapidly accelerate CR to PeV energies with a hard p 2 f (p) ∝p −1 spectrum (Bykov et al. 2013).Their diffusion to >30 pc distances in less than 10 4 years after they escape from the colliding shocks would produce a hard CR cocoon in Cygnus X that has not recently occurred in Orion-Eridanus where the stellar density is lower.
Future Gaia studies will help to better constrain the stellar history of Orion-Eridanus as, unfortunately, many of the Orion OB stars are too bright to be part of the present Gaia DR2 survey (Gaia Collaboration 2018).

Superbubble size
As mentioned in the introduction, the massive OB association G25.18+0.26 may power a cocoon of hard CR (Katsuta et al. 2017).At a distance of 7.7 kpc, its projected size of 210 × 170 pc compares with that of Orion-Eridanus, so the different CR behaviours of the two systems do not simply relate to their comparable sizes.The models indeed predict that the efficiency of the power transfer from turbulence to CR always exceeds 10% after the first Myr, irrespective of the superbubble size (from 50 to 220 pc).The larger systems are more efficient at late times because of a reduced escape rate (Bykov 2001(Bykov , 2014)), but the difference is small (factor of ∼2).The next generation of γ-ray telescopes, such as CTA, will play a key role in extending the present sample of superbubble observations and in assessing the influence of their size and of their stellar-cluster content on CR acceleration.With higher resolution imaging, the investigations for hard CR should be directed towards sources like Westerlund 2, barely resolved by H.E.S.S. (H.E.S.S. Collaboration 2011), towards the confused case of Westerlund 1 where the GeV and TeV emissions appear to have different origins (Ohm et al. 2013), towards NGC 3603 which was not resolved by Fermi-LAT in the 4FGL catalogue6 , or towards the bright Arches, Quintuplet, and Cl * 1806-20 stellar clusters near the Galactic Centre, as suggested by Ackermann et al. (2011), H.E.S.S. Collaboration (2012), H.E.S.S. Collaboration (2018), Aharonian et al. (2019).

Cosmic-ray re-acceleration
We found no γ-ray evidence of re-acceleration of Galactic CR in the Orion-Eridanus superbubble.The bubble expansion has swept-up, compressed, and wrapped the ambient ISM magnetic field to the point where the field lines appear to be mostly tangent to the outer rim and perpendicular to the shock normal (Soler et al. 2018, Paper I and Fig. 10).This was considered as an unfavourable configuration for CR acceleration by the outer shock wave (Ellison et al. 1995).Recent simulations indicate, however, that the magnetic amplification by streaming instability can rearrange highly inclined fields ( 70 • ) to quasi-parallel configurations that allow particle injection into the shock for diffusive re-acceleration (Caprioli et al. 2018).The downstream medium may also be rich in supra-thermal particles due to the series of past supernova remnants in the superbubble and they can get re-accelerated.Tolksdorf et al. (2019) have analytically studied stochastic re-acceleration in the turbulent interior of a spherical superbubble, including spatial and momentum diffusion as well as pion losses.The acceleration potential is captured by the ratio of the acceleration time scale, τ A ∼ 9D(p)/v 2 A , over the diffusive escape timescale, τ D ∼ R 2 SB /D(p), with R SB the superbubble radius, v A the Alfven velocity in the interior plasma, and D(p) the spatial diffusion coefficient estimated in Sect.5.2.1.The ratio indicates if CR have spent enough time in the turbulent region to be re-accelerated before diffusing out.Pion losses play a negligible role in the Orion-Eridanus interior as their timescale is orders of magnitude longer than τ A and τ D .
Fig. 11.Theoretical spectral energy distributions of the γ-rays resulting from pion decay inside a turbulent superbubble for different values of the acceleration efficiency parameter ξ (according to Tolksdorf et al. 2019).The efficient case (ξ = 0.1) is in blue, the intermediate one (ξ = 1) in green and the inefficient case (ξ = 10) in orange.All distributions are normalised to one at the pion mass.For a better readability, the SED for ξ = 1 and ξ = 10 are divided by a factor of two and four, respectively.
Our γ-ray emissivity measurements in the shell are close to the spectral index of 2.6 predicted for an inefficient system with ξ 1 (see Fig. 11).Yet, the superbubble characteristics (R SB ≈100 pc, B int ≈ 1−3 µG, v A ≈ 25−80 km s −1 , Ntormousi et al. 2017, Paper I) yield an efficiency 0.2 < ξ(10 GeV/c) < 1.2 suggesting that CR re-acceleration is possible up to energies of 10-100 GeV if the L * injection scale of magnetic turbulence is close to 10 pc and the turbulence index approaches b T ∼ 3/2. Figure 11 shows that a system with an intermediate efficiency ξ ∼ 1 would yield a hard spectrum with a spectral index of about 0.7 that is much harder than what we observe.A steeper, Kolmogorov-type of turbulence would be much less efficient in the same conditions (3 < ξ(10 GeV/c) < 12).A stronger magnetic field in the hot plasma would lower the ξ ratio.Assuming equipartition with the hot gas pressure (Paper I), the field strength of 14 µG would clearly favour re-acceleration for both types of turbulence scaling (0.02 ≤ ξ(10 GeV/c) ≤ 0.3 for 3/2 ≤ b T ≤ 5/3).Substantially reducing the superbubble radius to increase ξ is not supported by the data and the L * injection scale has only a limited impact on the diffusion coefficient (D ∝ L * 1/2 ).Therefore, in order to reach an efficiency ξ(10 GeV/c) ∼ 10 that yields a CR spectrum close to our observations, the required turbulence level inside the superbubble should be η T 30% for Kraichnan turbulence.Hence the apparent re-acceleration inefficiency of the superbubble points to a low η T amplitude of the large-scale turbulence and/or to a steep index close to 5/3 of the power density of the fluctuations.Figure 12 shows that, for inefficient systems like the present one, Galactic CR can penetrate deep into the rarefied medium of the superbubble.

Eridu cirrus
The Eridu cirrus, located at an altitude of about 200-250 pc below the Galactic plane, lies next to, but outside the superbubble.It presents a γ-ray emissivity 34% lower than the local ISM average, and 28% lower than the superbubble average.We have investigated several leads that could explain a lower CR flux in this cloud.
First, the loss could result from partial penetration into the cloud.The Eridu cirrus has, however, a modest column density (1−5 × 10 20 cm −2 ) compared to the nearby superbubble clouds that are pervaded by a larger cosmic-ray flux.Everett & Zweibel (2011) studied the evolution of the CR pressure in a diffuse cloud with gas and magnetic-field conditions similar to those in the Eridu cirrus.They used an inner gas density of 100 cm −3 and a uniform magnetic field of 3 µG compared to the density of ∼7 cm −3 and B sky values of 5 µG in Eridu (Paper I).They found that the Alfven-wave pressure increases by nearly two orders of magnitude at the cloud periphery because of the growth of the streaming instability.But the density of low-energy particles (1 GeV/nucleon) drops by about 5% only in the cloud interior.There would be even less contrast for higher-energy particles.A 30% drop as observed in Eridu would require a much larger magnetic-field strength than indicated by the Planck polarisation data in the Eridu cloud and than measured in other diffuse clouds.The CR probed in our analysis should therefore not be depleted because of a lower diffusivity on the cloud boundary.
Clouds of similar structure have already been probed by Planck Collaboration Int.XXVIII (2015), with the same analysis method.In their study of the Chamaeleon complex, they identified an intermediate-velocity arc of gas which compares with the Eridu cirrus in several ways: they both have an elongated structure and comparable column densities around 5 × 10 20 cm −2 ; they exhibit large velocity gradients spanning between −40 and −4 km s −1 in the arc and −15 to −2 km s −1 in Eridu.The arc lies between 50 and 100 pc below the Galactic plane, assuming the same distances as in Planck Collaboration Int.XXVIII (2015).The arc has, however, a γ-ray emissivity that is 8% above q LIS for a spin temperature of 140 K, thereby indicating that the lower CR flux in the Eridu cirrus is not linked to the gas filamentary structure and velocity shear.
In order to assess the propagation regime in this cloud, we used the model of Evoli et al. (2018) describing the vertical CR transport out of the Galactic disc through the competition of two processes: propagation is dominated by scattering on self-generated Alfvén waves within a few kiloparsecs from the Galactic plane and at sub-TeV energies, whereas ISM turbulence injected by sources near the plane and advected vertically away takes over when it has had time (hence altitude) to cascade down to low enough scales.This model neglects CR diffusion on the toroidal component of the Galactic magnetic field.The Eridu cirrus lies at a distance of 200-250 pc from the Galactic plane and we probed its CR with energies in the range of 2-500 GeV.According to Evoli et al. (2018), we expect them to preferentially diffuse on self-generated Alfvén waves and to be advected at the Alfvén speed.The estimates of the gas density and magneticfield strength obtained in Paper I allow to compare the Alfvén speed in the Eridu cirrus ( 4 km s −1 ) and in the West Rim of the superbubble (9-14 km s −1 ).We only have an estimate of the B sky component of the magnetic field in Eridu (hence the lower speed limit given above), but having an Alfvén velocity in the cirrus much larger than in the West Rim would require an unusually large total field strength in excess of >10−20 µG for this type of diffuse cloud (Crutcher 2012).So the slower transport by Alfvén waves in the Eridu cirrus cannot explain a loss of CR.
On the other hand, dust polarisation observations with Planck show both a low dispersion in polarisation orientations and a high polarisation fraction along the cirrus, as along the outer rim of the superbubble.This is visible in Fig. 2 of Soler et al. (2018).The polarisation data implies a low level of magnetic turbulence inside the cirrus, of order δB/B 0 ∼ 0.1, which compares with the low level found along the superbubble rim due to magnetic compression in the expanding shell.A closer view of the magnetic orientation towards the Eridu cirrus is displayed in Fig. 10.Its magnetic field appears to be ordered and oriented along the filamentary structure of the cirrus, as expected for low-column density atomic clouds (Clark et al. 2014;Planck Collaboration Int. XXXII 2016).An ordered magnetic field with a low level of turbulence favours a rapid transport of CR along the field lines.
Figure 10 further illustrates the change in polarisation orientation between the cirrus direction and the superbubble West Rim on the one hand, and between the cirrus and the rest of the gas at very low negative latitudes where the field is mostly parallel to the Galactic plane.Both changes in Stokes parameters are significant at the 40 scale and they happen in regions of high signal-to-noise ratio (S/N) in polarised intensity (P/σ P ≥ 3).The dust signal from the cirrus dominates the polarisation data along those lines of sight.This data therefore suggests that the Eridu cirrus is on a different bundle of field lines.It does not follow the circular rim of the superbubble shell, nor the planeparallel configuration of the diffuse H I layer, but the field lines rather point out towards lower Galactic latitudes, with a ∼45 • angle from the south Galactic pole.The lower CR flux pervading the Eridu cirrus could then be due to propagation along a magnetic "tube" pointing towards the Galactic halo.We are currently studying another cirrus with a similar orientation towards the halo to confirm this possibility.
We also note that the older part of the superbubble corresponds to a direction in the sky of low synchrotron intensity at 408 MHz (Haslam et al. 1982).This emission probes CR electrons with energies near 10 GeV in the magnetic fields of the cirrus and superbubble shells.The LIS electron spectrum has a spectral index of −3.5 around this energy (Boschini et al. 2018).Despite the enhanced field strengths in the superbubble shell, we observe the same synchrotron intensities in the Eridu and West Rim directions for which we have magnetic field measurements (Paper I).The ratios of synchrotron intensities between the Eridu and West Rim directions indicate that the CR electron flux in the cirrus is half that in the superbubble West Rim, thereby confirming that there are fewer CR in this diffuse cloud than along the superbubble shells.

Local cosmic-ray vertical gradient
Assuming there is no strong confinement in the superbubble, one could consider the simplest 1D model of pure diffusion In order to probe the relationship between altitude and CR density, we plotted in Fig. 13 the γ-ray emissivities measured in several clouds against the altitude to the Galactic plane.The sample includes the local values found by Remy et al. (2017) in anticentre clouds.We have discarded their Perseus cloud estimate as they note that strong confusion between H I and CO structures in this compact cloud leads to an unreliable emissivity measurement in each phase.The sample also includes measurements towards high-latitude clouds seen at various distances, including intermediate-and high-velocity H I clouds (Tibaldo et al. 2015).The distance to the nearby Ursa Major cloud included in this analysis has recently been constrained to be 360 ± 20 pc, so we have updated its height in Fig. 13.The distances to the high-velocity clouds are still poorly constrained.Figure 13 shows that the emissivities found in the superbubble shells are consistent with the spread shown in local clouds from previous studies.Systematic uncertainties in H I optical depth corrections can explain the 10-15% dispersion, but the Eridu cirrus stands at variance as the lowest black data point.
The overall trend in Fig. 13 suggests an emissivity drop with Galatic height.A linear fit to the sample of local clouds (|z| ≤ 300 pc) yields a 5.5σ detection of a slope.The resulting halo size of 490 ± 60 pc (including the Eridu data point or not) is too low compared to the 5.4 ± 1.4 kpc value inferred from GALPROP models (Trotta et al. 2011).It compares favourably with the few hundred parsec heights favoured by dynamic spiral arm models, which add ISM advection downstream of the spiral arms to diffusion (Nava et al. 2017).A linear fit to the set of highlatitude data points (only) yields a halo height of 2.2 ± 1.9 kpc (Tibaldo et al. 2015).Figure 13 also shows the best-fit f (z) curve obtained for the case combining vertical diffusion and advection (dotted curve).It yields a halo height of 2.6 ± 0.1 kpc in reasonable agreement with elemental constraints, and a small scaling factor ζ = (1.7 ± 0.3) × 10 −3 .The latter implies that, for a standard diffusion coefficient of 6.2 × 10 28 cm 2 s −1 (R/10 GV) 0.344 (Jóhannesson et al. 2019), vertical advection has little impact on the local propagation of 10 GeV CR.Adding 10% systematic uncertainties on the data points because of H I optical depth corrections significantly impacts the halo height which becomes 4.5 ± 0.2 kpc, but it does not change the conclusion about a negligible advection velocity.As discussed in the previous section, the loss of CR in the Eridu cirrus may correspond to its specific magnetic configuration, and not to large-scale transport properties.Removing this cloud from the sample does not significantly change the halo height indicated by the fit.Additional γ-ray analyses of nearby clouds at an altitude of about 500 pc are clearly necessary to confirm this trend.
Figure 13 also displays the expected CR profile from the halo model of Evoli et al. (2018) for 10-GeV protons.We have set their distribution to match the γ-ray emissivity at low altitude.Their model includes advection of ISM turbulence away from the plane, as well as CR scattering on self-generated waves induced by the streaming instability through the vertical CR gradient.This model successfully reproduces the hardening of primary CR near 300 GV, but Fig. 13 shows that the modelled drop in CR flux with height is shallower than suggested by the γ-ray observations of individual clouds.

Conclusion
In order to assess the impact of the turbulent environment of superbubbles on CR propagation, we probed γ-ray emissivities in atomic clouds of the Orion-Eridanus superbubble using ten years of Fermi-LAT data.Our analysis revealed a CR spectrum that is uniform within 6% between the superbubble shell and the average found in the rest of the local ISM within a few hundreds parsecs from the Sun.The CR spectrum in the superbubble is also consistent with the measurement in the Solar system.Uncertainties from the HI optical depth correction in those shells amount to 7% and can easily explain such a small difference.We found no spectral evidence of CR (re)-acceleration even though the stellar and gas content of the superbubble suggests that it could actively re-accelerate CR.Its inefficiency could be due to a lower frequency and space density of supernovae in the last Myr, compared to other starburst regions like Cygnus X.It could also be due to a weak average magnetic field ( 1 µG) inside the superbubble or to a steep, Kolmogorov-like spectrum of the magnetic turbulence if the average field reaches a few microGauss.
We have gathered a sample of nearby clouds that have been probed in γ rays with Fermi-LAT using the same analysis method A96, page 15 of 19 A&A 635, A96 (2020) to separate the clouds along the lines of sight and to separate the different gas phases in order to measure the γ -ray emissivity in the atomic phase.Their comparison shows emissivity variations as a function of height above the Galactic plane.Adding more clouds to this sample will allow for a local estimation of the CR halo size and will establish or rule out the importance of vertical advection, rather than inferring these important parameters from larger-scale propagation models.
A cirrus cloud located outside the superbubble presents a 34% lower CR flux compared to the local ISM average.Its filamentary structure lies along magnetic-field lines that appear to be highly inclined with respect to the Galactic plane.This result opens the possibility to investigate how CR flow along magnetic flux tubes in order to better understand the role of gas fountains and of their magnetic field loops in CR escape to the halo.

Fig. 1 .
Fig. 1.H α intensity map of the Orion-Eridanus superbubble based on VTSS, SHASSA and WHAM data.The dashed line traces the perimeter of the present analysis.The white labels on the left hand side mark the key H α features towards Orion.In the analysis region, the black contours delineate the main H I shells related to the superbubble, i.e the North Rim (N), South Loop (S), East Shell (E), and West Rim (W), and the Eridu cirrus (C) lying outside the superbubble (Joubaud et al. 2019).The molecular MBM 20 cloud (M) in front of the superbubble edge is outlined in white.The line of constant Galactic latitude b = −8 • in the upper left corner indicates the orientation of the Galactic plane.

Fig. 2 .
Fig. 2. Top: γ-ray counts of gaseous origin recorded in the 0.25-63 GeV band in a 0. • 25 pixel grid.γ-ray emissions other than due to CR interactions with the gas were subtracted.Middle: dust optical depth measured at 353 GHz and displayed at the Fermi-LAT angular resolution for comparison.Bottom: dust optical depth measured at 353 GHz at its original resolution of 5 .

Fig. 3 .
Fig. 3. Photon yields arising in the γ-ray model for the best-fit to the data in the 0.25-63 GeV band on a 0. • 25 grid, from the H I and CO phases of the North Rim (EriN), South Loop (EriS), East Shell (EriE), West Rim (EriW), MBM 20, Cetus and Taurus, Eridu cirrus, Main Taurus and Galactic disc clouds, from the DNM gas column density, from the H α template (HII), from the γ-ray point sources, from the Sun and Moon emissions (SM), from the isotropic background (iso) and from the Galactic IC emission.

Fig. 4 .
Fig. 4. Residual maps, expressed as data minus best-fit model, in sigma units for the dust optical depth (σ τ , top) and the γ-ray counts in the total energy band (σ γ , bottom).The residuals have been smoothed with a 0. • 15 Gaussian kernel for display.The left-and right-hand side maps show the residuals for models without and with the DNM clouds, respectively.

Fig. 5 .
Fig. 5. Ratio of the DNM to H I emissivity factors in the North Rim as a function of energy, for a spin temperature of 100 K.The grey band gives the ±1σ confidence range on the ratio obtained for the full energy band.The constant ratio indicates that the same CR spectrum permeates the H I envelope and DNM phase.

Fig. 6 .
Fig. 6.Ratio of the CO to H I emissivity factors as a function of energy, for a spin temperature of 100 K, for two sets of molecular clouds, in the North Rim of the superbubble (dark) and in the North-Taurus-Cetus clouds (light).The bands give the ±1σ confidence range on the ratios obtained for the full energy band.The constant ratios indicate that the same CR spectra permeate the atomic and molecular phases of the clouds.

Fig. 10 .
Fig. 10.N H I column density map in Galactic coordinates.The brown contour delineates the West Rim and the South Loop of the superbubble at 2.6 ×10 20 cm −2 and the magenta countour delineates the Eridu cirrus at 3.1 × 10 20 cm −2 .The black lines give the orientation of the plane-ofthe-sky magnetic field inferred from the dust polarisation observations at 353 GHz with Planck .

Fig. 12 .
Fig. 12. Theoretical steady-state proton distribution, integrated between 10 GeV and 10 TeV, for Galactic cosmic rays injected at the edge of the turbulent superbubble (according to Tolksdorf et al. 2019).The different curves correspond to different values of the efficiency parameter ξ.The efficient case (ξ = 0.1) is in blue, the intermediate one (ξ = 1) in green and the inefficient case (ξ = 10) in orange.R SB is the superbubble radius and R Esc the free escape boundary radius.

Fig. 13 .
Fig.13.Integrated (0.4-10 GeV) γ-ray emmissivity per gas nucleon as a function of height above the Galactic plane.The black data points are from this study, the dark grey points from nearby anticentre clouds(Remy et al. 2017) and the light grey ones from high-latitude clouds (adapted fromTibaldo et al. 2015, see text).The grey band gives the ±1σ dispersion around the mean emissivity found at low Galactic height (≤300 pc).The dotted curve shows the best fit for pure CR diffusion and advection with Galactic height (see text).The dash-dotted curve corresponds to the distribution of 10 GeV protons in the halo model ofEvoli et al. (2018), rescaled to the q LIS emissivity at low Galactic height.