© The Authors 2024

1 Introduction

The Sun is located in a valley of low gas density surrounded by two walls of dense molecular clouds (see Fig. 2 of Leike et al. 2020). Part of the valley, out to ~150 pc or ~250 pc depending on the direction, corresponds to the Local Bubble (see Fig. 1 of Pelgrims et al. 2020), but the elongated valley extends further away in opposite directions roughly parallel to the local spiral arm. Interactions of cosmic rays (CRs) with the interstellar gas and radiation fields present inside and around the valley produce γ rays that have been recorded at energies above 30 MeV by the Large Area Telescope (LAT) on board the Fermi Gamma-Ray Space Telescope. As the hadronic contribution to this γ-ray intensity integrates the product of the CR and gas volume densities along the lines of sight, the γ-ray map of nearby clouds can be used to probe the CR flux in different locations and in different types of environment (Lebrun et al. 1982). This method requires an accurate measurement of the gas mass in a cloud.

In the warm atomic phase of clouds, often referred to as the warm neutral medium (WNM), one can directly infer the gas mass from observations of the HI 21 cm lines (e.g. Dickey et al. 2000; Heiles & Troland 2003). In the denser, cold atomic phase of the clouds, where the HI brightness temperature probes the excitation (spin) temperature of the gas rather than the column density, the lines must be corrected for partial self-absorption in order to estimate the gas mass (Murray et al. 2015, and references therein). This cold atomic phase, referred to as the cold neutral medium (CNM), generally contributes a fraction between 5% and 28% of the total atomic mass (Murray et al. 2018, 2020). H₂ molecules radiate inefficiently in the cold conditions of molecular clouds, so we must rely on observations of surrogate lines, such as the rotational lines of ¹²CO and its isotopologues, in order to quantify the H₂ mass. This is only possible where the abundance and excitation temperature of CO molecules are sufficient. Environmental variations of the X_CO conversion factor between the ¹²CO(1–0) line intensity and the H₂ column density have been studied theoretically and observationally (Bolatto et al. 2013; Gong et al. 2018; Remy et al. 2017, and references therein), but the variations are too large and too uncertain to reliably use the CO data to measure the CR flux in the dense molecular phase of clouds.

Complementary information can be obtained from the thermal emission of large dust grains mixed with the gas. The emission is observed at submillimetre and infrared wavelengths (Planck Collaboration XI 2014; Planck Collaboration Int. XVII 2014). Dust grains and CRs are both tracers of the total gas in all chemical and thermodynamical forms (Grenier et al. 2005; Planck Collaboration Int. XXVIII 2015). The γ rays reveal the gas mass convolved by the ambient CR flux, whereas the dust optical depth shows the gas mass convolved by the ambient dust-to-gas mass ratio and by the ambient grain emissivity. Both the CR and dust properties vary from place to place, but the spatial correlation between the two tracers in a given cloud can be used to trace the dark gas that escapes HI and CO line observations (Grenier et al. 2005). This dark neutral medium (DNM) accumulates at the HI–H₂ transition in the clouds (Liszt et al. 2019). It is not seen in HI or CO emission because of self-absorption in the dense HI clumps and because CO molecules are efficiently photo-dissociated or too weakly excited in the diffuse H₂ gas.

Given the uncertainty in gas mass estimates in the dark and molecular phases of clouds, one can only infer the CR spectrum with adequate precision in the atomic phase of the interstellar medium (ISM). The derivation of the γ-ray emissivity, q_HI, per gas nucleon in the atomic gas relies on the spatial correlation between the γ-ray intensity and HI column density (N_HI), but the analysis must include maps of the other gas phases in order to avoid biasing the HI emissivity spectrum (Planck Collaboration Int. XXVIII 2015). These spectra are typically obtained at energies above a few hundred MeV in order to take advantage of the better angular resolution of the LAT in order to resolve the clouds and their different phases.

Such emissivity measurements have been performed in different clouds of the local ISM (e.g. Abdo et al. 2010; Ackermann et al. 2012a; Planck Collaboration Int. XXVIII 2015; Casandjian 2015; Tibaldo et al. 2015; Remy et al. 2017; Joubaud et al. 2019) and across the Milky Way (e.g. Ackermann et al. 2013; Acero et al. 2016). The local emissivity spectra appear to be often consistent with the local average, q_LIS(E), which has been obtained for the whole atomic gas seen at Galactic latitudes ranging between 7° and 70° (Casandjian 2015). This consistency implies that the CR flux is rather uniform within a few hundred parsecs of the Sun (Grenier et al. 2015).

At variance with this uniformity, a significantly lower γ-ray emissivity spectrum has been found in a nearby cloud named Eridu. The spectrum is 34% lower than the q_LIS average (Joubaud et al. 2020). The deviation is statistically significant (14 σ) and there is no confusion with other clouds along the lines of sight. The deviation cannot be attributed to uncertainties in the HI mass because the minimum amount of gas inferred from the HI lines (optically thin case) gives an upper limit to the emissivity that is still 25% below q_LIS. More CRs would outshine the observed γ rays. The γ-ray analysis has probed CRs with energies above several GeV. Their penetration in the cloud is likely unhampered given the modest column densities (1–5 × 10²⁰ cm⁻²) and low volume densities (around 7 cm⁻³) of the diffuse gas (Joubaud et al. 2019). The Eridu cloud therefore appears to be pervaded by a ~30% lower CR flux, but with the same particle energy distribution as in the local ISM.

The Eridu cloud is a diffuse, atomic, and filamentary cirrus. It is located outside the edge of the Orion-Eridanus superbub-ble, at high Galactic latitudes of around 50° and at an altitude of about 200–250 pc below the Galactic plane (Joubaud et al. 2019). Figure 1 shows that the filament is inclined with respect to the Galactic plane and the dust polarisation data at 353 GHz from Planck (Planck Collaboration Int. XIX 2015) show that the filament extends along a bundle of well-ordered magnetic-field lines, pointing within 45° of the southern Galactic pole.

Magnetic loops are pushed out of the Galactic disc by gas fountains. In order to investigate their impact on CR transport to the halo, we searched for another filament with gaseous and magnetic characteristics comparable to those of Eridu. The Reticulum cloud is also an atomic and filamentary cloud that lies at the edge of the Local Valley, at a distance we estimate to be about 240 pc (Appendix A). Figure 1 shows that, like Eridu, it is largely inclined with respect to the Galactic plane. It points towards the southern Galactic halo and has well-ordered magnetic field lines aligned with the filament axis. We therefore aim to measure the γ-ray emissivity spectrum per gas nucleon in the Reticulum filament and to compare it with the average in the local ISM and with the emissivity spectrum in the Eridu cloud.

The paper is structured as follows: The data are presented in Sect. 2 and the models are described in Sect. 3. The results and the detection of dark gas are detailed in Sect. 4. The CR content of the clouds is discussed in Sects. 5.1–5.5.

Fig. 1 Dust optical depth τ353 at 353 GHz in Galactic coordinates towards the Eridu (top) and Reticulum (bottom) filaments. The black lines give the orientation of the plane-of-the-sky Bsky magnetic field inferred from the dust polarisation observations at 353 GHz and at an angular resolution of 14′ (FWHM). The red contours delineate the filaments at HI gas column densities of 3 × 1020 cm−2 and 1.2 × 1020 cm−2, respectively.

2 Data

We selected an analysis region around the Reticulum filament extending from 251°. 5 to 302°. 5 in Galactic longitude and from −54°.5 to −30°. 5 in Galactic latitude, as shown in Fig. 2. We masked regions with large column densities from the Large Magellanic Cloud (LMC) and Small Magellanic Cloud (SMC), and from the edge of the Gum Nebula.

2.1 Dust data

In order to trace the dust column density, we used the all-sky map of the dust optical depth, τ₃₅₃, at 353 GHz (Planck Collaboration XI 2014, updated by Irfan et al. 2019). It was constructed at an angular resolution of 5′ by modelling the intensity, I_ν, of the dust thermal emission recorded in frequency v by Planck and the Infrared Astronomical Satellite (IRAS). The dust optical depth relates to the modified black-body intensity B_v(T_d) for a dust temperature T_d and spectral index, β, according to I_v = τ₃₅₃ B_v(T_d)(ν/353 GHz)^β. The dust opacity at 353 GHz is defined for a hydrogen column density N_H as the ratio τ₃₅₃/N_H.

Figure 2 shows the distribution of the dust optical depth in the analysis region at the original Planck resolution in the upper panel and at the angular resolution of Fermi-LAT for the 0.16– 63 GeV band in the middle one.

Fig. 2 Dust optical depth measured at 353 GHz from the Planck and IRAS data and displayed on a 10′.2×7′.5 grid at the angular resolution of Planck (top panel) and of Fermi-LAT (middle). The bottom panel shows the γ-ray photon counts produced by CR interactions with gas and recorded by Fermi-LAT in 0°25×0°25 bins in the 0.16–63 GeV energy band (other emissions of non-gaseous origin are subtracted). Three zones are excluded from the analysis region: the SMC on the left, LMC in the middle, and the edge of the Gum nebula in the upper right. The Reticulum cloud is the long filament crossing the region below the LMC.

2.2 Gamma-ray data

We used 14 yr of Fermi-LAT Pass 8 P8R3 survey data for energies between 0.16 and 63 GeV (Bruel et al. 2018) and we analysed the data in eight energy bands that are bounded by 10^2.2, 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.

The Fermi-LAT point spread function (PSF) strongly varies with energy. In order to retain sufficient photon statistics in this faint region and to preserve a good contrast in photon counts across small cloud structures, we combined photons of different PSF types in each energy band. Photons are labelled from the worst (PSF0) to best (PSF3) quality in angle reconstruction. We used PSF 3 in the three lowest-energy bands, PSF 2 and 3 in the next band, PSF 1 to 3 in the next two bands, and PSF 0 to 3 in the two highest energy bands. The 68% containment radius of the resulting PSF varies from 0°.9 above 630 MeV to 1°.4 at 250 MeV and to 2°.1 at 160 MeV. The resulting photon map is shown in the lower panel of Fig. 2 where we subtracted the point sources and the very smooth isotropic and inverse-Compton (IC) components to display the γ rays produced in the clouds.

In order to reduce to less than 10% the contamination of the photon data by residual CRs and by Earth atmospheric γ rays, we applied tight selection criteria such as the SOURCE class and photon arrival directions within an energy-dependent and PSF-type-dependent angle varying from 90° to 105° from the Earth zenith, namely 90° in the first energy band, 100° in the second, and 105° at higher energies. We used the associated instrument response functions P8R3_SOURCE_V3.

We modelled the interstellar emission in a broader region than the analysis one in order to take into account the extent of the Fermi-LAT PSF. The 5°-wide peripheral band around the analysis perimeter corresponds to the PSF >95% containment radius for our selections. It is displayed in Fig. 3.

The positions and flux spectra of the γ-ray sources in the field were provided by the Fermi-LAT 4FGL-DR3 12-yr source catalogue (Abdollahi et al. 2022). The analysis region and the 5°-wide peripheral band contain 220 point sources. Their flux spectra were calculated with the spectral characteristics given in the catalogue.

The observation also includes IC γ rays produced by CR electrons up-scattering the interstellar radiation field in the Milky Way. This was modelled by GALPROP¹ and we used the LRYusifovXCO5z6R30_Ts150_Rs8 model version (Ackermann et al. 2012b). The isotropic spectrum for the extragalactic and residual instrumental backgrounds is given by the Fermi-LAT science support centre². We used the isotropic spectrum provided for each PSF type of the SOURCE class.

2.3 Gas data

In order to map the neutral atomic hydrogen (HI) in the region, we used the 21-cm Galactic All Sky Survey GASS III data (Kalberla & Haud 2015). The survey has an angular resolution of 16′.2, a velocity resolution of 0.82 km s⁻¹ in the local standard of rest, and an rms sensitivity of 57 mK. Lines in the analysis region span velocities from −125 to +495 km s⁻¹.

Following the method developed by Planck Collaboration Int. XXVIII (2015), we decomposed the HI spectra into sets of individual lines with Pseudo-Voigt profiles. Based on the (l, b, v) distribution of the fitted lines, we identified seven nearby cloud complexes that are coherent in velocity and in position across the region. We named these clouds according to the constellation in their direction, namely Reticulum (Ret), Far Hydrus (HyiF), North Hydrus (HyiN), South Hydrus (HyiS), East Dorado (DorE), West Dorado (DorW), and Horologium (Hor). They are displayed in Fig. 3 and they span velocities from 3.3 to 20.6 km s⁻¹, −15.7 to −1.6 km s⁻¹, −0.8 to 13.2 km s⁻¹, −14.8 to 2.5 km s⁻¹, 2.5 to 4.9 km s⁻¹, 5.8 to 33.8 km s⁻¹, and −14.8 to 5.8 km s⁻¹, respectively. These clouds are local. We gathered lines from the Magellanic clouds in a single map. We added a ninth component gathering the very weak lines found at velocities below −8 km s⁻¹.

We produced N_HI column-density maps for each complex by integrating the fitted lines that have a central velocity within the chosen velocity interval for each (l, b) direction. This method corrects for potential line spill-over from one velocity interval (i.e. cloud) to the next. The line fits produce small residuals between the modelled and observed HI spectra. In order to preserve the total HI intensity, we distributed these residuals among the fitted lines according to their relative strength in each velocity channel.

In order to correct the HI intensity for potential self-absorption, we produced the N_HI maps for a set of uniform spin temperatures, T_s, of 100, 125, 150, 200, 300, 400, 500, 600, 700, and 800 K and also for the optically thin case. Nguyen et al. (2018) have shown that a simple isothermal correction of the emission spectra with a uniform spin temperature can well reproduce the more precise N_HI column densities that are inferred from the combination of emission and absorption spectra. Given the wide range of potential spin temperatures we could not test individual temperatures for each cloud, but we used the same temperature for all clouds and we let it vary when modelling the dust and γ-ray data.

The nearby clouds in the region are mainly composed of atomic gas. We looked for CO(1–0) emission at 115 GHz towards them, but did not find published observations in that area. We checked that Planck detected no CO(1–0) intensity up to 1 K km s⁻¹ in the region (Planck Collaboration XIII 2014).

Fig. 3 NHI column densities in the Reticulum (Ret), Far Hydras (HyiF), North Hydras (HyiN), South Hydras (HyiS), East Dorado (DorE), West Dorado (DorW), Horologium (Hor), and the Magellanic clouds for optically thin HI. The faint column densities in the negative-velocity map (NegVel) are multiplied by 4 for display. The maps include a 5º-wide peripheral band around the analysis region (black perimeter) to account for the wings of the Fermi-LAT PSF in the γ-ray model.

2.4 Cloud distances

Leike et al. (2020) used the stellar data from the Gaia-DR2, ALLWISE, PANSTARRS, and 2MASS surveys to model the dust extinction and to reconstruct the 3D map of dust clouds within ~400 pc from the Sun with a spatial resolution of 2 pc. Distances to compact clouds and their dust densities are expected to be well constrained by the precise locations of the stars. In Figs. A.1 and A.3, we display cuts of the 3D dust distribution in 20 pc wide intervals, spanning distances from 150 pc to about 350 pc. The shapes and locations of the clouds isolated in the HI data in Fig. 3 can be roughly recognised in the dust distribution. We also probed the dust distribution along lines of sight exhibiting large gas column densities in the different clouds (see Fig. A.2). From the dust and gas correspondence we infer average distances of 240 ± 30 pc, 220 ± 30 pc, 220 ± 30 pc, 170 ± 30 pc, 310 ± 30 pc, 340 ± 30 pc, and 280 ± 30 pc to the Reticulum, North Hydrus, South Hydrus, Far Hydrus, East Dorado, West Dorado, and Horologium clouds, respectively. The ± 30 pc uncertainty reflects the size of dust overdensities that can be isolated in the 3D dust cube.

We applied the same method to revisit the distance to the Eridu cloud. The rim of the Orion-Eridanus superbubble is clearly detected at a distance around 200 pc in the dust 3D data. The Eridu cloud lies beyond this edge. Figure A.3 shows a second dust front at a distance of 300 ± 30 pc with an inclination and extent in the sky that is compatible with that of the elongated Eridu cloud seen in HI. This distance places the cloud at an altitude of 260 pc below the Galactic plane, near the edge of the ~220 pc thick HI layer of the local ISM (Nakanishi & Sofue 2016). A much larger distance would be disfavoured. The HI line velocities of the cloud are also typical of the local ISM and not of intermediate- or high-velocity clouds outside the disc. Extending the 3D dust analysis to larger distances than available in the Leike et al. (2020) cube would help confirm the distance to Eridu.

3 Models and analyses

Figure 2 shows that the dust optical depth and the γ-ray intensity produced by CRs in the different clouds are strongly correlated. The dust and γ rays both trace, to first order, the total gas column density N_H. In order to detect neutral gas not observed in the HI data, we used the fact that it is permeated by both CRs and dust grains. We therefore extracted γ-ray and dust residuals between the observation and the model expected from the N_HI distribution. We then used the spatial correlation between the positive residuals found in dust and in γ rays to extract the additional gas column densities (see Sect. 3.3). This additional gas corresponds to the DNM at the HI–H₂ transition where the dense HI becomes optically thick to 21-cm line emission and CO emission is absent or too faint to trace the H₂ gas. The analysis results yield small amounts of dark gas in the Reticulum and North Hydrus clouds (see Sect. 4.2).

3.1 Dust model

Provided that the dust-to-gas mass ratio and the mass emission coefficient of the grains are uniform (but not necessarily the same) in each gas phase of a particular cloud, the dust optical depth τ₃₅₃(l, b) can be modelled as a linear combination of the different gas column densities , each with free y normalisation to be fitted to the data, as in Planck Collaboration Int. XXVIII (2015). A free isotropic term, y_iso, was added to the model to account for the residual noise and the uncertainty in the zero level of the dust maps (Planck Collaboration XI 2014; Casandjian et al. 2022).

The τ₃₅₃(1, b) model can be expressed as: (1)

where N_HI(l, b) denotes the nine N_HI maps depicted in Fig. 3 and stands for the column density in the two DNM components constructed from the γ-ray data (see Sect. 3.3). We degraded the optical-depth map from its original resolution of 5′ down to 16′.2 to match that of the GASS HI survey.

The y_HI,i coefficients give the average values of the dust τ₃₅₃/N_H opacity in the different clouds. The y_DNM,i parameters can probe opacity changes in the denser DNM phases. They were fitted to the data using a least-square χ² minimisation. We expect the uncertainties of the model to exceed those of the dust map because of our assumption of uniform grain distributions through the clouds and because of the limited capability of the HI data to trace the total gas. As we cannot precisely estimate the model uncertainties, we set a fractional error of 11% that yields a reduced χ² value of one in the dust fit. This fraction is slightly larger than the 1% to 9% measurement uncertainties in the τ₃₅₃ values across this region.

3.2 Gamma-ray model

Earlier studies have indicated that the Galactic CRs radiating at 0.1–100 GeV have diffusion lengths much larger than typical cloud dimensions and that they propagate through the different gas phases (Grenier et al. 2015). The observed interstellar γ-ray emission can therefore be modelled, to first order, by a linear combination of the gas column densities present in the different clouds seen along the lines of sight, in both the HI and DNM gas phases. The model also includes a contribution from the Galactic IC intensity, I_IC(l, b, E), point sources of non-interstellar origin with individual flux spectra S_j(E), and an isotropic intensity, I_iso(E), accounting for the extragalactic γ-ray background and for residual CRs misclassified as γ rays. We verified that the soft emission from the Earth limb is not detected in the present energy range for the choice of maximum zenith angles.

The γ-ray intensity, I(l, b, E), in γ cm⁻² s⁻¹ sr⁻¹ MeV⁻¹ in each (l, b) direction, is modelled in each energy band E as: (2)

where N_HI,i(l, b) denotes the HI column density map of each cloud and the extraction of the dust DNM maps is explained in Sect. 3.3. The q coefficients of the model are determined from fits to the Fermi-LAT data.

The q_HI,i, q_IC, and q_iso parameters are simple normalisation factors to account for possible deviations from the input spectra taken for the local gas γ-ray emissivity q_LIS(E) and for the isotropic and IC intensities. The input q_LIS spectrum was based on the correlation found between the γ radiation and the N_HI column densities derived from the Leiden-Argentine-Bonn (LAB) survey for a spin temperature of 140 K, at latitudes 10° ≤ |b| ≤ 70° (Casandjian 2015). The q_HI,i parameter in the model can therefore compensate for differences in the HI data (calibration, angular resolution, spin temperature) and for cloud-to-cloud variations in CR flux. Such differences will affect the normalisation equally in all energy bands whereas a change in CR penetration inside a specific cloud will show as an energy-dependent correction, primarily as a loss at energies below one GeV (Schlickeiser et al. 2016).

The source flux spectra Sj(E) were calculated with the spectral characteristics given in the 4FGL-DR3 catalogue. The q_{S j} (E) normalisation factors in the model allow for possible changes due to the longer exposure used here and to the use of a different interstellar background for source detection in the 4FGL-DR3 catalogue. The q_{S j} (E) normalisation factors of the sources inside the analysis region or in the peripheral band were fitted individually, except for those inside the SMC and LMC regions and for those found in the lowest dust zones (τ₃₅₃ < 3 × 10⁻⁶), which were treated collectively in three separate sets with a single normalisation factor per set.

In order to compare the model 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 and energy resolution. The convolution with the response functions was done for each PSF event type independently in small energy bands. We also calculated the isotropic map for each PSF event type in the same small energy bands. We then summed the photon count maps in energy and for the different combinations of PSF types to get the modelled count maps in the eight energy bands of the analysis for the same photon selection as in the data.

Figure 4 gives the photon yields that were obtained in the entire energy band. It shows that the emission originating from the gas can be distinguished from the much smoother emission backgrounds and that the LAT angular resolution is sufficient to separate the various clouds, as well as the different gas phases within the clouds. The distinct spatial distributions of the clouds and of their phases allow rather independent measurements of their individual γ-ray emissivities. On the contrary, we do not expect a firm separation of the isotropic and IC spectra across this region. We used a binned maximum-likelihood with Poisson statistics to fit the model coefficients q to the LAT data in each of the eight energy bands and in the overall one.

Fig. 4 Photon yields on a 0º25 pixel grid from the main components of the γ-ray model in the 0.16–63 GeV band. The yields come from the HI and DNM column densities derived for the optically thin case in the different clouds, and from the isotropic background, Galactic IC emission, and γ-ray point sources.

3.3 Analysis iterations

In order to extract the DNM gas present in both the dust and γ-ray data, we iteratively coupled the γ-ray and dust models. We started by fitting the dust optical depth and the γ-ray data independently according to Eqs. (1) and (2) with only the HI maps as gaseous components and with ancillary (other than gas) components. We then built maps of positive residuals (see how below) between the data (dust or γ rays) and the best-fit contributions of the model components. We separated these excesses into two DNM maps according to their location, one in the Reticulum cloud and the other one in the North Hydrus area.

The next loops in the iterative process provided the DNM templates estimated from the dust to the γ-ray model ( in Eq. (2)) and conversely provided the DNM column density maps derived from the γ rays to the dust model in Eq. (1)). We repeated the iteration between the γ-ray and dust models until the log-likelihood of the fit to γ-ray data saturated.

The estimates of the q and y model coefficients and of the DNM templates improve at each iteration since there is less and less need for the other components, in particular the HI ones, to compensate for the missing gas. They can still miss gas at some level because the DNM maps provided by the γ rays or dust emission have limitations (e.g. dust emissivity variations, limited γ-ray sensitivity).

Care must be taken in the extraction of the positive residuals because of the noise around zero. A simple cut at zero would induce a positive offset bias, and so we denoised the residual maps using the multi-resolution support method implemented in the MR filter software (Starck & Pierre 1998). For the dust residuals, we used a multi-resolution thresholding filter, with six scales in the B-spline-wavelet transform (à trous algorithm), Gaussian noise and a 2-σ threshold to filter. For the γ-ray residuals, we used the photon residual map obtained in the 0.16–63 GeV band and we filtered it with six wavelet scales, Poisson noise, and a 2-σ threshold. The stability of the iterative analysis and the results of the dust and γ-ray fits are discussed in the following section.

4 Results

Before presenting the γ-ray emissivity spectra obtained in the Reticulum cloud in Sect. 4.4, we discuss how the fits to the γ-ray and dust data improve with optical depth corrections to the N_HI column densities (Sect. 4.1) and with the addition of DNM gas templates (Sect. 4.2). We also assess the robustness of the best-fit models with jackknife tests in Sect. 4.3.

4.1 HI optical depth correction

We do not know the exact level of HI opacity in the different clouds. This opacity induces systematic uncertainties in the N_HI column densities, and therefore also in the HI contributions to the γ-ray and dust models. Nguyen et al. (2019) showed that a uniform isothermal correction of the HI emission spectra across a cloud provides a good approximation to more precise estimates based on emission and absorption HI spectra. We therefore produced the dust and γ-ray models for a wide range of optical depth corrections, from a rather optically thick case (T_S = 100 K) to the optically thin case that yields the minimum amount of gas in the cloud. Given the moderate N_HI column densities of the present clouds, the difference between those extreme corrections reaches 25% on a few pixels and less than 1% on the mean column density of a cloud.

Even though the spin temperature is uniform, the HI opacity correction depends on the HI intensity recorded in each direction and velocity channel, so it varies across the cloud and slightly modifies the N_HI spatial distribution of a cloud. Figure 5 shows that the goodness of the γ-ray fit (after the DNM iteration) significantly improves as we increase the HI spin temperature. The N_HI maps best match the spatial structure of the γ-ray flux emerging from the atomic gas for T_S > 300 K; therefore for small HI opacities. Given the saturation of the log-likelihood ratio towards the optically thin case, we present the analysis results obtained for this thin case in the rest of the paper, unless otherwise mentioned. The residual maps between the γ-ray data and best-fit model show that using the same T_S for all the clouds can reproduce the γ-ray data without any significantly deviant cloud.

As the spin temperature T_S increases, the optical depth correction and N_HI column densities decrease, and the γ-ray emissivity q_HI and dust opacity y_HI of the HI clouds increase. They do so respectively by 4–11% and 3–7% as the spin temperature rises over the entire range. HI optical depth corrections have therefore a small effect in this region because the clouds are not massive, so they are fairly transparent to HI line radiation.

Fig. 5 Evolution of the log-likelihood ratio of the γ-ray fit in the 0.16– 63 GeV band as a function of the HI spin temperature used to calculate the NHI column densities in the clouds.

4.2 DNM detection

In order to check for the presence of DNM gas not traced by HI line intensities, we first performed the dust and γ-ray fits with only the N_HI maps as gaseous components. The residuals are displayed in the upper panels of Fig. 6. Positive residuals indicate that the data exceed the model predictions. Coincident excesses are seen in dust and in γ rays even though the two fits are independent. They cross the region in diagonal, from the North Hydrus cloud (upper left corner) down along the Reticulum filament. A coincident increase in CR density and in dust opacity is very improbable. Hadronic CR interactions with dust grains produce undetectable γ-ray emission (Grenier et al. 2005), so this correlated excess is likely due to additional gas. The distribution of positive γ-ray excesses along the Reticulum cloud is sparse, but the large extent of negative (blue) excesses across the whole analysis region indicates that the γ-ray fit has unduly pushed the HI emissivities of the clouds upward to compensate for some missing gas.

The DNM gas column density is built by correlating and filtering the dust and γ-ray excess signals in the iteration between both model fits (see Sect. 3.3). We built separate DNM maps for the North Hydrus and Reticulum clouds. Adding them improves the fits quality (see the middle maps of Fig. 6). A log-likelihood ratio increase by 161 between the γ-ray fits with and without the DNM templates implies a very significant DNM detection. The spectrum of the γ-ray emission associated with both DNM components is consistent with that found in the HI gas phase over the whole energy range (see below). This gives further support to a gaseous origin of the γ-ray emission in the DNM structures.

Positive residuals remain in the final dust fit, even after the addition of the DNM maps. Those in spatial coincidence with the DNM phase are partially due to the non-linear increase in dust opacity as the gas density rises towards the dense molecular phase (Remy et al. 2017). They also result from the limited γ-ray photon statistics that limit the accuracy of the DNM template provided to the dust fit. The negative and positive dust residuals seen at longitudes l < 270° reflect limitations in the τ₃₅₃ data map due to the difficult separation of the dust emission from other submillimetre emission in these very faint zones (Planck Collaboration XI2014; Irfan et al. 2019, and references therein). These regions were excluded from the DNM analysis.

4.3 Best fits and jackknife tests

The best dust and γ-ray fits with the models described by Eqs. (1) and (2) include the nine HI and two DNM maps for the gaseous components. Tables 1 and 2 give the best-fit coefficients obtained for the HI and DNM phases of the main Reticulum and North Hydrus clouds. The results found for the other small local clouds will be discussed in a subsequent paper (their emissivity spectra are rather close to the q_LIS expectation). We focus here on the Reticulum filament to compare it with the Eridu one and we include the emissivity measurement in the North Hydrus cloud as it partially overlaps the northern end of Reticulum. The γ-ray emissivities are given in each of the eight energy bands and for the whole 0.16–63 GeV band.

The final iteration of the dust and γ-ray fits leads to the residual maps presented in the middle panels of Fig. 6. The γ-ray residuals are small and randomly distributed, so our linear model provides an excellent statistical fit to the γ-ray data in the 0.16–63 GeV energy band across the whole region. The same conclusion is reached in the eight separate energy bands that are not shown here.

The statistical errors on the best-fit coefficients are inferred from the information matrix of the fits (Strong 1985). They include the correlation between parameters. They are generally small (2–9% for the HI gas γ-ray emissivities, 15% and 10% for DNM emissivities for Ret and HyiN, 13% for q_IC, 9% for q_iso, 0.4–1.3% for the HI dust opacities, 8% and 3% for y_DNM,Ret and y_DNM,RyiN and 2% for y_iso) because the clouds are spatially well resolved in the analysis and because tight correlations exist between the gas, dust and γ-ray maps (see Fig. 2).

We also checked the magnitude of potential systematic uncertainties due to the linear-combination approximation underlying the models. In order to evaluate the impact of possible spatial changes in the model parameters and/or in the mean level of HI self- absorption, we repeated the last iteration of the dust and γ-ray fits a thousand times over random subsets of the analysis region, each time masking out 20% of the pixels with 1°-wide randomly selected squares. In γ rays, the jackknife tests were performed for the total 0.16–63 GeV energy band and for each of the eight bands. Figure 7 shows the distribution obtained for the best-fit coefficients relating to the Reticulum and North Hydrus clouds and to the large-scale IC and isotropic components. The Gaussian-like distributions indicate that the emissivities and opacities presented in Tables 1 and 2 are statistically stable and that they are not driven by subset regions in each cloud complex. The standard deviations found in the jackknife distributions amount to 1–6% for the HI gas γ-ray emissivities, 9% and 4% for the Ret and HyiN DNM emissivities, 8% for q_IC, 5% for q_iso, 0.6–3% for the HI dust opacities, 12% and 3% for y_DNM,Ret and y_DNM,HyiN and 4% for y_iso. Both statistical errors (1σ) and systematic uncertainties (jackknife standard deviations) are given in Tables 1 and 2. We quadratically added them to construct the final uncertainties on the q and y coefficients.

Figure 8 shows the stability with γ-ray energy of the best-fit coefficients relating to the Reticulum and North Hydrus clouds and to the IC and isotropic intensities. We find no significant deviation with respect to the average coefficient fitted for the entire energy band (dashed lines). This implies that the HI and DNM gas emissivity spectra in the two clouds have the same shape as the average q_LIS (E) spectrum measured in the local ISM. We further discuss their relative normalisation with respect to q_LIS in the next section.

Figure 8 and the value of the relative intensity ratio q_IC = 1.18 ± 0.15 in the whole energy band further show that the IC spectrum detected in these directions of the sky is fully consistent in shape and in intensity with the GALPROP prediction based on the complex evolution of the interstellar radiation field and of the CR lepton spectra across the Milky Way. The isotropic spectrum found in this small region of the sky is also close to the larger-scale estimate obtained between 10º and 60º in latitude³. The relative intensity ratio is q_iso = 0.84 ± 0.08 in the 0.16–63 GeV energy band. The interstellar foregrounds are faint enough to separate the IC and isotropic components despite rather similar spatial distributions across the analysis region (see Fig. 3).

Planck Collaboration XI (2014) and Remy et al. (2017) found marked gradients in dust opacity across the sky. These variations reflect a systematic opacity increase as the gas becomes denser and the grains evolve in size and chemistry (Remy et al. 2017; Ysard et al. 2019, and references therein). Diffuse atomic cirri at high latitudes exhibit low opacities around a mean value of (7.0 ± 2.0) × 10⁻²⁷ cm² (Planck Collaboration XI 2014) or (7.1 ± 0.6) × 10⁻²⁷ cm² (Planck Collaboration Int. XVII 2014) whereas larger values ranging from 10 ×10⁻²⁷ to 16 ×10⁻²⁷ cm² are found in the denser atomic envelopes of nearby clouds (Remy et al. 2017), and values from 15 ×10⁻²⁷ to 17 ×10⁻²⁷ cm² characterise the DNM phase (Planck Collaboration Int. XXVIII 2015; Remy et al. 2017; Joubaud et al. 2020). Dust opacities further increase in the dense, CO-bright, molecular phase, reaching values as high as 60 × 10⁻²⁷ cm² (Remy et al. 2017). The dust opacities we found in the HI phase of the Reticulum and North Hydrus clouds are consistent with this trend. They are close to the low diffuse cirrus values and, in particular, to the dust opacity of (6.4 ± 0.1) × 10⁻²⁷ cm² measured in the Eridu filament (Joubaud et al. 2020). The dust opacity we found in the DNM phase of the North Hydrus cloud fully agrees with other estimates. We caution the reader that the unusually low y_DNM value found for the Reticulum cloud may not be reliable because of the difficulty discussed in Sect. 4.2 to fully extract the column density along Reticulum from the γ-ray data. The final dust residual map in Fig. 6 shows a leftover positive excess along the Reticulum filament. Taking this excess into account would increase the y_DNM value.

Fig. 6 Residual maps (data minus best-fit model, in sigma units) in γ-ray counts in the 0.16–63 GeV band (σγ) and in dust optical depth , obtained when including different sets of gaseous components in the models: without (top) and with (middle) the DNM components in addition to the HI templates. The column-density maps (bottom) extracted from the γ-ray residuals (left) and from dust residuals (right) are shown in the bottom row. The models are derived for an HI at the optically thin case.

Fig. 7 Distributions of the best-fit coefficients obtained in the thousand jackknife fits to the 0.16-63 GeV γ-ray data (two left columns) and dust data (two right columns) for the optically thin HI case. The vertical lines and overlaid numbers give the mean values of the distributions. The units are 1026 cm–2 for qDNM, 10–27 cm2 for yHI and yDNM, and 10–7 for yiso. qHI,i, qIC and qiso are dimensionless normalisation factors.

4.4 Cosmic rays in the Reticulum and Eridu filaments

Figure 9 compares the γ-ray emissivity spectrum we found per gas atom in the Reticulum filament to that in the similar filament of Eridu and in the more massive North Hydrus cloud. We also compare it to the average q_LIS spectrum measured in the local atomic gas at latitudes between 10º and 70º across the sky (Casandjian 2015). We compare them for the spin temperatures that best fitted their respective γ-ray data; therefore at T_S = 100 K for Eridu, 140 K for q_LIS, and for optically thin HI for Reticulum and North Hydrus. The results indicate that the CR population that permeates the HI and DNM phases of the Reticulum and North Hydrus clouds has the same energy distribution and same flux as the average population in the local ISM. This population applies to the 100–150 pc long sector subtended by the Reticulum and North Hydrus clouds along the rim of the Local Valley. As expected for the low N_HI column densities characterising these clouds, we find no significant spectral deviation at low energies that would mark a loss or concentration of CRs when penetrating into the denser DNM clumps (Schlickeiser et al. 2016; Bustard & Zweibel 2021).

Figure 9 illustrates that the CR flux in the Reticulum filament is 1.57 ± 0.09 times that pervading the similar Eridu filament. This discrepancy decreases to 1.51 ± 0.09 if we compare them for the same HI spin temperature of 100 K, to 1.38 ± 0.08 if we assume them both optically thin. The minimum ratio is 1.32 ± 0.08 if Eridu is optically thin and Reticulum thicker with T_s = 100 K, a case yielding a worse γ-ray fit to the data for both clouds. The loss mechanism that affects CRs in Eridu should be energy independent to preserve the same emissivity spectral shape in the three clouds across the entire energy band.

H_α hydrogen recombination line emission has been recorded along the Reticulum cloud in the Southern Halpha Sky Survey Atlas (Gaustad et al. 2001), suggesting that the cloud be a cooling filament or cooling shell surrounded by recombin-ing warm ionised gas. The spatial correlation between the dust and γ rays detects any gas in addition to the HI gas, regardless of its chemical or ionisation state. The inferred DNM column densities are dominated by the dense CNM and diffuse H₂, but they also include the much more dilute ionised gas. The γ-ray emissivity we obtained in the HI phase of Reticulum is therefore not strongly biased by untraced ionised gas surrounding it. As a further check, we used the H_α map reprocessed by Finkbeiner (2003) as an additional gas template into the dust and γ-ray fits, even though the H_α emission measure does not scale with the column density, but it integrates the product of the electron and ion densities along the line of sight, . With this crude additional template, the Reticulum γ-ray emissivity remains consistent with the local-ISM q_LIS average and is still 1.51 ± 0.11 times larger than in Eridu over the whole energy band for the choice of best fitting spin temperatures.

Assuming that the same CR flux pervades the recombin-ing and neutral hydrogen, the γ-ray signal related to the H_α template yields column densities around 10¹⁹ cm⁻² in Reticulum, an order of magnitude below those in the neutral phase. Column densities of the Galactic warm ionised medium, integrated to the Reticulum distance, are also an order of magnitude below the N_HI ones and they would not spatially correlate with the HI clouds to impact the values of the q_HI emissivities (Gaensler et al. 2008).

Fig. 8 Spectral evolution relative to the local interstellar spectrum qLIS of the γ-ray emissivities found for the HI and DNM phases of the Reticulum and North Hydrus clouds (top and middle panels). Spectral evolution of the isotropic and IC intensities found in the region relative to the input spectra (bottom panel), qDNM/qLIS is given in 1026 cm−2. The results are shown for the optically thin HI case. The error bars combine the statistical and systematic uncertainties.

Fig. 9 Spectral evolution of the γ-ray emissivities per gas nucleon in the Eridu (green), Reticulum (red), and North Hydrus (cyan) clouds. The orange curve gives the average qLIS emissivity measured in the atomic gas at latitudes 10° ≤ |b| ≤ 70° (Casandjian 2015). The emissivities are shown for the HI optical-depth correction that best fits the γ-ray data (Ts = 100 K for Eridu, 140 K for qLIS, and optically thin for the others). The shaded bands include the statistical and systematic uncertainties.

5 Discussion

The Reticulum and Eridu clouds are both light and filamentary neutral clouds, with ordered magnetic fields on a few parsec scale, and a large inclination with respect to the Galactic plane. Yet, the first has 57 ± 9% more CRs than the second. In the next sections, we try to characterise the dynamical and magnetic states of those clouds at a parsec scale in order to search for differences that can influence CR transport inside them at the much smaller scales of their gyro-motion. CR momenta parallel to the local magnetic field randomly vary by interaction with magnetic perturbations. The Alfvén and slow modes of interstellar magnetohydrodynamic (MHD) turbulence are thought to be inefficient scatterers because of their large anisotropy at small scales (Chandran 2000; Xu et al. 2016). We therefore first focus on gas properties that can impact CR diffusion in the self-confinement scenario where the particles scatter on Alfvén waves they generate by the streaming instability (Kulsrud & Pearce 1969). Both interstellar and self-excited MHD turbulence should be severely damped by ion-neutral interactions in these neutral clouds, but CR gyration along tangled magnetic-field lines as well as the perpendicular (super)diffusive wandering of the lines also induce an effective CR diffusion in space, with a mean free path close to the magnetic-field coherence length (Lazarian 2016; Hu et al. 2022; Sampson et al. 2023). We therefore compare the dispersion in gas velocity and in magnetic bending in the two clouds. All CR transport modes strongly depend of the state of the MHD turbulence, so we derive the ß_sk_y magnetic-field strength from dust polarisation observations and we compare the sonic Mach numbers and gas-to-magnetic pressure ratios along the cores of both filaments.

5.1 Gas densities in the atomic phase

Given the elongated aspect of the Reticulum and Eridu filaments, we used two different assumptions for the cloud geometry, namely a cylindrical fibre or a shell seen edge on. These geometries bracket the unknown depth along the line-of-sight (LOS) in order to estimate HI gas volume densities, n_HI. In the fibre configuration, the LOS depth is close to the observed width in the plane of the sky. In order to estimate the latter, we used the Gaussian full widths at half max (FWHM) of the N_HI column-density profiles observed along cuts perpendicular to the filament axis. Figure 10 shows the eight and seven cuts respectively chosen along the Reticulum and Eridu filaments. The resulting depths are shown in the middle panel of Fig. 11 with an uncertainty reflecting that in cloud distance. The mean depth and rms dispersion are 5.0 ± 0.6 pc in Reticulum and 10.5 ± 1.0 pc in Eridu. In the shell shaped configuration, the LOS depth is closer to half the chord length through the thick shell. Given the mean widths measured above and a shell radius of curvature of about 100 pc, we estimated the LOS depths to be roughly 20 pc for Reticulum and 30 pc for Eridu. These limits are represented as triangles in the middle plot of Fig. 11. They provide reasonable upper limits to the cloud depth.

We derived the mean gas volume density in a lº-wide square centred on the peak of each N_HI profile cutting perpendicularly through the Reticulum and Eridu filaments. The results are shown in the bottom plot of Fig. 11. The shell geometry provides lower limits to the gas densities.

Figure 11 illustrates that the larger N_HI column densities recorded in the Eridu filament are likely due to longer LOS depths and that the Reticulum filament has denser substructures than Eridu. In both geometries, the density estimates along Reticulum are typical of the CNM and they are consistent with the detection of comparable column densities in the DNM. For both geometries, the Eridu volume densities are more typical of the unstable lukewarm neutral medium (LNM) that transitions between the stable WNM and CNM phases. They indicate that CNM clumps are likely present inside both clouds, but should be more numerous and more densely packed in Reticulum than in Eridu.

Fig. 10 NHI maps of the Reticulum (top) and Eridu (bottom) filaments, smoothed by a Gaussian kernel of 0.3° and 0.2°, respectively, to highlight the cloud edges. The white lines show the perpendicular cuts used to estimate the cloud widths. The cyan squares show lº-wide squares adopted to estimate gas volume densities, Bsky field strengths and parallel diffusion coefficients ĸ‖ along the filament core.

Fig. 11 Distributions of HI column densities (top), LOS depths (middle), and mean HI volume densities (bottom) in the Reticulum (magenta) and Eridu (green) clouds. The offset distance is the separation between a 1°-wide square along a filament and the first one (of largest latitude or declination). The different symbols refer to the fibre or shell cloud geometry. The latter provides upper/lower limits to the depth/density. The LOS depth shading is due to the distance uncertainty. The data points and error bars in volume density respectively give the mean value in each ~5 pc-wide square and the rms dispersion at 0.50.7 pc scale inside each square.

5.2 Self-confined diffusion

CRs gyrate and flow along magnetic field lines. They scatter off magnetic fluctuations with sizes close to their gyro-radius (e.g. ~2 µpc for a 10 GeV CR in a B = 5 µG field). For instance, in the frame moving with Alfvén waves along the field, the Lorentz force acting on a CR during a gyro-resonant interaction changes the particle’s pitch angle (tan θ = p Λ B/p • B), but not its energy, so random phase shifts between the gyro-motion and the encountered waves cause a random walk in pitch angle and reduce the CR effective speed along the field line. See Ruszkowski & Pfrommer (2023) for a nice pictural description of gyro-resonant scatterings in pitch angle. Whereas very high energy CRs are likely scattered by interstellar MHD turbulence, the lower-energy CRs that we sample here in γ rays are more likely scattered by Alfvén waves that they self-generate via the resonant streaming instability (Kulsrud & Pearce 1969). The anisotropy of interstellar MHD turbulence increases with decreasing scale in the cascade (Goldreich & Sridhar 1995). The eddies become elongated along the field because of turbulent reconnection. As a result, during the resonant encounter of a CR with a gyro-radius-long eddy, the gyro-orbit perpendicular to the field encompasses other uncorrelated eddies that cancel out the net change in pitch angle (Chandran 2000). Unlike Alfvén and slow modes, fast compressive turbulent modes remain much more isotropic, and are therefore efficient scatterers on small scales. In the neutral atomic gas, however, they mostly affect CRs of much larger energies (≫ 100 TeV) than those of interest here (Xu et al. 2016).

In the self-confinement scenario, a slight anisotropy in the CR pitch-angle distribution can excite Alfvén waves with the right geometry to induce further pitch-angle scattering (Kulsrud & Pearce 1969). The required level of anisotropy to trigger the instability corresponds to a bulk drift velocity v_D larger than the ion Alfvén speed (the drift velocity is that of the frame where the CR distribution is isotropic). The wave growth rate by this streaming instability depends on the integrated number density of exciting CRs and on the mass density of interstellar gas ions, ρ_i, as: (3)

where B is the mean field strength and we add the contributions of different CR nuclei with charge Z_je, rigidity 𝓡_j, and volume number density spectrum n_j(𝓡) with a power-law spectral index α_j in the range of interest (Kulsrud & Pearce 1969; Kulsrud & Cesarsky 1971). The CR density is integrated in rigidity above the threshold corresponding to the resonance condition 𝓡_j = B/k for the wave number k. We considered the spectra of CRs with Z_j ≤ 28 that reach the heliopause (Boschini et al. 2020). They correspond to the q_LIS emissivity discussed above. We integrated the spectra over 1 GV in rigidity and fitted the α_j spectral indices in the 1–20 GV range that contributes most to the growth rate.

The actual level of CR anisotropy that causes the streaming instability and the growth of scattering waves is extremely difficult to estimate from observations and from theory. We therefore follow Jiang & Oh (2018) and Bustard & Zweibel (2021) in expressing the net drift velocity with respect to the Alfvén wave frame as a diffusive flux: (4)

with the characteristic length for a change in CR flux, ℓ_CR, and the diffusion coefficient κ_‖ parallel to the mean B field.

Wave amplitudes and scattering rates are reduced by wave damping processes that depend on the properties of the background MHD turbulence and on ambient plasma conditions. The main damping mechanism in the neutral gas phases that we have probed in γ rays is the ion-neutral one. Following Hopkins et al. (2021) and Xu & Lazarian (2022), we checked that the linear Landau, non-linear Landau, and turbulent damping rates are two to three orders of magnitude smaller than the ion-neutral one in the typical environments of the Reticulum and Eridu clouds.

Momentum exchange between the ions and hydrogen or helium neutral atoms, all with the same equilibrium temperature T, leads to the damping rate given by Soler et al. (2016): (5)

where k_B is the Boltzmann constant, m_H is the H atom mass, x_HI and x_He are the neutral hydrogen and helium abundance by number, relative to the total number of hydrogen nucleons, and A_He=4 is the atomic weight of helium. The ion-neutral momentum-transfer cross sections are σ_iH = 10^–14 cm² and σ_iHe = 3 × 10⁻¹⁵ cm² from Vranjes & Krstic (2013).

The ion abundance, x_i, and mean ion atomic weight, A_i, vary with the gas state (Draine 2011). The predominant form of ions gradually varies from H⁺ in the WNM to about half H⁺ and half C⁺ in the CNM, and their total abundance decreases from about 2% in the WNM down by two orders of magnitude in the CNM. We used the heating, cooling, chemistry, and ionisation model of Kim et al. (2023) to calculate the equilibrium gas temperature T, ion abundance x_i, and average ion weight A_i as a function of the gas densities observed in the neutral phase (n_HI) of the two clouds, and for the interstellar conditions of the local ISM (x_He = 0.1, x_C = l·6 × 10⁻⁴ in the gas phase, solar metallicity, primary CR ionisation rate of 2 × 10⁻¹⁶ s⁻¹, and local FUV intensity). This allowed us to estimate the damping rates of Alfvén waves in different locations along both clouds.

Another formalism for ion-neutral damping is given by Kulsrud & Pearce (l969): (6)

where n_n is the neutral number density, m_n is the mean mass of neutrals, and m_i = A_im_H is the mean mass of ions. The ion-neutral momentum transfer rates are 〈σv〉_in = 3.25 × 10⁻⁹ cm³ s⁻¹ for collisions between H and H⁺, 2.39 × 10⁻⁹ cm³ s⁻¹ between H and C⁺, and 1.42 × 10⁻⁹ cm³ s⁻¹ between He and H⁺ (from Table 2.1 of Draine 2011). This formulation does not depend on the gas temperature.

We adopted the Soler et al. (2016) damping rate as the baseline for this work. Assuming a uniform LOS depth equal to the mean value of the estimates shown in Fig. 11 in each geometry, we find rather uniform damping rates close to 6 × 10⁻⁹ s⁻¹ across both clouds, in both geometries. The two damping rate formulations yield comparable values across the clouds in the shell geometry, but not in the denser fiber configuration where the Kulsrud & Pearce (1969) formulation yields a few times larger rates than the Soler et al. (2016) ones (one to three times larger along Eridu and one to five times larger along Reticulum). This difference is due to the drop in temperature and increase in ion atomic weight in the denser fibre gas in the Soler et al. (2016) formulation.

We further note that these rates likely over-predict the actual damping rate in the partially ionised plasma because they do not take the charged dust grains into account even though each grain carries a charge that varies from O.le (small grains) to 2Oe or 3Oe (large grains) in the atomic gas phases of interest here (Draine 2004). The presence of dust charges reduces the dissipation of Alfvén waves by at least an order of magnitude (Hennebelle & Lebreuilly 2023). Charged grains also lower the drift velocity required for CR excitation of the waves. The gas-to-dust mass ratio being of order 100, the ion Alfvén speed is limited to about ten times the gas one whereas, if dust is ignored, one expects much larger speeds in all gas densities above 1 cm⁻³ where the abundance of gaseous ions drastically falls below 1% (Kim et al. 2023). Charged dust therefore allows larger wave intensities for gyro-resonant scattering in the CNM and denser gas phases (Hennebelle & Lebreuilly 2023).

A quasi-steady state is reached when the growth rate Γ_grow of the gyro-resonant waves is balanced by the damping rate Γ_in, giving a parallel diffusion coefficient: (7)

We note that, under the steady-state self-confinement transport assumption, the diffusion coefficient depends on two opposite trends that limit its potential range inside a cloud since both the gas temperature and ionisation fraction rapidly drop as the gas density increases. We also note that the κ_‖/ℓ_CR formulation for the net drift velocity cancels the κ_‖ dependence on the ion Alfvén speed, and therefore also the dependence on the mean magnetic field strength.

The distance that separates the Reticulum and Eridu filaments ranges from 170 pc to 270 pc given their respective length and distance uncertainty to the Sun. They are also about 150 pc away from other clouds exhibiting γ-ray emissivities close to the average q_LIS. Given the 40–60% change in CR flux between Eridu and the other clouds, we have adopted ℓ_CR=400 pc as a typical CR gradient length scale for the CR energies probed with the LAT. This macroscopic length scale is motivated by differences in CR flux inside various clouds and not by the actual CR gradients that prevail in the immediate vicinity of the Eridu and Reticulum clouds, or those that are induced by magnetic bottlenecks, changes, or CR losses in the cloud envelopes (Schlickeiser et al. 2016; Bustard & Zweibel 2021). Given its uncertainty, we explicitly quote the ℓ_CR dependence in our results.

Figure 12 presents maps of the parallel diffusion coefficient obtained for both clouds when using a uniform LOS depth across the whole cloud (the mean one in the fibre configuration or the upper limit in the shell geometry). We further relaxed the uniformity of the LOS depth by calculating the diffusion coefficient in each of the 1° × 1° squares chosen along both filaments, using the fibre LOS depth estimated in each square. The results are shown in Tables B.1 and B.2 and as a function of the mean HI gas density in each square in Fig. 13. For a characteristic CR gradient scale length ℓ_CR of 400 pc, we find κ_‖ values that typically range from 6 × 10²⁸ to 14 × 10²⁸ cm² s⁻¹ across the Reticulum cloud in both geometries, and over a smaller range from 5 × 10²⁸ to 8 × 10²⁸ cm² s⁻¹ for Eridu because of its lower contrast in gas column density. The median κ_‖ value in all the pixels composing the squares along the core of each filament differs by less than 20% between the two clouds in each geometry.

The maps show that the diffusion coefficient does not vary substantially as the gas density increases from the cloud envelopes to the cores because the increase in gas density is tempered by the loss in ion abundance and because the damping rise is tempered by the drop in gas temperature and increase in ion mass. The overall variation from edge to core is of the order of two to three. We stress that these κ_‖ estimates are based on the cloud properties sampled at parsec scales rather than at the actual micro-parsec scales of gyro-resonance, so in-situ variations may well be much more severe.

Figure 12 illustrates that, when using the same CR flux to excite Alfvén waves in the Reticulum and Eridu clouds, both filaments exhibit very similar diffusion properties over most of their area in each of the fibre or shell configurations. The larger κ_‖ values found along the spine of Reticulum in Fig. 13 correspond to the densest cores where CRs can flow more rapidly due to the enhanced damping. In the self-confinement theory, CRs create the magnetic fluctuations they interact with, so fewer CRs have less power to generate perturbations to modify their pitch angle. The CR flux reduction in Eridu with respect to Reticulum leads to an increase in diffusion length by the flux ratio of 1.57 ± 0.09. Taking this ratio into account, Fig. 12 shows that CRs should diffuse at about the same speed in the densest parts of Eridu and Reticulum, but twice as fast at the periphery of Eridu than in the envelope of Reticulum, in both geometries. However, the diffusivity difference between the two clouds is too small and the uncertainties in the growth and damping rates are too large to relate the loss of CRs in Eridu to different self-confinement efficiencies.

Cosmic rays transfer energy to the magnetic fluctuations if they have an anisotropic distribution with a degree of anisotropy larger than , therefore if their drift velocity v_D along the mean magnetic field exceeds (Kulsrud & Pearce 1969). If the CRs are well scattered, there is a direct relationship between their anisotropy and spatial density gradient (Zweibel 2013), so a difference in CR flux inside Eridu and Reticulum may originate from different CR gradients in their immediate environment (different ℓ_CR values). Both clouds lie near the dense walls of the Local Valley, but Eridu stands in the vicinity of the Orion-Eridanus superbubble where steep CR gradients are possible even though we have found no evidence for fresh acceleration or reacceleration of CRs in the compressed outer rim of the superbubble (Joubaud et al. 2020). Steeper gradients, however, would excite more streaming instabilities and decrease the diffusion coefficient rather than increase it compared to Reticulum.

On the other hand, Reticulum may be a cooling fragment of an old shell (possibly from a supernova remnant) surrounded by steep CR gradients to efficiently confine the particles inside the shell, but a local boost of the wave growth rate inside this cloud is hard to reconcile with the fact that its CR flux matches that found in other local clouds with different histories and environments. The detection of Hα recombination emission towards Reticulum reveals the presence of ionised gas close to the cloud. Efficient streaming and gyro-scattering in the ionised gas may have reduced the CR anisotropy before entering the cloud, thereby preventing the further excitation of scattering waves inside the cloud. But this scenario would lead to faster diffusion in Reticulum than in Eridu since we find no ionised shell around the latter.

We note that κ_‖ linearly scales with ℓ_CR. This is the reason why the present estimates drastically differ from those found by Armillotta et al. (2021) for the same gas density. They find very smooth CR densities in the CNM with about two orders of magnitude larger gradient scale lengths in their simulations (Lucia Armillotta, priv. comm.).

The κ_‖ values we find in these individual clouds in the self-confinement scenario compare well with the range of recent estimates obtained for the entire Milky Way for isotropic diffusion. The latter vary from 4 × 10²⁸ to 20 × 10²⁸ cm² s⁻¹ for 10 GeV/n CRs (Jóhannesson et al. 2016, 2019; Weinrich et al. 2020; De La Torre Luque et al. 2022). They are based on direct measurements of CR nuclear spectra measured inside and near the heliosphere (by Voyager, ACE-CRIS, HEAO-3, CREAM, AMS-02, and PAMELA) and on different propagation codes to model spallation reactions (namely GALPROP, USINE, and DRAGON-2). The different estimates vary by a factor five even though the models merely try to constrain an isotropic, homogeneous, environment-independent coefficient in the Galaxy. The convergence between the Galactic averages and the values presented in Figs. 12 and 13 is nice, but it primarily relies on the choice of 400 pc for the ℓ_CR gradient scale and on the assumption of CR self-confinement in these clouds. As discussed above, variations in the cloud properties (temperature, ionisation fraction, and gas density) tend to compensate each other to reduce the range of κ_‖ values in this transport scenario. The choice of ℓ_CR was not tuned to match Milky Way values, but motivated by the 3D separation between clouds exhibiting different γ-ray emissivities. Nevertheless, while we can compare the diffusion coefficients presented in Figs. 12 and 13 for the two clouds because they rely on the same assumptions and on the same spatial scale and methods to measure the relevant interstellar parameters, we should consider the absolute value of these coefficients with great care because of the unknown state of the plasma and unknown level of CR anisotropy at the gyro-resonance scale.

Fig. 12 Parallel diffusion coefficients obtained for the Reticulum (top) and Eridu (bottom) clouds in the fibre (left) and shell (right) geometries, assuming the damping rate of Eq. (5), a characteristic CR flux variation length ℓCR = 400 pc, and a uniform LOS depth equal to the shell one or to the average of the fibre ones. For each geometry, the Eridu plots show ĸ‖ values for the same CR flux exciting the Alſvén waves as in Reticulum (left) or for the reduced CR flux corresponding to the cloud lower γ-ray emissivity (right).

Fig. 13 Estimates of the parallel diffusion coefficient in the 1°-wide squares chosen along the core of the Reticulum (magenta) and Eridu (green) filaments, for the fibre (top) and shell (bottom) geometries, for the same CR flux pervading the two clouds (filled symbols) and for the reduced CR flux in Eridu (open symbols). The data points and error bars respectively give the mean value in each ~5 pc-wide square and the rms dispersion at 0.5–0.7 pc scale inside each square. Different nHI and κ‖ scales are used for each geometry.

5.3 Velocity dispersion and sonic Mach numbers

As indicated when opening the discussion section, other CR transport modes relate to the turbulent wandering and bending of magnetic-field lines, in particular in super-Alfvénic turbulence (Lazarian & Xu 2021; Hu et al. 2022; Xu & Lazarian 2022; Sampson et al. 2023, and references therein). We therefore explored the turbulence level in gas velocity and in B_sky orientation in each cloud, and we estimated the sonic Mach numbers of the turbulence to gauge the plasma compressibility.

We took advantage of the pseudo-Voigt decomposition of the HI spectra we used for the cloud separation to follow changes in LOS velocity for the main HI line across each cloud. Figure 14 presents the maps of the σ_υ line widths of the main line across Eridu (top panel) and Reticulum (bottom panel). We also mapped the position-to-position velocity dispersion at the smallest scale resolved by the GASS using a two-point structure function in LOS velocity: (8)

where the velocity dispersion compares the line central velocity in a given direction x and in a displaced direction v(x + δ_i). The average is taken over displaced directions spanning an annulus δ/2 ≤ δ_i ≤ 3δ/2 around the x direction, with δ equal to the effective GASS FWHM angular resolution of 16′.2 (Kalberla et al. 2010). The resulting S_v maps are shown for the Reticulum and Eridu clouds in the top and bottom panels of Fig. 14, respectively.

The δ_υ line widths are commensurate with the thermal broadening expected from the temperature map inferred in the previous section. The velocity resolution of GASS and the uncertainty in the modelled gas temperature do not allow a reliable subtraction of the thermal broadening to measure the dispersion in turbulent velocities. The σ_υ line widths therefore provide upper limits to the turbulent motions along the lines of sight. On the other hand, because of the average over the δ ring, the S_v structure function is closer to a lower limit to the turbulent motions along the line of sight. Both S_v and σ_υ are one-dimensional measures of the turbulent velocity field.

Figure 14 shows that the line width and velocity dispersion values decrease from the cloud edges to the cores of both clouds, as expected from CNM regions. It also indicates that the dispersion in the velocity field of the Eridu cloud is about twice larger than in the more compact Reticulum filament (both in S_v and σ_v). This is expected for a more diffuse cirrus (Audit & Hennebelle 2010).

We calculated the sonic Mach numbers of the turbulence in the 1°-wide squares sampling the core of each filament in the case of isotropic turbulence . Figure 15 shows that these estimates are consistent with the distribution of Mach numbers obtained in cold CNM clumps by Heiles & Troland (2003) when comparing the HI line width with the thermal broadening inferred from spin temperatures. Figure 15 indicates that the turbulence is mildly supersonic along both Eridu and Reticulum in the fibre configuration and it is trans-sonic in the warmer shell geometry. The Mach numbers compare well in both clouds and the turbulence appears to be moderately compressible in both of them.

Fig. 14 Maps of the velocity line width, συ, and of the velocity dispersion function, Sv(δGASS), of the main HI emission line across the Eridu (top) and Reticulum (bottom) clouds. We note the different colour scale for each cloud.

Fig. 15 Sonic Mach numbers of turbulent motions in the 1°-wide squares chosen along the Eridu (green) and Reticulum (magenta) filaments, in the fibre (top) and shell (bottom) geometries, assuming isotropic turbulence and using the Sv dispersion function (dark filled symbols) or the σv line width (light open symbols) to trace the one-dimensional turbulent velocity dispersion. The data points and error bars respectively give the mean value in each ~5 pc-wide square and the rms dispersion at 0.5–0.7 pc scale inside each square. The histogram shows the Mach number distribution in cold CNM clumps from Heiles & Troland (2003).

5.4 Dispersion in the B_sky magnetic field orientation

In order to study the geometrical complexity of the magnetic field lines in the Reticulum and Eridu clouds we used the dust polarisation data recorded by Planck at 353 GHz. We degraded the original 4′.8 resolution to an angular resolution of 14^'(FWHM) in order to increase the signal to noise ratio. The data come from the third public data release (PR3), which is described in the Planck 2018 Release Explanatory Supplement⁴. The Q and U Stokes parameters provided in the Planck database are defined with respect to the Galactic coordinates, but in a direction opposite to the IAU convention along the longitude axis, so we define the polarisation angle: (9)

where ψ is zero at the north Galactic pole and it follows the IAU convention by increasing counter-clockwise with increasing longitude. The arctan(sin, cos) is the arc-tangent function that solves the π ambiguity taking into account the sign of the cosine.

The nearby Reticulum and Eridu clouds are the dominant dust clouds along their respective lines of sight, so the dust polarisation measurements give us access to the magnetic field in the plane of the sky, B_sky, in those clouds. In order to quantify the regularity of the B_sky direction, we used the polarisation angle dispersion function (Planck Collaboration Int. XIX 2015): (10)

where Δψ_xi = ψ(x) − ψ(x + δ_i) is the difference between the polarisation angle in a given direction x and in a displaced direction x + δ_i. In practice, we calculated Δψ_xi from the Stokes parameters as (Planck Collaboration Int. XIX 2015): (11)

where the x and i indices stand for the central and displaced directions, respectively. The average in S_ψ is taken over displaced directions spanning an annulus around the x direction with an inner radius δ/2 and an outer radius 3δ/2. We took 15′ for the δ lag.

The S_ψ dispersion values range from 0° for a well ordered B_sky magnetic field to S_ψ = 90° for a disordered one. When the polarisation signal is dominated by noise, S_ψ converges to . The results on the B_sky field orientation and on its dispersion are displayed as field lines in Fig. 1, as ψ angle maps in Fig. 16, and as Sψ angle dispersion maps in Fig. 17.

Atomic filamentary clouds have mean magnetic fields often aligned with the filament axis (Clark et al. 2014; Clark & Hensley 2019) whereas the field orientation turns to perpendicular or to disorder in the dense molecular parts (Soler et al. 2013; Planck Collaboration Int. XXXV 2016). The atomic Eridu and Retic-ulum clouds follow the expected trend. Figures 1, 16, and 17 jointly show that the B_sky magnetic field is well ordered along the axis of both filaments, with a stable direction along nearly their whole length, and with a dispersion lower than 10° on a 1-pc scale over large parts of the filament extent. Moreover, these ordered fields are almost perpendicular to the orientation of the Galactic magnetic field inferred in the solar neighbourhood from dust polarisation in local clouds (Planck Collaboration Int. XLIV 2016) or by the field draping around the heliosphere (Zirnstein et al. 2016).

The Reticulum field appears to be slightly more ordered in its dense part than in Eridu, likely because of the compression that originally formed the compact fibre or shell. The dispersion in field orientation appears to be rather uniform across the whole of Eridu and is a little more contrasted between the filament axis and edges of Reticulum. The parsec-scale data, however, do not suggest that CRs have slowed down in Reticulum because of a much more tangled field. To the contrary, compared to Eridu, they should flow along more ordered fields over a large fraction of the filament length.

Because of the lateral super-diffusion of turbulent magnetic field lines, CRs with large enough pitch angles can bounce back on magnetic mirrors and leap onto adjacent lines (see Fig. 1 of Lazarian & Xu 2021) whereas the low pitch-angle CRs (in the loss cone for mirroring) keep flowing on and undergo gyro-resonant pitch-angle scattering. Mirror diffusion efficiently slows down the effective CR speed along the mean field, but only for the small fraction of particles with large pitch angles, so the overall reduction in κ_‖ is small. It is typically less than a factor of 2 and fades out for CR energies below 100 GeV (Lazarian & Xu 2021). The fact that we find comparable dispersion levels in B_sky orientation in the two clouds does not suggest widely different levels of mirror diffusion, unless the field gets much more tangled in the more numerous CNM clumps of Reticulum at unresolved scales below one parsec. We are investigating this possibility with interstellar MHD simulations.

Fig. 16 Map of the polarisation angle ψ at 353 GHz with 14′ resolution in the Eridu (top) and Reticulum (bottom) clouds. The angle gives the orientation of the Bsky field. The black contours delineate the filaments at NHI column densities of 3 × 1020 cm−2 and 1.2 × 1020 cm−2, respectively.

Fig. 17 Map in the Eridu (top) and Reticulum (bottom) clouds of the polarisation angle dispersion function Sψ at 353 GHz with 14′ resolution and δ=15′ lag. The black contours delineate the filaments at NHI column densities of 3 × 1020 cm−2 and 1.2 × 1020 cm−2, respectively.

Fig. 18 Plane-of-sky magnetic-field strengths in the 1°-wide squares chosen along the Reticulum (magenta) and Eridu (green) filaments, in the fibre (top) and shell (bottom) geometries, using the Sv dispersion function in the ST (dark filled symbols) or DCF method (light open symbols). The data points and error bars respectively give the mean value in each ~5 pc-wide square and the rms dispersion at 0.5–0.7 pc scale inside each square. The offset distance of the DCF data points has been slightly enlarged for clarity. The grey band shows the mean CNM Bsky value from Heiles & Troland (2005) assuming .

5.5 Magnetic-field strength

The total magnetic field is composed of a mean component, B₀, and a small fluctuating one, |δB| ≪ |B₀|, so that the magnetic energy density is . To infer the mean field strength from dust polarisation, the Davis-Chandrasekhar-Fermi (DCF) method assumes that isotropic turbulent motions initiate incompressible, transverse, Alfvén waves (δB • B₀ = 0) in the frozen field (Davis 1951; Chandrasekhar & Fermi 1953). They also assume that the kinetic energy of the turbulence is equal to the fluctuating magnetic energy density, and that the dispersion in polarisation angle ψ traces the dispersion in magnetic-field orientation. From the equipartition condition, and from S_ψ ≈ δB/B₀ one obtains: (12)

where the f correction factor varies with authors between 0.3 and 0.7, for instance for isotropic turbulence (Chandrasekhar & Fermi 1953). Other values are reviewed by Skalidis & Tassis (2021, Sect. 2.2.3). We adopted their f = 0.5 recommendation for this work.

The DCF method tends to overestimate the actual field strength because of line-of-sight and beam averaging in the data and because the method overlooks the strong anisotropy in ISM turbulence and the role of compressible modes in the cross-term δB.B₀ ≠ 0 (Skalidis & Tassis 2021). To include the compressible modes, these authors propose an alternative formulation that agrees within 20–30% with MHD simulation data: (13)

We applied both methods to calculate the plane-of-sky magnetic-field strengths in the 1º-wide squares chosen along the Reticulum and Eridu filaments, with the velocity dispersion function S_v as a proxy for the turbulent velocity δv. The results are shown in Fig. 18 and in Tables B.1 and B.2. Following the discussion by Skalidis & Tassis (2021), the and estimates likely bracket the true B_sky value, the latter being about 40% (50%) larger than the former in the case of Reticulum (Eridu). We also calculated upper limits to B_sky using the σ_v line width as an upper limit to the turbulent velocity. The results are provided in Tables B.1 and B.2. They are about three times larger than the S_v -based estimates for each method.

Figure 18 shows that the B_sky field strength is rather uniform along the core of each cloud and that it is about one and a half time larger along Eridu than along Reticulum: around 3.1 (1.7) µG for Eridu and 1.9 (0.9) µG for Reticulum in the fibre (shell) configuration if we use the S_v velocity dispersion to trace the turbulent motions. Three times larger upper limits to the fields are found when using the σ_v line width instead. The difference between Eridu and Reticulum does not depend much on the cloud geometry, but one should keep in mind that the magnetic-field component along the line of sight may significantly differ in both clouds. Figure 18 also illustrates that the B_sky strengths in Eridu and Reticulum compare reasonably with the median total field B_tot = 6 ± 1.8 µG inferred from Zeeman splitting of HI lines and Monte Carlo simulations (Heiles & Troland 2005) if one assumes for a random field orientation with respect to the line of sight (Heiles & Troland 2005).

We cannot estimate whether the MHD turbulence is super-Alfvénic or not since both B_sky derivations assume equipartition between the fluctuating kinetic and magnetic energy densities, whether the latter is dominated by δB²/8π in the case of incompressible turbulence (DCF) or by δBB/4π when compressible modes are present (ST). In both cases, the Alfvén Mach numbers will be close to one. The large-scale order we observe in the Bsky orientation in both clouds suggests sub- or trans-Alfvénic turbulence as often found for the atomic gas phases in interstellar MHD simulations.

We can estimate the gas magnetisation level from the squared ratio of the sound to Alfvén speed, (Soler et al. 2013), which is close to the gas to magnetic pressure ratio. To derive the Alfvén velocity of the total gas, we assumed and the ST formulation for B_sky. Figure 19 shows the normalised distributions of the ratios found in the 1 ^º-wide squares sampling the gas along the core of both filaments. If we use the B_sky(S_v) estimate, the gas and magnetic pressures compare well across both clouds in their fibre geometry and the gas pressure slightly exceeds the magnetic one in the shell configuration. The magnetisation is skewed to values below one (i.e. to dominant magnetic pressure) when using the B_sky(σ_v) upper limits to the field, but not severely, so equivalent gas and magnetic pressures may prevail over large parts of the clouds. Therefore, whether in terms of B_sky orientation dispersion or in terms of magnetic stiffness against turbulent gas motions, we find no important difference between the two clouds at a parsec scale that would suggest more magnetic-field tangling and slower CR diffusion in one of them.

Fig. 19 Histograms of the gas magnetisation levels in the 1º-wide squares chosen along the Eridu (green) and Reticulum (magenta) filaments, for each geometry, and using the Sv velocity dispersion (filled bars) or the σv line width (markers) in the ST formulation to infer Bsky and the total magnetic-field strength.

6 Conclusions

We compared the CR flux in the atomic phase of two nearby clouds that have many properties in common. They are both elongated filaments, which cross the same hollow sector of the Local interstellar Valley at comparable distances from the Sun. Both are threaded by ordered magnetic fields that are mostly aligned with the long axis of the filament and that are largely inclined with respect to the Galactic plane and to the large-scale magnetic-field direction in the local ISM. Both clouds have magnetic pressures of the same magnitude as the thermal ones. Nevertheless, we find that the γ-ray emissivity per gas nucleon in one cloud, Reticulum, is 1.57 ± 0.09 times larger than in the other cloud, Eridu, over the whole energy range of analysis, from 0.16 to 63 GeV. The difference cannot be attributed to uncertainties in gas mass because the correlation between the dust and γ-ray data allows an efficient detection of any additional gas that is not traced by HI line emission (as we find in Reticulum) and the γ-ray contribution of this additional gas is included in the model. Furthermore, we cannot reduce the CR flux ratio below 1.32 ± 0.08 when allowing for strong differences in HI spin temperature between the two clouds (at the expense of a poorer γ-ray fit).

Beyond characterising the global magnetic orientation and average field strength of a few microGauss in the clouds, we studied turbulence characteristics at a parsec scale that might hint at different environmental conditions for CR diffusion, namely the dispersion in gas velocity and in magnetic-field orientation, differences in the sonic Mach numbers (compressibility) of the turbulence, and differences in magnetisation level (magnetic stiffness against gas turbulent motions). We find that the two clouds compare well in all these aspects. The only notable difference is that the core of the Reticulum cloud has collapsed (or has been compressed) to larger densities typical of the CNM and DNM phases, whereas a large fraction of the more diffuse gas in Eridu may be transitioning through the unstable, lukewarm phase. MHD simulations will help in studying the level of magnetic-field tangling in the different atomic and DNM gas phases, or in shell compression behind a slow shock, in order to search for changes that may impact the collective drift of the CR population, such as changes in the turbulent energy fraction of fast compressive modes or in the density of magnetic mirrors (Lazarian et al. 2023).

Alternatively, the 1 –500 GeV/n CRs we probed in γ rays may primarily scatter off MHD waves they have excited through streaming instability. If so, one expects a stronger damping in the denser Reticulum cloud, but we show that the difference in expected steady-state diffusion coefficient between the two clouds can be small (less than a factor of two) because the increase in diffusion length with the ambient gas density is moderated by the corresponding drop in gas temperature and ionisation fraction. We plan to further investigate the impact of charged dust grains, which should reduce the damping rates along the core of the two clouds (Hennebelle & Lebreuilly 2023). However, we only probed the cloud properties at a par-sec scale, whereas the growth and damping of the relevant MHD waves occur at wavelengths of a few micro-parsecs where the interstellar gas properties cannot be explored. Moreover, in the self-confinement mode, the wave growth rate in each cloud relates to the degree of anisotropy in CR pitch angle or, equiv-alently, to the steepness of the local CR spatial density gradient (Zweibel 2013), which could differ in the vicinity of the clouds or in bottlenecks along their edge (different ℓ_CR).

The damping rates we estimate suggest that the magnetic perturbations propagating along the field lines should dissipate over typical lengths of tens to hundreds of micro-parsecs (, Spangler et al. 2011). The median of their distribution is twice as large in Eridu as in Reticulum. The dissipation lengths correspond to proton gyro-radii for energies in the ~10 GeV to ~ 1 TeV band, with a median energy three times larger in Eridu than in Reticulum. The waves must be constantly replenished or we would observe a spectral change in γ-ray emissivity in the corresponding ~1 to ~100 GeV band as the higher-energy particles would leak out more rapidly. The fact that we find a change in CR flux but not in spectrum between the two clouds is therefore an important constraint. The same overall spectrum has been observed in several clouds in the local ISM (Grenier et al. 2015), but with improved photon statistics the Fermi-LAT data can now serve to precisely test the uniformity of this spectrum in different types of atomic clouds and independently of the CR flux pervading them.

Acknowledgements

We thank Jeong-Gyu Kim for sharing his interstellar model for the solar neighbourhood conditions, Lucia Armillotta for sharing information from her work, the members of the Interstellar Institute (I2) CNRS International Research Network for enlightening discussions, and the referee for helping us clarify several points in the discussion. We acknowledge the financial support by the LabEx UnivEarthS (ANR-10-LABX-0023 and ANR-18-IDEX-0001) for this work. The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K.A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.

Appendix A Distance estimates to the nearby clouds

We used the 3D dust distribution reconstructed from stellar surveys by Leike et al. (2020) out to distances of about 400 pc from the Sun to place constraints on the distance to different cloud complexes seen in the analysis region. Figure A.1 displays the dust extinction in the G band in contiguous 20-pc-wide intervals. The lack of dust information at latitudes below about 50° for the three most distant intervals is due to the limited size of the 3D data cube perpendicular to the Galactic plane (+ 270 pc). In addition to those maps, one can plot the extinction profile as a function of distance along sight lines of large N_H column density in each cloud. Examples are shown in Fig. A.2 when probing directions along the Reticulum filament. We used the peaks found in the dust profiles and the correspondence in shape found between atomic clouds in Fig. 3 and dust clouds in Fig. A.1 to estimate the cloud distances that are given in Sect. A.

Fig. A.1 Maps of the G-band extinction optical depth integrated over 20 pc-wide intervals in distance from the Sun across the analysis region. The data come from the 3D distribution of Leike et al. (2020). The distance ranges are labelled in each map. The black contours delineate the Reticulum HI cloud as in Fig. 1.

Fig. A.2 Dust extinction density profiles as a function of radial distance from the Sun for selected lines of sight along the Reticulum filament. The data come from the 3D distribution of Leike et al. (2020).

Fig. A.3 Maps of the G-band extinction optical depth integrated over 20 pc-wide intervals in distance from the Sun, towards the Orion-Eridanus region. The data come from the 3D distribution of Leike et al. (2020). The dust cloud that spans the region nearly horizontally around b = 50º and at distances between 160 and 260 pc marks the rim of the Orion-Eridanus superbubble. At larger distances, the black contours delineate the Eridu HI cloud as in Fig. 1. The distance ranges are labelled in each map.

We repeated the analysis to revisit the distance to the Eridu cloud. The latter, being rather diffuse, is more difficult to identify in the dust data, but Fig. A.3 suggests that it lies at about 300 pc from the Sun, beyond the edge of the Orion-Eridanus super-bubble that is clearly detected at closer distances around 200 pc as the long arc crossing the maps horizontally at latitudes around −50º.

Appendix B ISM properties in the 1°-wide squares along the cloud filaments

Tables B.1 and B.2 list the mean value and rms dispersion of the interstellar properties found in the 1^º × 1^º squares that are depicted in Fig. 10 and that span the length of the Reticulum and Eridu filaments, along their core. The properties include the LOS depth, the N_HI column density and n_HI volume density in the atomic gas, the HI velocity dispersion in the plane of the sky (S_v) and along the line of sight (σ_v line width), the gas temperature T and the abundance (x_ion) and atomic weight (A_ion) of ions inferred from the gas density and the Kim et al. (2023) model, the B_sky magnetic-field strengths estimated according to Eqs. 12 and 13 with the S_v velocity dispersion or σ_v line width as a proxy for the turbulent velocity δv, the sonic Mach numbers M_S of the turbulence for the same options in δv, the Alfvénic Mach numbers M_A of the MHD turbulence inferred with the ST and DCF methods, the ion-neutral damping rate according to Eq. 5, and the steady-state κ_‖ diffusion coefficient for self-confinement with the CR flux reaching the heliosphere in the 1–20 GV rigidity range.

References

All Tables

All Figures

Fig. 1 Dust optical depth τ353 at 353 GHz in Galactic coordinates towards the Eridu (top) and Reticulum (bottom) filaments. The black lines give the orientation of the plane-of-the-sky Bsky magnetic field inferred from the dust polarisation observations at 353 GHz and at an angular resolution of 14′ (FWHM). The red contours delineate the filaments at HI gas column densities of 3 × 1020 cm−2 and 1.2 × 1020 cm−2, respectively.

In the text

Fig. 2 Dust optical depth measured at 353 GHz from the Planck and IRAS data and displayed on a 10′.2×7′.5 grid at the angular resolution of Planck (top panel) and of Fermi-LAT (middle). The bottom panel shows the γ-ray photon counts produced by CR interactions with gas and recorded by Fermi-LAT in 0°25×0°25 bins in the 0.16–63 GeV energy band (other emissions of non-gaseous origin are subtracted). Three zones are excluded from the analysis region: the SMC on the left, LMC in the middle, and the edge of the Gum nebula in the upper right. The Reticulum cloud is the long filament crossing the region below the LMC.

In the text

Fig. 3 NHI column densities in the Reticulum (Ret), Far Hydras (HyiF), North Hydras (HyiN), South Hydras (HyiS), East Dorado (DorE), West Dorado (DorW), Horologium (Hor), and the Magellanic clouds for optically thin HI. The faint column densities in the negative-velocity map (NegVel) are multiplied by 4 for display. The maps include a 5º-wide peripheral band around the analysis region (black perimeter) to account for the wings of the Fermi-LAT PSF in the γ-ray model.

In the text

	Fig. 4 Photon yields on a 0º25 pixel grid from the main components of the γ-ray model in the 0.16–63 GeV band. The yields come from the HI and DNM column densities derived for the optically thin case in the different clouds, and from the isotropic background, Galactic IC emission, and γ-ray point sources.
In the text

	Fig. 5 Evolution of the log-likelihood ratio of the γ-ray fit in the 0.16– 63 GeV band as a function of the HI spin temperature used to calculate the NHI column densities in the clouds.
In the text

Fig. 6 Residual maps (data minus best-fit model, in sigma units) in γ-ray counts in the 0.16–63 GeV band (σγ) and in dust optical depth , obtained when including different sets of gaseous components in the models: without (top) and with (middle) the DNM components in addition to the HI templates. The column-density maps (bottom) extracted from the γ-ray residuals (left) and from dust residuals (right) are shown in the bottom row. The models are derived for an HI at the optically thin case.

In the text

Fig. 7 Distributions of the best-fit coefficients obtained in the thousand jackknife fits to the 0.16-63 GeV γ-ray data (two left columns) and dust data (two right columns) for the optically thin HI case. The vertical lines and overlaid numbers give the mean values of the distributions. The units are 1026 cm–2 for qDNM, 10–27 cm2 for yHI and yDNM, and 10–7 for yiso. qHI,i, qIC and qiso are dimensionless normalisation factors.

In the text

Fig. 8 Spectral evolution relative to the local interstellar spectrum qLIS of the γ-ray emissivities found for the HI and DNM phases of the Reticulum and North Hydrus clouds (top and middle panels). Spectral evolution of the isotropic and IC intensities found in the region relative to the input spectra (bottom panel), qDNM/qLIS is given in 1026 cm−2. The results are shown for the optically thin HI case. The error bars combine the statistical and systematic uncertainties.

In the text

Fig. 9 Spectral evolution of the γ-ray emissivities per gas nucleon in the Eridu (green), Reticulum (red), and North Hydrus (cyan) clouds. The orange curve gives the average qLIS emissivity measured in the atomic gas at latitudes 10° ≤ |b| ≤ 70° (Casandjian 2015). The emissivities are shown for the HI optical-depth correction that best fits the γ-ray data (Ts = 100 K for Eridu, 140 K for qLIS, and optically thin for the others). The shaded bands include the statistical and systematic uncertainties.

In the text

Fig. 10 NHI maps of the Reticulum (top) and Eridu (bottom) filaments, smoothed by a Gaussian kernel of 0.3° and 0.2°, respectively, to highlight the cloud edges. The white lines show the perpendicular cuts used to estimate the cloud widths. The cyan squares show lº-wide squares adopted to estimate gas volume densities, Bsky field strengths and parallel diffusion coefficients ĸ‖ along the filament core.

In the text

Fig. 11 Distributions of HI column densities (top), LOS depths (middle), and mean HI volume densities (bottom) in the Reticulum (magenta) and Eridu (green) clouds. The offset distance is the separation between a 1°-wide square along a filament and the first one (of largest latitude or declination). The different symbols refer to the fibre or shell cloud geometry. The latter provides upper/lower limits to the depth/density. The LOS depth shading is due to the distance uncertainty. The data points and error bars in volume density respectively give the mean value in each ~5 pc-wide square and the rms dispersion at 0.50.7 pc scale inside each square.

In the text

Fig. 12 Parallel diffusion coefficients obtained for the Reticulum (top) and Eridu (bottom) clouds in the fibre (left) and shell (right) geometries, assuming the damping rate of Eq. (5), a characteristic CR flux variation length ℓCR = 400 pc, and a uniform LOS depth equal to the shell one or to the average of the fibre ones. For each geometry, the Eridu plots show ĸ‖ values for the same CR flux exciting the Alſvén waves as in Reticulum (left) or for the reduced CR flux corresponding to the cloud lower γ-ray emissivity (right).

In the text

Fig. 13 Estimates of the parallel diffusion coefficient in the 1°-wide squares chosen along the core of the Reticulum (magenta) and Eridu (green) filaments, for the fibre (top) and shell (bottom) geometries, for the same CR flux pervading the two clouds (filled symbols) and for the reduced CR flux in Eridu (open symbols). The data points and error bars respectively give the mean value in each ~5 pc-wide square and the rms dispersion at 0.5–0.7 pc scale inside each square. Different nHI and κ‖ scales are used for each geometry.

In the text

	Fig. 14 Maps of the velocity line width, συ, and of the velocity dispersion function, Sv(δGASS), of the main HI emission line across the Eridu (top) and Reticulum (bottom) clouds. We note the different colour scale for each cloud.
In the text

Fig. 15 Sonic Mach numbers of turbulent motions in the 1°-wide squares chosen along the Eridu (green) and Reticulum (magenta) filaments, in the fibre (top) and shell (bottom) geometries, assuming isotropic turbulence and using the Sv dispersion function (dark filled symbols) or the σv line width (light open symbols) to trace the one-dimensional turbulent velocity dispersion. The data points and error bars respectively give the mean value in each ~5 pc-wide square and the rms dispersion at 0.5–0.7 pc scale inside each square. The histogram shows the Mach number distribution in cold CNM clumps from Heiles & Troland (2003).

In the text

	Fig. 16 Map of the polarisation angle ψ at 353 GHz with 14′ resolution in the Eridu (top) and Reticulum (bottom) clouds. The angle gives the orientation of the Bsky field. The black contours delineate the filaments at NHI column densities of 3 × 1020 cm−2 and 1.2 × 1020 cm−2, respectively.
In the text

	Fig. 17 Map in the Eridu (top) and Reticulum (bottom) clouds of the polarisation angle dispersion function Sψ at 353 GHz with 14′ resolution and δ=15′ lag. The black contours delineate the filaments at NHI column densities of 3 × 1020 cm−2 and 1.2 × 1020 cm−2, respectively.
In the text

Fig. 18 Plane-of-sky magnetic-field strengths in the 1°-wide squares chosen along the Reticulum (magenta) and Eridu (green) filaments, in the fibre (top) and shell (bottom) geometries, using the Sv dispersion function in the ST (dark filled symbols) or DCF method (light open symbols). The data points and error bars respectively give the mean value in each ~5 pc-wide square and the rms dispersion at 0.5–0.7 pc scale inside each square. The offset distance of the DCF data points has been slightly enlarged for clarity. The grey band shows the mean CNM Bsky value from Heiles & Troland (2005) assuming .

In the text

	Fig. 19 Histograms of the gas magnetisation levels in the 1º-wide squares chosen along the Eridu (green) and Reticulum (magenta) filaments, for each geometry, and using the Sv velocity dispersion (filled bars) or the σv line width (markers) in the ST formulation to infer Bsky and the total magnetic-field strength.
In the text

	Fig. A.1 Maps of the G-band extinction optical depth integrated over 20 pc-wide intervals in distance from the Sun across the analysis region. The data come from the 3D distribution of Leike et al. (2020). The distance ranges are labelled in each map. The black contours delineate the Reticulum HI cloud as in Fig. 1.
In the text

	Fig. A.2 Dust extinction density profiles as a function of radial distance from the Sun for selected lines of sight along the Reticulum filament. The data come from the 3D distribution of Leike et al. (2020).
In the text

Fig. A.3 Maps of the G-band extinction optical depth integrated over 20 pc-wide intervals in distance from the Sun, towards the Orion-Eridanus region. The data come from the 3D distribution of Leike et al. (2020). The dust cloud that spans the region nearly horizontally around b = 50º and at distances between 160 and 260 pc marks the rim of the Orion-Eridanus superbubble. At larger distances, the black contours delineate the Eridu HI cloud as in Fig. 1. The distance ranges are labelled in each map.

In the text