Gas shells and magnetic fields in the Orion-Eridanus superbubble

Aims. The Orion-Eridanus superbubble has been blown by supernovae and supersonic winds of the massive stars in the Orion OB associations. It is the nearest site at which stellar feedback on the interstellar medium that surrounds young massive clusters can be studied. The formation history and current structure of the superbubble are still poorly understood, however. It has been pointed out that the picture of a single expanding object should be replaced by a combination of nested shells that are superimposed along the line of sight. We have investigated the composite structure of the Eridanus side of the superbubble in the light of a new decomposition of the atomic and molecular gas. Methods. We used H I 21 cm and CO ( J =1 − 0) emission lines to separate coherent gas shells in space and velocity, and we studied their relation to the warm ionised gas probed in H α emission, the hot plasma emitting X-rays, and the magnetic ﬁelds traced by dust polarised emission. We also constrained the relative distances to the clouds using dust reddening maps and X-ray absorption. We applied the Davis–Chandrasekhar–Fermi method to the dust polarisation data to estimate the plane-of-sky components of the magnetic ﬁeld in several clouds and along the outer rim of the superbubble. Results. Our gas decomposition has revealed several shells inside the superbubble that span distances from about 150–250 pc. One of these shells forms a nearly complete ring ﬁlled with hot plasma. Other shells likely correspond to the layers of swept-up gas that is compressed behind the expanding outer shock wave. We used the gas and magnetic ﬁeld data downstream of the shock to derive the shock expansion velocity, which is close to ∼ 20 kms − 1 . Taking the X-ray absorption by the gas into account, we ﬁnd that the hot plasma inside the superbubble is over-pressured compared to plasma in the Local Bubble. The plasma comprises a mix of hotter and cooler gas along the lines of sight, with temperatures of (3–9) and (0 . 3 − 1 . 2) × 10 6 K, respectively. The magnetic ﬁeld along the western and southern rims and in the approaching wall of the superbubble appears to be shaped and compressed by the ongoing expansion. We ﬁnd plane-of-sky magnetic ﬁeld strengths from 3 to 15 µ G along the rim.


Introduction
The Orion-Eridanus superbubble is the nearest site of active high-mass star formation. It has been studied in particular to investigate the feedback of these stars on the interstellar medium. Together with supernovae and stellar winds, the intense UV stellar radiation has carved a 200 pc wide cavity that is filled with low-density ionised gas and expands into the interstellar medium. The material that was initially present was swept up and compressed to form a shell of neutral gas, the near wall of which may be as close as 180 pc from the Sun (Bally 2008). The boundary between the ionised gas and the neutral shell can be seen in H α emission that traces the recombination of ionised hydrogen. The H α emission that is detected in the region is displayed in Fig. 1. This figure features the H II region λ Orionis, the bright crescent of Barnard's Loop close to the Orion A and B molecular clouds, and the arcs of the Eridanus Loop to the west. The westernmost arc nicely delineates the outer rim of the superbubble, but the origin of the brightest vertical arc, called Arc A by Johnson (1978), is still debated. It may come from a shell inside the superbubble or from the superposition of separate structures along the line of sight. Pon et al. (2014) summarised arguments in favour and against the association of Arc A with the superbubble based on several tracers and methods, but did not firmly conclude. They mentioned that the molecular cloud identified as MBM 18 (for Magnani, Blitz and Mundi, near α J2000 = 60. • 6, δ J2000 = 1. • 3) by Magnani et al. (1985) coincides in direction and tentatively in velocity with Arc A. The detection of dust reddening fronts in the PanSTARRS-1 stellar photometric data place MBM 18 within 200 pc from the Sun (Zucker et al. 2019). No significant reddening is found beyond 500 pc in this direction; this strengthens the association between Arc A and the superbubble.
The H α emission has also been used to determine the outer boundaries of the superbubble. Whereas it appears rather clearly defined on the western side, it had been thought that the eastern boundary was Barnard's Loop until Ochsendorf et al. (2015) pushed the limit farther out to a fainter feature beyond the extent of Fig. 1. The authors argued that the dusty shell of Barnard's Loop was not optically thick enough to absorb all ionising photons from the Orion OB1 association. These photons leak through to a more distant wall that is identified in gas velocity and as a faint H α filament closer to the Galactic plane. Barnard's Loop is part of a closed bubble that may be driven by a supernova and expands inside the superbubble. The other shells (GS206-17+13, Orion Nebula) and the complex history of λ Orionis (possibly a supernova remnant cavity that has later been filled by the H II region around the OB association, Dolan & Mathieu 2002) led Ochsendorf et al. (2015) to change the simple picture of a single expanding object to a more complex combination along the line of sight of evaporating clouds and nested shells filled with X-ray emitting hot gas at different temperatures .
On the plane of the sky, the elongated bubble extends from beyond Orion up to the Eridanus Loop, but its overall orientation is still uncertain. Absorption line measurements (Welsh et al. 2005) and dust extinction data from Gaia and the Two Micron All-Sky Survey (2MASS; Lallement et al. 2018) support a close end towards Eridanus (southwest) and a far end beyond Orion (east), whereas models of a single shock front expanding in an exponentially stratified interstellar medium (ISM) fit the H α data better, with the Orion side being closer than the Eridanus side (Pon et al. 2016). The models cannot accommodate the composite structure of the superbubble, however, which has likely evolved in space and time from near to far along the blue stream of massive stars identified by Pellizza et al. (2005) and Bouy & Alves (2015) between us and the Orion clouds.
We have performed a multi-wavelength analysis of the gas and cosmic-ray content of the Eridanus side of the superbubble, which we present in a series of three papers. This paper focuses on the composite structure of atomic and molecular gas shells and on their relation to the hot ionised gas, the recombination fronts, and the magnetic field structure. The second paper (Joubaud et al., in prep., hereafter Paper II) explores the cosmic-ray flux that pervades the gas shells. The third paper (Joubaud et al., in prep.) studies the gas mass found at the transition between the atomic and molecular phases and in the CO-bright parts, as well as the evolution of the dust emission opacity (per gas nucleon) across the different gas phases. The third paper also studies the small molecular cloud of MBM 20, also called Lynd dark nebula (LDN) 1642, that might be compressed at the edge of the Local and/or Eridanus bubbles, near α J2000 = 68. • 8, δ J2000 = −14. • 2.
Here we investigate the composite structure of the Orion-Eridanus superbubble in the light of a new separation of the H I gas in position and velocity. We find new shells and study their relation to dust reddening fronts, H α and X-ray emissions, and the magnetic field orientation derived from dust-polarised emission. We also use the H I line information to estimate the magnetic field intensity in several parts of the superbubble.
The paper is structured as follows. Data are presented in Sect. 2. The separation of the H I and CO clouds and their relations in velocity, distance, and to known entities are discussed in Sect. 3. We study the X-ray emission and derive the emitting plasma properties in Sect 4. We probe the magnetic field orientation and strength in Sect. 5, before we conclude in Sect. 6. We present X-ray optical depth maps in Appendix A and the two methods we used to derive the angular dispersion of the magnetic field in Appendix B.

Data
We have analysed the eastern part of the superbubble towards the Eridanus constellation. It extends from 43 • to 78 • in right ascension and from −29 • to 21 • in declination, as shown in Fig. 1. We masked two 5 • wide areas on the western and eastern sides of the region to avoid complex gas distributions in the background. All maps are projected onto the same 0.25 • spaced Cartesian grid.

Ionised gas
Warm ionised gas is visible through H α emission. It is displayed in Fig. 1 using the data of Finkbeiner (2003). This is a composite map of the Virginia Tech Spectral line Survey (VTSS), the Southern H-Alpha Sky Survey Atlas (SHASSA), and the Wisconsin H-Alpha Mapper (WHAM). The velocity resolution is 12 km s −1 and the spatial resolution is 6 (full width at half maximum, FWHM).

H I and CO emission line data
We used the 16 2 resolution HI4PI survey (HI4PI Collaboration 2016), with a velocity resolution of 1.49 km s −1 in the local standard of rest (LSR). We selected velocities between −90 and 50 km s −1 to exclude the H I emission from the high-velocity clouds that lie in the hot Galactic corona far behind the local medium we are interested in Wakker et al. (2008).
In order to trace the molecular gas, we used the 8. 5 resolution 12 CO (J = 1−0) observations at 115 GHz from the momentmasked CfA CO survey (Dame et al. 2001;Dame & Thaddeus 2004). We completed this dataset with the CO observations of the MBM 20 cloud obtained with the Swedish-ESO Submillimetre Telescope (SEST) that were kindly provided by Russeil et al. (2003).

X-ray data
In order to study the hot-gas content of the superbubble, we used data from the Roentgen satellite (ROSAT) X-ray all-sky survey (Snowden et al. 1994). The spatial resolution is about 30 onaxis. We combined the original energy bands into three bands: the 0.25 keV band, spanning from 0.11 to 0.284 keV at 10% of peak response (R1+R2), the 0.75 keV band from 0.44 to 1.2 keV (R4+R5), and the 1.5 keV band from 0.7 to 2 keV (R6+R7). The exposure over the analysed region varies between 0 and 600 s, therefore we masked out the underexposed zones with exposure times below 80 s. They appear as white stripes in our intensity maps. Even though we used count rates that take striped variations in the exposure map into account, we caution that systematic biases remain in the survey calibration. They appear as striped enhancements or diminutions in count rates that cross the analysis region roughly parallel to the underexposed stripes. We re-sampled the maps onto our 0. • 25 grid and have compared the two low-energy count rate maps with those of Snowden et al. (1994).

Dust data
We used the 3D dust reddening maps of Green et al. (2018). These maps span three quarters of the sky (δ −30 • ) and are based on the stellar photometry of 800 million stars from Pan-STARRS 1 and 2MASS. The authors divided the sky into pixels containing a few hundred stars each. This results in a map with pixel sizes that vary from 3. 4 to 55 .
The polarisation data, that is, the total intensity and the Q and U Stokes parameters, come from the Planck 353 GHz 2018 polarisation data with an original angular resolution of 5 (Planck Collaboration XII 2019). To enhance the signal-tonoise ratio, we smoothed the maps and their covariances to 40 resolution (FWHM) and downgraded the healpix maps to N side = 512.

Velocity decomposition
Following the method developed by Planck Collaboration Int. XXVIII (2015), we decomposed the H I and CO velocity spectra into individual lines and used this information to identify and separate eight nearby cloud complexes that are coherent in position and in velocity. We thus fitted each H I or CO spectrum as a sum of lines with pseudo-Voigt profiles that combine a Gaussian and a Lorentzian curve that share the same mean and standard deviation in velocity. The Lorentzian part of the profile can adapt, if necessary, to extended line wings that are broadened by velocity gradients within the beam. The prior detection of line peaks and shoulders in each spectrum provided a limit A&A 631, A52 (2019) on the number of lines to be fitted as well as initial guesses for their central velocities. We improved the original method by taking the line information of the neighbouring spectra into account for each direction in order to better trace merging or fading lines. The fits yield small residuals between the modelled and observed spectra, but in order to preserve the total intensity that is recorded in each spectrum, we distributed the residuals among the fitted lines according to their relative strength in each channel.

H I and CO components
The 3D (right ascension, declination, and central velocity) distribution of the peak temperatures of the fitted lines showed that the data can be partitioned into several entities that are depicted in N H I column densities and in W CO intensities in Fig. 2. Some were detected in both H I and CO emission lines, such as the north rim, south loop, east shell, Cetus, North Taurus, and MBM 20. The others were only seen in the H I data, such as the west rim and Eridu cirrus. To construct these maps, we defined 3D boundaries in right ascension, declination, and velocity for each component. The velocity range for each cloud is presented in Table 1. We selected the fitted lines with central velocities falling within the appropriate velocity interval, depending on the (α, δ) direction, and we integrated their individual profiles in velocity. The resulting maps were then re-sampled into the 0. • 25 spaced Cartesian grid of the analysed region.
In order to investigate the effect of the unknown H I optical depth, we derived all the N H I maps for a set of uniform spin temperatures (100, 125, 200, 300, 400, 500, 600, 700, and 800 K) and for the optically thin case. Nguyen et al. (2018) have recently shown that a simple isothermal correction of the emission spectra with a uniform spin temperature reproduces the more precise H I column densities that are inferred from the combination of emission and absorption spectra quite well. Their analysis fully covers the range of column densities in the Eridanus region. All the figures and results presented here were derived for a spin temperature of 100 K. This choice comes from the best-fit models of the γ-ray data and of the dust optical depth at 353 GHz; they are detailed in Paper II. This choice has no effect on cloud separation. Compared with the optically thin case, the correction increases the highest column densities (>1.5 × 10 21 cm −2 ) of the thickest cloud of the north rim by 40%. The other clouds are more optically thin.
As we discuss in this paper, some structures are likely associated with the superbubble (north rim, south loop, east shell, and west rim), while the others are foreground (MBM 20), background (Cetus, North Taurus) or side (Eridu) clouds. The north rim component gathers elements that have previously been identified as independent features, such as the MBM 18 (L1569) molecular cloud at a distance of 155 +3 −3 pc (Zucker et al. 2019) and the G203-37 cloud towards (66 • , −9 • ) that Snowden et al. (1995) placed midway through the hot superbubble cavity. Other MBM clouds of interest are listed in Table 1. We did not push the partition of the north rim farther because the lines are confused in space and velocity, despite the sparse decomposition of the spectra into individual lines. We reached the same conclusion using another decomposition method called the regularized optimization for hyper-spectral analysis (ROHSA 1 ). The analysis kindly provided by Antoine Marchal recovered the same structures for the east shell and west rim. This supports our separation. However, it also struggled to identify individual entities in the complex distribution in the northern part. We therefore preferred to keep the north rim as a single component for our analyses. 1 https://antoinemarchal.github.io/ROHSA/index.html Part of the analysed region overlaps the anticentre region studied by Remy et al. (2017), who used a similar method to separate gas clouds. We recovered the same structures in velocity and space in the overlap region: their North Taurus and Cetus components, with minor differences due to their use of data with higher velocity resolution from the Galactic Arecibo L-Band feed array (GALFA) in those directions. Their South Taurus complex mainly contributes to the north rim in our analysis. Their main Taurus complex only marginally overlaps with our region. The bulk of this cloud lies beyond our analysis boundary and is not discussed here.

Relative motions in the superbubble
The distribution in position and velocity of the H I line cores that result from the line decomposition allows us to resolve the relative motions that subtend this complex environment. In the original data cubes, the bulk motions are often  buried in the overlap of the broad wings of the H I lines. We display slices through the H I line distributions for specific velocities in Fig. 3. The central velocities of the H I lines associated with the superbubble range from −15 to +18 km s −1 . This range compares with the [−25, +5] km s −1 range derived by Reynolds & Ogden (1979) from the H α emission. The map at +9 km s −1 exhibits a coherent circular ring that suggests that we intercept the outer rim of the expanding shell tangentially along these lines of sight. Most of the gas at higher (receding) velocities is indeed seen within the boundary of the outer rim. The slice at +4 km s −1 shows another circular ring in the south that corresponds to the rim of the south loop. The expansion velocities of these two rings are in the plane of sky because we see the rims tangentially. The radial line velocities thus indicate that the bulk of the outer rim and of the south loop recede at +9 and +4 km s −1 , respectively, with respect to us. This means that the south loop approaches us with respect to the rest frame of the superbubble. The slice at −9 km s −1 shows the east shell, which approaches fast. We give evidence in Sect. 4 that it moves inside the superbubble.
When we exclude the internal motion of the east shell to follow the outer expansion, the H I line range is reduced to [−5, +19] km s −1 . With a rest frame velocity of +9 km s −1 , the line data suggest an expansion velocity of about 10-15 km s −1 , which is lower than the 40 km s −1 inferred by Brown et al. (1995) from the full extent of the H I spectra, including the line wings. It is consistent with the 15-23 km s −1 velocity found in H α (Reynolds & Ogden 1979). It also compares with the low velocity of 8 km s −1 expected from simulations (Kim et al. 2017) for a superbubble that expands in a medium with a mid-plane gas density of 1 cm −3 and a scale height of 104 pc that is powered by a supernova rate of 1 per Myr; this is close to the rate of Orion-Eridanus superbubble (Voss et al. 2010). The simulations also show that the compressed swept-up gas behind the shock expands at slightly lower velocities than the shock itself and than the hot phase.

Dust reddening distances
In order to study the 3D gas distribution in the superbubble and to substantiate the cloud separation in velocity, we used the 3D reddening maps of Green et al. (2018). They inferred a reddening profile as a function of distance modulus, µ, for each direction in the sky using the photometric surveys of Pan-STARRS 1 and 2MASS in a Bayesian approach. We located gas concentrations in distance in these profiles by searching for reddening fronts that can be identified as sharp increases in the amount of dust reddening per unit length. To do so, we took the derivative of their bestfit profiles for each direction in the sky. We display the derivative maps in Fig. 4 for three distances at 200, 320, and 400 pc. The west rim and MBM 20 clouds are clearly identified around 200 pc, and the north rim appears between 200 and 320 pc. The CO-bright molecular parts of the east shell are visible at 320 pc. The G203-37 cloud appears between 320 and 400 pc.
The 3D partitioning of the dust fronts was not induced by the Bayesian priors that were used to build the reddening profiles. Green et al. (2018) have used broad-band photometric measurements of each star and have computed a probability distribution over the stellar distance and foreground dust column. They constrained their prior on the Schlegel et al. (1998) reddening map, which is a scaled version of an optical depth map of the dust thermal emission integrated along the lines of sight, without information along the radial dimension. The spatial correlation we find between the dust maps in distance and the morphology of the clouds we isolated strengthens our confidence that the different complexes are separated in velocity as well as in distance.
We constrained the distance range to each cloud complex using the reddening front information towards several directions in each complex. Their locations are marked in Fig. 5. For each direction, we selected derivative profiles within a radius of 1 • , and we averaged those whose peak dE(g−r) dµ value exceeded 80% of the maximum. The results are shown in Fig. 5. Distance ranges based on these profiles are given in Table 1 together with the velocity ranges observed in the H I data. Table 1 also lists distances to MBM molecular clouds that appear to be associated in space and velocity with our cloud complexes. Their distances have been derived by Zucker et al. (2019) using the same reddening method, but choosing individual stars. Lallement et al. (2018) constructed another 3D dust map using inversion methods and the stellar surveys of Gaia, 2MASS, and the APO Galactic Evolution Experiment (APOGEE) DR14. Three vertical cuts across A&A 631, A52 (2019) Fig. 5. Mean reddening per distance unit profiles towards several directions for each H I cloud. The mean is computed for each direction over selected profiles in a circle with a radius of 1 • . The map shows the pointed directions overlaid on the H I column density. The different markers correspond to the different cloud complexes, and the colours show the different directions for a given cloud. The cloud labels are north rim (EriN), west rim (EriW), east shell (EriE), south loop (EriS), and North Taurus (TauN). For clarity, we divided the north rim into its northern part (EriN1, upwards-pointing triangle), which marks the superbubble boundary, and its southern part (EriN2, downwards-pointing triangle), which overlaps the MBM 18 cloud. The white markers on the map show the MBM clouds. Their shape corresponds to the complex they are associated with. Their distances were derived by Zucker et al. (2019) and are reported in the plots with grey dotted lines.
the Galactic disc towards the Galactic longitudes of 192 • , 205 • , and 212 • yield distances to our gas complexes that are listed in Table 1.
The different measurements in Table 1 indicate that most of the clouds associated with the superbubble, such as the north rim and west rim, the east shell, south loop, and MBM 20, lie at distances ranging between 150 and 250 pc. The distances and velocities we obtain for the different complexes are consistent with those found with NaI and CaII interstellar absorption lines towards early-type stars throughout the region (Welsh et al. 2005). We note, however, that distances based on the reddening information from Green et al. (2018) tend to yield higher values than towards individual stars. This is partly due to their binning in space and in log(µ), which was adapted to large-scale 3D mapping in the Galaxy, and to the lack of stars at high Galactic latitudes. The lower bound of our distance estimates should also be considered with care as they approach the shortest range allowing large enough stellar samples to infer the average reddening and to recognise the shape of the cloud. Adding stars from Gaia alleviates the biases. All sets of values nevertheless give weight to a superbubble orientation with its close end towards Eridanus.
A series of MBM clouds (MBM 16, 106, 107, 108, and 109) is associated with the north rim and more specifically with the outer rim because of their location and velocity. They gather at distances between 140 and 170 pc. The MBM 18 cloud, at a comparable distance, is also associated with the north rim in H I and CO. Together they suggest that the broad north rim complex combines gas from the outer rim seen edge-on and face-on (see the middle right panel of Fig. 3). We therefore conclude from the synergy in distance and velocity of the north rim, south loop, and west rim that they represent the front face of the superbubble with respect to us. Figure 4 also shows the G203-37 cloud at a distance between 300 and 400 pc, in agreement with Snowden et al. (1995). These authors had placed it midway through the hot superbubble cavity.
The MBM 22 cloud is seen in a direction close to the rim of the south loop, and its velocity is consistent with the south loop range, but its distance of 266 +30 −20 pc places it behind the south loop. This is corroborated by the maps in Fig. 4, which show that the south loop is most conspicuous in the left (200 pc) panel and fades out at 320 pc, whereas MBM 22 appears at 320 pc and fades out at a larger distance. MBM 9 and the clouds of Cetus and North Taurus also lie outside and behind the superbubble, at distances between 250 and 400 pc. We note that the Cetus and North Taurus clouds are part of larger complexes that extend well beyond the western side of the analysed region (Remy et al. 2017).

MBM 20 and the edge of the superbubble
Inspection of relative distances in Figs. 4 and 5 suggests that MBM 20 and the closest clouds of the south loop may be interacting with the Local Bubble at about 150 pc in the picture of Burrows et al. (1993). This is consistent with the velocity configuration discussed above, where the south loop approaches us compared to the rest of the superbubble.
A clear H I and CO component in our separation corresponds to the well-studied MBM 20 (L1642) cloud (Russeil et al. 2003).    (Lallement et al. 2018) places it near the boundary between the superbubble and the Local Bubble. Snowden et al. (1995) used X-ray absorption to place MBM 20 inside the Local Bubble, thus pushing the closest superbubble wall beyond 140 pc. The cometary tail of MBM 20 is visible in H I. It points to the north-east of the compact CO cloud (see Fig. 2) and extends over 10 pc (at 140 pc). It presents a velocity gradient whose tail tip moves faster away from us than the molecular head, in agreement with having been swept up by the expanding medium of the Local Bubble rather than by the south loop. MBM 20 may therefore lie just in front of the edge of the superbubble. We note that despite similar orientations and some overlap in direction with the east shell, the two clouds are distinct and have different velocities, as is shown in Table 1.

Cloud relations to the H α filaments
The H I and CO structures we isolated can be compared to the H α intensity map. H α line emission at 656.3 nm arises from the recombination of an electron onto ionised hydrogen, so that the lines trace the ionised gas at the interface with colder atomic hydrogen, where the ionisation fraction is still high, but the density is also high enough for recombination to occur. Figure 6 displays H α line intensity maps integrated over [−60, −11] and [−11, 50] km s −1 . We overlaid H I contours of specific clouds. The 12 km s −1 resolution of the WHAM H α data prevents detailed velocity comparisons between the H I, CO, and H α lines, but it can still capture interesting relations. As the neutral H I structures surround the recombining shells with a small angular offset, Fig. 6 outlines a mosaic of boundaries between hot sub-bubble interiors and their colder outer shells, which are compressed by the hot-gas expansion (see Figs. 3 and 4 of Krause et al. 2013).
The east shell is related to a specific H α filament. Both structures move inside the superbubble towards us (negative velocities) at a distance constrained by the dust to be about 300 pc. The southern part of the south loop might be associated with the H α filament called Arc C (Johnson 1978). The round structure of the loop in H I and H α and its radius of ∼35 pc at a distance of 200 pc are reminiscent of an old supernova remnant, but they may also reflect global oscillations of the hot gas in the superbubble. Such oscillations are induced by off-centred supernova explosions in a hot cavity (Krause et al. 2014). We further discuss the origin of the loop in the next sections. The west rim is related to the Arc B filament in H α (Johnson 1978). We see an elongated shock wall almost edge-on, with velocity excursions of −10 to +7 km s −1 about the rest velocity of +9 km s −1 . The closest part approaches us and is superimposed on receding parts that lie farther away. The H α emission from the rest of the outer rim is likely absorbed by the high gas column densities pertaining to the north rim (typical visual extinctions of 0.4-1.8 magnitudes A52, page 7 of 19 A&A 631, A52 (2019) (top) and [−11, 50] km s −1 (bottom), overlaid with N H I column density and W CO emission contours outlining the different cloud components: the east shell in H I (red, 1.5 × 10 20 cm −2 ), the west rim in H I (cyan, 4 × 10 20 cm −2 ), the south loop in H I (green, 2.5 × 10 20 cm −2 ), and the north rim in CO (yellow dashed line, 1 K km s −1 ). The N H I and W CO maps are smoothed with a Gaussian kernel of 0. • 28 to display the contours. A specific north rim contour is shown in H I for central line velocities of +9 km s −1 (dark blue dash-dotted line, 18 K). It was smoothed with a Gaussian kernel of 1. • 5 before the contours were computed.
up to a distance of 300 pc along the northern part of the outer rim; this is outlined in dotted contours in Fig. 6. Towards the west rim, the visual extinctions range from 0.1 to 0.4 magnitudes).
The brightest H α arc, Arc A, extends vertically in the centre of the positive-velocity map in Fig. 6. Its relation to the superbubble is still debated; see Pon et al. (2014) for a discussion that takes the proper motion and radial velocities, as well as X-ray and H α emission and absorption studies, into account. Based on the width and intensities of the H α lines, these authors argued that the H α filaments are either elongated sheets seen edgeon or objects that are out of equilibrium as a result of ionised clouds that are compressed by shocks that currently recombine. While Arc B seems to comply with the edge-on sheet hypothesis, the situation is not clear for Arc A because of the confusion along the line of sight, the spread of H α and H I emission over a wide range of velocities, and absorption from atomic and dark neutral gas (see Paper II) and from molecular gas in MBM 18. As reported by Pon et al. (2014), the coincidence in space and velocity between Arc A and MBM 18 tends to favour an association of the arc with the superbubble. Figure 6 shows that the northern part of Arc A nicely follows an arc of molecular clouds belonging to the north rim component and that Arc A is partially absorbed by the dense gas in MBM 18 near α J2000 = 60. • 6, δ J2000 = 1. • 3. The dust associated with the molecular arc is visible at close distance in Fig. 4, but not at 400 pc, nor beyond. The distance to Arc A is thus constrained to be between that of MBM 18 (155 +3 −3 pc, Zucker et al. 2019) and an upper bound of about 400 pc. Arc A therefore likely represents an inner shell of the superbubble seen edge-on, as for the east shell.

X-ray emission
The Eridanus X-ray enhancement arising from diffuse hot gas has been extensively studied to constrain the properties of the hot superbubble interior and the distribution of foreground gas seen in absorption against the bright X-rays (Burrows et al. 1993;Snowden et al. 1995;Heiles et al. 1999). We now compare the spatial distribution of the different clouds we isolated with that of the X-rays recorded in three energy bands with ROSAT, and we take advantage of our study of the total gas column densities in the atomic, dark, and molecular gas phases (see Paper II) to quantify the X-ray absorption. Figure 7 displays the ROSAT maps recorded in the 0.11-0.284 keV (R1+R2), 0.44-1.21 keV (R4+R5), and 0.73-2.04 keV (R6+R7) energy bands (Snowden et al. 1994). We refer to these bands below as the 0.25, 0.75, and 1.5 keV bands, respectively. Two main features have originally been identified (Burrows et al. 1993). The first, called Eridanus X-ray enhancement 1 (EXE1), appears in the 0.75 and 1.5 keV maps as a crescent bounded to the west by the west rim clouds, to the south by the east shell, and to the east, it extends towards Orion. Its northern extent is unknown because it is obviously heavily absorbed, even at 1.5 keV, by the 1.4 ×10 5 solar masses of gas associated with the north rim (see Fig. 7). EXE1 was associated with the interior of the superbubble by Burrows et al. (1993). The brightest part of the crescent at 0.75 keV is located between the H α arcs A and B. This region was found to be a small independent expanding shell in H α by Heiles et al. (1999).
The second main feature (EXE2) is visible only at 0.25 keV, below δ = −10 • . This soft excess is obviously barred by absorption from the gas in the south loop and by the compact MBM 20 cloud, which creates a marked dip in X-rays around α = 68. • 8 and δ = − 14. • 2. The east shell, however, leaves no obvious absorption imprint. The bright part of EXE2 that partially overlaps the south loop was first interpreted as hot gas originating from a supernova remnant or a stellar wind bubble (Burrows et al. 1993), but the absence of an OB star in this direction led Heiles et al. (1999)   propose that hot gas leaks out of EXE1 at the rear of the superbubble through a break in its shell and that the gas cools in the external ISM. Because the soft emission extends below −25 • in declination, Snowden et al. (1995) instead interpreted EXE2 as soft background emission arising from the million-degree gas in the Galactic halo and being shaped into a roundish feature by intervening gas absorption. We revisit the properties of EXE1 and EXE2 below.

Subtracting foreground and background emission
We followed the simple slab model of Snowden et al. (1995) and modelled the observed X-ray intensity in each energy band and direction as I obs (α, δ) = I Local Bubble + I Eri (α, δ) e −τ X (α,δ) + I halo e −τ tot X (α,δ) , (1) where I Local Bubble , I Eri , and I halo represent the X-ray intensities arising from the Local Bubble, the Orion-Eridanus superbubble, and the halo, respectively. The foreground and background emissions from the Local Bubble and the halo are assumed to be uniform. τ X denotes the X-ray optical depth of the absorbing gas located in front of the superbubble, and τ tot X denotes the total optical depth along the line of sight.
In order to compute the X-ray optical depths, we used the N H gas column densities inferred in the different gas phases A&A 631, A52 (2019) and different cloud complexes from our combined H I, CO, dust, and γ-ray study, following the same method to map and quantify the amount of gas as applied earlier in other nearby regions (see Paper II, Planck Collaboration Int. XXVIII 2015; Remy et al. 2017). Snowden et al. (1995) have instead used the dust emission map at 100 µm to trace N H , but they did not correct the map for biases that have since then been found to be significant (e.g. spatial variations in dust temperature and spectral index (Planck Collaboration XI 2014) and the non-linear increase of the dust emissivity per gas nucleon with increasing N H (Planck Collaboration XI 2014; Planck Collaboration Int. XXVIII 2015; Remy et al. 2017). We derived the optical depths from our column density maps using with E X the X-ray band and σ the energy-band-averaged photoelectric absorption cross section from Snowden et al. (1994). For the 0.25 and 0.75 keV bands, we took an incident thermal spectrum with a temperature of 1.3 MK close to the mean colour temperature observed towards the Eridanus X-ray excesses (see below). For the 1.5 keV band, we assumed a power-law incident spectrum with an index of −2. We verified that the choice of temperature and spectral index around these values has a small effect on the result. The optical depth maps obtained for the total gas column densities τ tot X are presented in Fig. A.1 for the three energy bands. We find that the north rim is optically thick in the lower energy bands and is still partially thick (0.5 τ tot X 1) at 1.5 keV, which explains the northern shadows in all plots of Fig. 7. Inspection of the optical depth maps for individual cloud complexes indicates that the west rim is thick at 0.25 keV, nearly so at 0.75 keV, and mostly thin (0.1 τ X 0.2) at high energy. The south loop is thick at low energy and mostly thin in the two upper bands. The east shell is thick at low energy (1 τ X 2.5), partially thick at 0.75 keV (0.1 τ X 0.4), and thin at high energy. The molecular head of MBM 20 is optically thick at low energy, nearly so at medium energy, but partially transparent at high energy.
In order to estimate upper limits to the emission from the Local Bubble in the three ROSAT bands, we studied the distribution of X-ray intensities versus gas column density towards a set of heavily absorbed and nearby clouds: towards MBM 20 (which is likely located at the interface between the Local Bubble and the Eridanus superbubble), MBM 16, MBM 109, and MBM 18. To limit the effect of noise fluctuations towards these faint regions, we calculated the minimum intensities averaged over the high-N Htot pixels in each region, at N H > 7 × 10 20 cm −2 in the low band and N H > 1.5 × 10 21 cm −2 in the upper bands, and we excluded pixels with X-ray point sources. We find an upper limit to the Local Bubble intensity of 4.0 × 10 3 cts s −1 sr −1 at 0.25 keV, which coincides with the estimate of Snowden et al. (1995). We find upper limits of 0.8 × 10 3 cts s −1 sr −1 at 0.75 keV and 1.1 × 10 3 cts s −1 sr −1 at 1.5 keV, which are consistent towards the four clouds.
In order to estimate the background halo intensities in directions away from the superbubble, we corrected the observed X-ray maps for the Local Bubble emission and for the total gas absorption using I halo = (I obs − I Local Bubble ) e +τ tot X . We averaged the intensities over the low-absorption region at 50 • ≤ α ≤ 60 • and δ ≤ −25 • to avoid the underexposed survey stripes. We obtained intensities of 22.0 × 10 3 cts s −1 sr −1 at 0.25 keV, 0.5 × 10 3 cts s −1 sr −1 at 0.75 keV, and 0.2 × 10 3 cts s −1 sr −1 at 1.5 keV. The bottom row of Fig. 7 shows the residual intensities obtained in the three ROSAT bands after the uniform Local Bubble emission and the absorbed halo background were removed from the original data, These maps show the distributions of X-rays that are likely produced by the hot gas that pervades the superbubble, but the spatial distributions are still modulated by photoelectric absorption by internal and/or foreground clouds. The residual intensities in the two upper bands are consistent with zero across the reference region we used to estimate the halo background. As our estimate of the halo intensity at 0.25 keV is 45% lower than that of Snowden et al. (1995), the data are not as over-subtracted at low declinations as in their analysis. Only a few patches of small negative residuals remain at δ < −20 • . They may be due to anisotropies in the Local Bubble and/or in the halo at low X-ray energies. Figure 8 presents a colour image of the superbubble emission that remains after the Local Bubble foreground and halo background are subtracted according to Eq. (3). It reveals several colour gradients in the X-ray emission: a marked hardening to the east towards Orion, and a more moderate hardening in the expanding region that is enclosed between the H α arcs A and B. As pictured by Ochsendorf et al. (2015), the hardening towards Orion traces hotter gas that is energised near Orion and flows into the colder Eridanus region. This hot gas may come from the supernova that is responsible for Barnard's Loop, the eastern side of which is visible in H α , but the western side of which has long been in the process of merging with the superbubble hot gas. The two bright emission regions seen at δ < −20 • in the original 0.25 keV map (top left plot of Fig. 7) are well accounted for by the Local Bubble and halo intensities, whereas significant emission remains visible inside the south loop after the foreground and background emissions are removed. In addition, the east shell that bounds the hard EXE1 emission leaves no absorption shadow at 0.25 keV despite its strong optical depth in this band. Together, these facts suggest that hot gas inside the superbubble spans distances beyond the east shell in hard X-rays and in front of the east shell towards the south loop in soft X-rays. By correcting the soft "Eridanus-born" X-ray intensity (bottom left plot of Fig. 7) for the amount of absorption expected from south loop gas, in other words, by plotting (I obs − I Local Bubble − I halo e −τ tot X )/e −τ X South , where τ X South is the optical depth of the south loop, we verified that the X-ray gap that separates the south loop from the hard emission shining between arcs A and B can be fully accounted for by absorption in the south loop rim. The in situ emission may thus be continuous.

Hot-gas properties
We studied the average properties of the hot X-ray emitting gas in a sample of ten directions spanning different parts of the EXE1 and south loop regions, as displayed in Fig. 9. The 2. • 25 width of each square area in the sample is close to the field of view of the ROSAT position sensitive proportional counter (PSPC) detector, whose diameter is 2 • . For each area, we calculated the mean X-ray intensity that remains in the three energy bands after the Local Bubble foreground and halo background are subtracted (according to Eq. (3)). We considered two hardness ratios: the 0.25-0.75 keV and the 0.75-1.5 keV intensity ratios. We also calculated the mean gas column density, N H , in each area, using specific cloud components: only the nearby south loop for the first three areas that sample the soft EXE2 emission without an indication of strong absorption by the optically thick east shell; using the south loop and the east shell for regions 4-8 because they sample the hard EXE1 emission behind the east shell; and using the whole gas for regions 9 and 10 because they might intercept the faint edge of the north rim. The resulting column densities are given in Fig. 10. We used XSPEC v12.9 (Arnaud 1996) and the ROSAT PSPC C response (appropriate for the survey) to model the expected flux in the three energy bands from the Mewe-Kaastra-Liedahl (mekal) model of thermal emission arising from a tenuous hot plasma with solar abundances. The mekal model includes bremsstrahlung radiation and the relevant atomic shell physics for line emissions (Mewe et al. 1985;Liedahl et al. 1995). We selected the Tübingen-Boulder (tbabs) model of interstellar absorption and scanned temperatures from 0.01 to 10 keV in log steps of 0.1. The absorption column was fixed to the mean value in each area. The free parameters of the model are thus the gas temperature, T , and the plasma emission measure (i.e. the volume integral n e n H dV of the hydrogen volume density, n H , times the electron volume density, n e , the latter becoming 1.2 n 2 H dV for a fully ionised plasma with 10% helium by number).
We considered a single-temperature model and a mixed model with two temperatures. In the former case, hardness variations across the field would primarily be due to absorption features. The latter case would imply two pockets of hot gas at different temperatures that partially overlap in direction, but are separated in distance and position relative to the east shell. The model can accommodate more complex interleaved geometries along the lines of sight as long as the absorbing N H slab stands in front of both hot gases. We first compared the modelled and observed hardness ratios in each area and found that a single temperature cannot account for the three-band colour distribution in any of the ten areas. The temperatures inferred from the 0.25 to 0.75 keV intensity ratios vary from 1.1 to 1.5 MK around a mean of 1.3 MK, which corroborates our choice of temperature for the derivation of optical depth maps in the former section and agrees with the former estimates of Guo et al. (1995). In particular, the temperature of 1.5 MK we find towards field 7 is consistent with the values of 1.75 +0.35 −0.34 and 1.64 +0.42 −0.17 found by Guo et al. (1995) in directions enclosed in this field. In all fields, the inferred temperatures do not predict enough emission in the 1.5 keV band, however.
In order to solve the two-temperature model in each area (which requires four parameters against three data points), we furthermore assumed that the two plasma pockets along the line of sight are in pressure equilibrium, in agreement with simulations (Kim et al. 2017). This is substantiated by the 0.7 Myr sound crossing time of the 80 pc wide region that emits X-rays (for 1.3 MK). We adopted a typical 80 pc size for the pockets in view of the 70 pc diameter of the south loop at a distance of 200 pc and of the apparent EXE1 size of 100 pc at the further distance of ∼250 pc of the east shell. The temperature, gas density, and mass we found for the two plasmas in each area are displayed in Fig. 10, together with their common pressure and the total thermal energy content of the two plasmas. The mass and energy estimates assume the same spherical geometry, 80 pc in diameter, for the two plasmas. We did not attempt to propagate the original count rate uncertainties to the results shown in Fig. 10 because of the striped calibration problems across the field, the simplicity of a uniform foreground and background subtraction, and the simplicity of a single slab that absorbs the emission from the two plasmas (as opposed to a more complex distribution of the material in front, in between, and inside the emitting regions). The model uncertainties dominate those in the X-ray data.  Fig. 10. Temperatures, gas densities and masses, and common pressure obtained from the superposition of two hot plasma bubbles (mekal models, see Sect. 4.2), 80 pc in diameter, towards each of the areas sampled in Fig. 9. Blue circles (left scale) and black crosses (right scale) mark the hotter and colder media, respectively. The lower right plot gives the mean N H column density we used for absorption in each area.
We nevertheless find a limited dispersion in the properties of the two plasmas across the different areas. Temperatures range from 3 to 9 MK in the hotter plasma and between 0.3 and 1.2 MK in the colder plasma. They compare well with the asymptotic values of 1.6-2.1 MK that were found in simulations for a superbubble resembling that of Orion-Eridanus (model n1-t1 of Kim et al. 2017). The gas enclosed between the H α arcs A and B (areas 8-10) appears to be globally hotter and at a higher pressure than in the rest of the region, which supports the Heiles et al. (1999) view of a younger and faster-expanding zone between the arcs. The remaining data points are roughly consistent with two 80 pc wide bubbles that encompass all of them.
The low plasma masses we find indicate that most of the hot gas inside the superbubble has already cooled down radiatively. According to simulations (Krause et al. 2014), soft X-rays trace the current thermal energy content of a superbubble, which we find to be about (1−2.4) × 10 50 ergs in Eridanus. They do not trace the cumulative energy injected during its lifetime. The thermal energy in the hot gas represents only a few percent of the current kinetic energy of the expanding H I gas (3.7 × 10 51 ergs, Brown et al. 1995). The latter must be revised downward to 5 × 10 50 ergs for an expansion velocity of 15 km s −1 . If this is confirmed, the thermal energy amounts to 20-50% of the current kinetic energy in the expanding gas. Pon et al. (2016) have modelled the expansion of the overall superbubble assuming an ellipsoidal geometry. Expansion in an exponentially stratified ISM can match the H α data for a variety of inclinations onto the Galactic plane (Pon et al. 2016), but all models require an internal (uniform) pressure of about 10 4 cm −3 K and an internal temperature of 3-4 MK. We find slightly higher values in the X-ray emitting hot gas across the whole Eridanus end of the superbubble, with pressures in the range of (3.2−7.8) × 10 4 cm −3 K. The pressure at the closest end of the superbubble towards the south loop exceeds that in the Local Bubble, which is filled with 1 MK gas with a pressure of 10 4 cm −3 K (Puspitarini et al. 2014). We formally find pressures that are at least three times higher in the south loop, but we recall that the measurements in the Local Bubble and in south loop have large uncertainties. An overpressure is consistent with the expansion of the Eridanus tip of the superbubble moving towards us; it will push against the Local Bubble if they are in contact.
The gas in both phases exhibits ten times higher volume densities than were found in the simulations of Kim et al. (2017) mentioned above. The 80 pc size of the X-ray emitting regions implies gas column densities of (0.6−1.8) × 10 18 cm −2 in the hotter plasma and (0.4−2.4) × 10 19 cm −2 in the colder plasma. These values compare well with the superbubble simulations of Krause et al. (2014) during the short stage when a recent supernova shock approaches the outer shell and drives a peak in X-ray luminosity. At other times, that is, during the megayear-long period that separates supernovae in Eridanus (Voss et al. 2010), the simulated column densities are ten times lower than what we measure in Eridanus. However, Krause et al. (2014) noted that their simulated X-ray luminosity fades too rapidly after a supernova explosion compared to various observations. A longer delay after the last supernova might therefore still be consistent with the Eridanus X-ray data. Half a million years after an offcentred supernova, the simulations show that the hot gas globally oscillates across the superbubble. The onset of such oscillations could explain the need for a multi-phase hot plasma towards all the sampled areas as well as the presence of hot gas in the south loop despite the lack of young massive stars. It is so far unclear whether the three hot regions of the south loop, the enclosure between arcs A and B, and the rest of EXE1 above the east shell represent distinct bubbles that are in the process of merging, or whether they represent the dynamical response of the overall superbubble plasma to the sequence of stellar winds and supernovae that have occurred along the stream of blue stars that extends to Orion (Pellizza et al. 2005;Bouy & Alves 2015). A joint modelling of the H α and X-ray emissions is required to follow the gas cooling inside the superbubble in order to distinguish the origin of the hardness variations across the different zones.

Magnetic field in the superbubble
The magnetic field in the Orion-Eridanus superbubble has recently been analysed by Soler et al. (2018, Soler18 for short,) using Planck 2015 polarisation observations at 353 GHz. Their study highlighted the interaction between the superbubble and the magnetic field, as revealed by the strong polarisation fraction and the low dispersion in polarisation angle along the outer shell. Figure 11 shows the orientation of B sky , the planeof-sky magnetic field that overlies the H I column density for negative-velocity structures with −30 < v LSR < −4 km s −1 and for positive-velocity structures with −4 < v LSR < 25 km s −1 . The bottom panel shows that the magnetic field, frozen in the gas, has been swept up by the superbubble expansion. The orientation of B sky indeed follows the shape of the shell along the west rim and north rim. A smaller magnetic loop corresponds to the south loop. Inside the superbubble, the magnetic field appears to be more disordered, as expected from the activity of stellar winds and past supernovae. Outside the superbubble, at higher latitudes below the Galactic plane (e.g. in the lower right corner of Fig. 11), the projected field lines appear to be smoothly ordered and nearly parallel to the Galactic plane. The expansion of the shock front along the west rim and part of the south loop is therefore in a favourable configuration to efficiently compress the external magnetic field.

B sky estimation method
Soler18 probed the magnetic field strength in three directions along the west rim by combining B sky estimates they derived from polarisation data and the Davis-Chandrasekhar-Fermi method (hereafter the DCF method, Davis 1951;Chandrasekhar & Fermi 1953) together with estimates of the B los field strength along the line of sight obtained by Zeeman splitting (Heiles 1989). The three directions are visible in the bottom panel of Fig. 11, together with other Zeeman measurements from Heiles (1989). They found B sky values of 25-87 µG that far exceed both B los and the average field strength found in cold atomic clouds in the solar neighbourhood (∼6 µG, Heiles & Crutcher 2005). The DCF method, however, tends to overestimate the field strength because of line-of-sight and beam averaging. A correction factor of about ∼0.5 was advised by Ostriker et al. (2001). It was not used in Soler18 because the probed regions had different physical conditions than in the MHD simulations that were employed to estimate the correction. Instead of applying this factor, we used the modified DCF method described by Cho & Yoo (2016) to probe B sky in the superbubble. Their method takes into account the reduction in angular dispersion of the magnetic field orientation that is due to averaging effects, which in turn arise from the pile-up of independent eddies along the line of sight. Assuming that the dispersion in B sky orientation is small and is due to isotropic Alfvénic turbulence, their expression to derive the plane-of-sky magnetic field strength is  11.2 ± 1.9 North rim (72.5, −13.2) 2.6 ± 0.9 West rim (b) (49.6, −8.3) 0.69 3 ± 1 5.5 ± 0.1 8.0 ± 0.8 4.8 ± 1.0 −6.6 ± 0.6 West rim (b) (51.2, −10.3) 0.60 3 ± 1 6.6 ± 0.1 9.7 ± 0.8 3.3 ± 0.7 −11.8 ± 1.3 West rim (b,c) (52.3, −12.3) 1.15 3 ± 1 6.4 ± 0.1 9.4 ± 0.8 6.1 ± 1.2 −9.5 ± 0.7 West rim (55.1, −17.2) 0.49 5 ± 2 5.3 ± 0.1 8.0 ± 0.8 4.5 ± 0.9 West rim (50.9, −13.5) 0.82 4 ± 2 2.9 ± 0. where ρ is the gas mass density, δV c is the standard deviation of H I emission line velocities, and ς ψ is the angular dispersion of the local magnetic-field orientations. ξ is a correction factor derived from simulations with values between ∼0.7 and ∼1. The results presented here are calculated for ξ = 0.7 because the DCF method tends to overestimate field strengths. We applied this method towards several directions in the superbubble. They were chosen to exhibit a single dominant velocity structure in the probed area for a reliable evaluation of δV c . This criterion rejected B sky estimates in broad regions of the north rim or towards Cetus where two velocity components gather along the lines of sight. We still found two directions that intercept the approaching front of the north rim without much contamination from other structures. For all directions in our sample, we verified that more than 70% of the dust thermal emission arises from the cloud of interest to ensure that the measured dispersion in magnetic field orientations strongly relates to that cloud and to limit the overestimation of the field strength because the Planck data are averaged over the entire line of sight. Because of this selection and because of velocity information of the H I lines, we can associate each magnetic field estimate with a single gas complex, as presented in Sect. 3. The probed directions were also chosen to have a signal-to-noise ratio (S/N) higher than 3 in polarised intensity, P/σ P ≥ 3, to derive the polarisation angle dispersion ς ψ . This second criterion rejected estimates of B sky in the east shell and in one of the directions analysed by Soler18 towards (α, δ) = (52.3, −13.3), which they considered at a lower angular resolution to maintain P/σ P ≥ 3. The resulting B sky value is displayed in Table 2 but is not used in the following analysis. We note that the measurements of Soler18 correspond to regions where 25-45% of the dust emission comes from the north rim and not from the west rim that they wished to probe. We therefore added directions towards denser parts of the west rim.
The H I line velocities come from the decomposition analysis presented in Sect. 3: they are the central velocities of individual pseudo-Voigt lines fitted against the H I spectra. Soler18 estimated B sky in circles with 3 • diameters. In order to avoid velocity gradients or crowding across the test regions, we reduced the diameters to 2 • . We selected the H I lines pertaining to the cloud complex to be probed, with peak brightness temperatures above 7 K in order to avoid noisy and poorly detected lines. We computed the standard deviation of the resulting distribution, δV c , weighted by the line brightness temperatures. The results are displayed in Table 2.
We used two methods to derive the angular dispersion of the local magnetic field orientations, ς ψ . The first directly uses the Stokes parameters as defined in Planck Collaboration Int. XXXV (2016). The second relies on the structure function of the polarisation angles. Both methods yield consistent results (see Fig. A.2), and they are detailed in Appendix B.
We derived the average mass volume density ρ from the H I column density assuming a line-of-sight depth. The south loop forms a nearly complete ring with a thickness of ∼10 pc and a radius of 36 pc at a distance of 200 pc. The west rim cloud lies along a roughly spherical cap of comparable thickness and with a radius of 90 pc at a distance of 200 pc. Upper limits to the line-ofsight depth across these shells are given by the crossing lengths that graze the inner radius. They correspond to 60 and 100 pc for the south loop and the west rim, respectively, but they yield volume densities of ∼1−2 cm −3 and B sky strengths of ∼2−4 µG, implying unrealistic Alfvén velocities that exceed the expansion velocity of the superbubble. When we take the thickness of the cap as a lower limit to the depth, we obtain average volume densities of ∼20 cm −3 , which is rather high for such diffuse atomic clouds. The low H I emission brightness temperatures do not support large amounts of cold optically thick gas along the lines of sight. We therefore assumed mean depths of 40 ± 20 pc and 50 ± 20 pc for the south loop and the west rim. For the other clouds, we assumed that the transverse size applies in the third dimension, as suggested by the compactness of MBM 20 and G203-37 and the filamentary structure of Eridu. We assumed a mean weight µ = 1.36m H per hydrogen atom in the gas. The resulting gas number densities, n H , are listed in Table 2. The values range between 3 and 10 cm −3 in the shells, in agreement with superbubble simulations (Kim et al. 2017). Twice higher values are found in the north rim, as expected from the optically thicker conditions in H I and the dark neutral gas (see Paper II). The highest values correspond to the two compact molecular clouds, MBM 20 and G203-37.

B sky estimates
The derived values of B sky are listed in Table 2 and displayed in Fig. 11. Our assumption of a low dispersion in polarisation orientation angle is valid: it is below the 25 • limit proposed by Ostriker et al. (2001) from simulations, so that the ordered component of the magnetic field is much larger than the turbulent component. We obtain B sky strengths ranging from ∼3 to ∼15 µG towards the superbubble clouds. We can compare these values with the Sofue et al. (2019) results, which are based on synchrotron emission and Faraday rotation measures. They have derived all-sky maps of B los and B tot = (B 2 sky + B 2 los ) 1/2 , assuming equipartition between the magnetic and cosmic-ray energy densities. In the total magnetic field map, the Orion-Eridanus superbubble stands out with field strengths of 10 to 14 µG compared to a local average of 7 µG. For the two directions in our sample where B tot is available, we obtain values of 8.2 ± 1.2 and 12.3 ± 1.5 µG, in good agreement with the 12 µG value found at these locations in the Sofue et al. (2019) map. We further discuss magnetic compression along the outer shell of the superbubble in the next section. We find B sky strengths that are a factor of 17 to 18 times lower than the Soler18 estimates in the same directions along the west rim. The difference stems from both higher values of the line-of-sight depth and the Cho & Yoo (2016) correction. The effect of the modified DCF method alone would lower their magnetic strengths by a factor of 5 in these directions.
We have an estimate of the magnetic field in, or close to the Local Bubble wall, towards MBM 20. The B los value from Heiles (1989) is consistent with that of Xu & Han (2019). Our estimate of B sky of 2.3 ± 0.5 µG is lower than the = 8 +5 −3 µG derived by Andersson & Potter (2006) using the DCF method, but this difference may originate from the fact that we have probed the magnetic field in the atomic gas phase, whereas they performed their measurement in the molecular gas, which generally exhibits stronger fields (Crutcher 2012).

Outer shock velocity and compression ratio
We find that the fields along the outer shell of the superbubble (along the west rim and the south loop) span B sky values from 5 to about 15 µG, which means that the total field downstream of the expanding shock wave exceeds the mean total field B tot = 7 µG that prevails outside the superbubble, at high Galactic latitudes (Sofue et al. 2019). Xu & Han (2019) reported a similar enhancement by a factor of 2 in the shell of the Gum nebula, using rotation measures from pulsars. The B sky strengths towards the two north rim directions mostly probe the magnetic field in the front part of the outer shell. Their values of 5.3 ± 1.3 µG and 8.3 ± 2.1 µG are much higher than the ∼1 µG field strengths found in B los in the same directions (Sofue et al. 2019), therefore the field lines are highly inclined to the lines of sight and to the shock velocity. A similar configuration has been observed in the wall of the Local Bubble, where Xu & Han (2019) have found a value of B los = 0.5−2 µG much lower than the B sky = 8 +5 −3 µG derived by Andersson & Potter (2006) using the DCF method.
These low strength B los ( 5 µG) are also consistent with estimates from rotation measures of pulsars. The Australia telescope national facility (ATNF) pulsar catalogue lists 11 pulsars in the region of analysis, 2 of which are in the distance range of the Orion-Eridanus superbubble: J0452-1759 at 400 pc towards (α, δ) = (73.1 • ,−18 • ) and J0459-0210 at 160 pc towards (α, δ) = (75 • ,−2.2 • ). The ratio of the rotation measure to the dispersion measure gives the line-of-sight integrated B los , weighted by the thermal electron density along the line of sight. The B los strengths of 0.34 ± 0.01 and 1.05 ± 0.53 µG towards the two pulsars are consistent with the Zeeman measurements in the atomic phase, but the latter mostly samples the magnetic field in the Local Bubble because of the close distance to the J0459-0210 pulsar.
As indicated by the dust polarisation orientation that closely delineates the outer shell in Fig. 11 and by the large inclination of the field lines to the shock velocity in the front wall, the measured field strengths likely result from magnetic compression of the external field by the expanding shock wave. The atomic gas shells also result from compression and rapid cooling of the swept-up gas. For a shock velocity of ∼20 km s −1 , a compression ratio of 4, a mean gas density of 5 cm −3 , and the local interstellar cooling function (Richings et al. 2014), the radiative cooling length scale downstream of the shock is smaller than 0.1 pc. We can therefore consider the shock as isothermal to infer its compression ratio, r, and its velocity, v sh , from the downstream data.
The expansion velocity, v exp , of the downstream H I shells in the local standard of rest is related to the shock velocity as v exp = v sh (1 − 1 r ). Mass conservation and the isothermal condition relate the mass densities, ρ, gas velocities in the shock frame, u, and gas pressures, p, on both sides of the shock as r = ρ d /ρ u = p d /p u = u u /u d , where the subscripts u and d note the upstream and downstream media, respectively. Magnetic flux conservation relates the magnetic field strengths perpendicular to the shock velocity as B d⊥ /B u⊥ = r. Using these jump relations to eliminate the upstream variables that are not measured, we can write the momentum conservation equation as We assumed a temperature T u = T d = 8000 K typical of the warm atomic gas to relate the downstream pressure p d to the mean downstream gas densities listed in Table 2. We numerically solved Eq. (4) for the compression ratio r in each direction.
The magnetic field component perpendicular to the shock velocity is close to B sky towards the two north rim directions that probe the front wall, whereas both B sky and B los contribute to B d⊥ along the western and southern rims. We took the corresponding B sky and B los values from Table 2 Fig. 12. Estimates of the velocity and compression ratio of an 8000 K isothermal shock constrained by the downstream gas densities and magnetic field strengths listed in Table 2 and by the LSR expansion velocity of 15 km s −1 of the downstream H I shells. The different directions are noted with their truncated right ascension, and the letters N, W, or S for the north rim, west rim, and south loop rims, respectively. 0 B los 3 µG along the west rim and 3 B los 5 µG along the south loop. We note that the resulting shock velocities and compression ratios vary by less than 5 and 10%, respectively, when we include or exclude the Sofue et al. (2019) B los values in the calculations. The results obtained for v exp = 15 km s −1 are shown in Fig. 12.
The data point with little dispersion to an outer shock velocity close to 20 km s −1 and a compression ratio near 4. The downstream Alfvén velocities of 3-14 km s −1 are compatible with magnetic compression of external fields with Bu⊥ strengths ranging from 1 to 5 µG around the superbubble. The data require upstream pressures of 6000-11 000 K cm −3 along the western and southern rims that are compatible with the Local Bubble environment. Higher upstream pressures of 13 000 and 30 000 K are required to explain the higher gas densities seen towards the two north rim directions, but the downstream densities in these directions have larger uncertainties because the thickness of the front shell is unknown.

Cloud shells and outer shock expansion
We have studied the distribution in position and velocity of H I and CO emission lines towards the Eridanus part of the Orion-Eridanus superbubble. We decomposed the gas distribution into individual cloud complexes that appear as rings or shells of neutral gas that flank the H α recombination arcs. Arc B extends alongside the west rim, Arc C along a part of the south loop, Arc A is related to the molecular part of the north rim, and the east shell is parallel to an unnamed H α arc that is visible at negative velocities. Visual extinctions of 1-2 magnitudes along the northwest part of the outer rim hinder the detection of the H α recombination arcs.
The clouds also appear as separate entities in distance in 3D dust reddening maps. Their location, ∼150 to ∼250 pc away from the Sun, confirms that the superbubble is oriented with its far end towards Orion and its close end towards Eridanus. The data suggest that the south loop is the closest and approaching end of the superbubble. The central velocities of the H I lines trace the bulk motions of the gas. Their distribution highlights another coherent ring at about +9 km s −1 that likely marks the outer shell of the superbubble. We find expansion velocities of 10-15 km s −1 about this rest velocity. These motions are slower than the 40 km s −1 proposed by Brown et al. (1995) from the highest velocities recorded in the wings of the H I lines, but they are consistent with the values of 15-23 km s −1 found in H α (Reynolds & Ogden 1979).
The gas shells likely result from compression and rapid cooling of the gas swept up by the outer shock wave of the superbubble. The average gas densities, the short radiative cooling length scale, and the enhanced strength and orientation of the magnetic fields along the shells support this interpretation. The measured expansion velocity of 15 km s −1 therefore corresponds to the downstream gas velocity, which differs from the shock velocity itself. We have inferred a shock velocity of 18-23 km s −1 for an isothermal shock constrained by the downstream data on gas and magnetic fields and for a gas temperature of 8000 K. The compression ratio of 3-6 implies magnetic and gas pressures upstream of the shock that are consistent with the Local Bubble values.

Hot gas interior
We have used the ROSAT data in the three energy bands around 0.25, 0.75, and 1.5 keV and calculated the optical depth in these bands for the individual gas shells and for the total gas column densities measured along the lines of sight. We used this information to study the hot gas inside the superbubble and to locate the clouds relative to the hot gas. We also modelled the X-ray spectra in ten specific zones using mekal emission models and the gas absorption to infer the thermal properties of the hot gas.
We find that the east shell nicely bounds the emission at 0.75 and 1.5 keV. This boundary is not entirely due to absorption because the cloud is optically thin to 1.5 keV X-rays and partially thin (τ X of 0.1-0.4) at 0.75 keV. The cloud is optically thick at 0.25 keV, but it does not heavily absorb the softest X-rays, meaning that it lies inside the superbubble rather than in front of the hot gas. In contrast, the south loop is closer and absorbs the 0.25 keV emission. We used our gas data to revisit the origin of the soft X-ray emission towards the south loop. As suggested by Snowden et al. (1995), we find that significant 0.25 keV emission remains inside the loop after the foreground emission from the Local Bubble and the background intensity are subtracted. The background intensity is the intensity that is expected from the Galactic halo and is modulated by the total gas absorption.
Taking the south loop absorption into account, we find no clear spectral difference between the EXE1 region north of the east shell and the EXE2 region towards the south loop. The spectral dichotomy that is visible on the maps is largely induced by the south loop absorption. A single-temperature plasma cannot match the X-ray spectra across the three bands in any of the sampled directions, however. Two hot phases are required, with respective temperatures of 3-9 MK and 0.3-1.2 MK under the assumption of a uniform pressure across both phases. The derived pressure of (3−8) × 10 4 cm −3 K exceeds the pressure of the Local Bubble of 10 4 cm −3 K (Puspitarini et al. 2014). The need for at least two plasma phases might indicate global plasma oscillations across the superbubble that are caused by its complex interior history (Krause et al. 2014). In addition, the significant emission hardening seen towards the region that is enclosed between the H α arcs A and B corresponds to a higher pressure and hotter temperature that support the hypothesis of a younger age and faster expansion of this sub-region (Heiles et al. 1999).
The origin of the south loop ring is unclear. A wind-blown bubble seems unlikely because of the lack of massive stars in this direction (Burrows et al. 1993). We find no gradient in temperature or pressure between the south loop and adjacent regions that would support a cooling flow of hot gas that leaks at the rear of the superbubble (Heiles et al. 1999). The gas ring, filled with plasma that is a few million degrees hot, is consistent with an old supernova remnant. The progenitor star may have belonged to the Cassiopeia-Taurus OB association, which extends from 130 to 300 pc (de Zeeuw et al. 1999;Pellizza et al. 2005) in these directions. The more sensitive and refined X-ray maps that will soon be obtained by the extended Roentgen survey with an imaging telescope array (e-ROSITA) will spatially and spectrally constrain the heterogeneity of the hot gas phases that fill the superbubble better. This will enable us to reconstruct the possible sequence of supernovae that have blown this structure.

Magnetic field in the superbubble
We have used Planck polarised dust emission and the gas dynamics information from the H I line decomposition to study the magnetic field in the superbubble. We probed the plane-of-sky magnetic field, B sky , using the modified Davis-Chandrasekhar-Fermi method (DCF method) proposed by Cho & Yoo (2016). We have obtained magnetic strengths 17 times lower than previous estimates along the west rim (Soler et al. 2018). The difference stems from (i) the selection of regions with little H I confusion along the lines of sight, (ii) refined estimates of the gas shell depths along the lines of sight, and (iii) the correction of the dispersion in polarisation angles due to the pile-up of independent eddies along the line of sight.
The expanding superbubble has likely swept up and compressed the external magnetic field component perpendicular to the shock velocity. The polarisation alignment with the gas shell along the west rim and south loop, the enhanced B sky strengths that reach up to 15 µG, and the estimates of the Alfvén velocities in the downstream gas support this interpretation. The high ratio found between the B sky and B los strengths towards approaching gas structures in the north rim is also consistent with magnetic compression in the front wall of the superbubble. The DCF method only provides approximate field strengths, however. The 0.7 ξ 1 span in correction factor and the uncertainty on the line-of-sight depths lead to large uncertainties in the resulting B sky . Furthermore, large fractions of the superbubble could not be probed in B because of the limited S/N in polarisation. Future dust polarisation observations at higher sensitivity and angular resolution will enable a finer sampling of B sky and will allow additional corrections, for instance of the polarisation smoothing in the telescope beam and the effect of small-scale ordered density gradients across the aperture used to measure the magnetic dispersion Houde et al. 2009).  In order to analyse the X-ray emission from the superbubble, we have computed the X-ray optical depth maps using τ X (E X ) = σ(E X , N H ) × N H , (A.1) where E X is the X-ray band, σ the energy-band-averaged photoelectric absorption cross section from Snowden et al. (1994) and N H the gas column densities inferred in the different gas phases from our combined H I, CO, dust, and γ-ray study. The optical depth maps obtained for the total gas column densities τ tot X and for the three ROSAT energy bands are displayed in Fig. A.1.

Appendix B: Angular dispersion of the magnetic field
To derive the angular dispersion of the magnetic field, we have first used the direct method described in Eqs. (D.5) and (D.11) of Planck Collaboration Int. XXXV (2016), and where ψ(x) denotes the polarisation angle at the position x in the sky, ... is the average over the circle with a diameter of 2 • , arctan(sin, cos) is the arc-tangent function that solves the π ambiguity taking into account the sign of the cosine, and Q and U are the Stokes parameters. The second method relies on the structure function of polarisation angles, S 2 (δ), defined in Planck Collaboration Int. XIX where ∆ψ xi is the angle difference between the polarisation at the position in the sky x (central pixel) and the polarisation at a position displaced by the vector δ i . The sum is over an annulus around the central pixel x of radius δ = |δ| (the lag), width ∆δ and containing N pixels. We chose a width ∆δ = δ as advised in Planck Collaboration Int. XIX (2015). The result was then averaged over all the positions x in the circle of 2 • diameter, leaving only a dependence on δ. Hildebrand et al. (2009) considered a magnetic field composed of a field with a large-scale structure, B 0 , and a turbulent component, B t . Assuming that both contributions are statistically independent, they made the approximation S 2 2 (δ) = s 2 + mδ 2 , the first term being the contribution of B t and the second term the linear contribution of the smoothly varying B 0 . We derived the intercept s by fitting the structure function with this model on the range δ = [1 • , 2 • ], above the 40' resolution of the data. s was then linked to the plane-of-sky magnetic field, taking into account the Cho & Yoo (2016)  We consider here small regions compared to the Planck beam size, however, which implies few independent vectors in our apertures with 2 • diameter. As a consequence, the structure function method that probes the magnetic field at different angular scales gives very similar results to the direct method that only probes the average, as shown in Fig. A.2. Only the magnetic field strengths from the direct method are displayed in Table 2.