© The Authors 2023

1. Introduction

The response of interstellar gas to energy and momentum deposition from supernovae and stellar winds is the growth of a hot bubble surrounded by a dense shell. Star-forming populations can inject energy and momentum long enough to sweep the shell up and out to lower density gas, where Rayleigh-Taylor instabilities eventually deform it before breaking it up. The hot gas then vents into the circum-galactic halo of the galaxy as a wind that also entrains warm and cold gas. If the wind is powerful enough, the gas may escape the gravitational potential of a galaxy and thereby enrich the intergalactic medium (IGM) with metals. The gas from less powerful outflows then rains back on the galaxy and is available to form the next generation of stars. Outflows and winds – next to inflows of fresh gas – are a cornerstone of galaxy formation models and are vital for regulating cosmic chemical evolution (see reviews by Veilleux et al. 2005, 2020; Collins & Read 2022; as well as Schneider & Robertson 2018; Nelson et al. 2019; Mitchell et al. 2020; Schneider et al. 2020; Pandya et al. 2021 for state-of-the-art computer models).

Observationally, signatures of galactic outflows in star-forming galaxies are ubiquitous, both in the nearby (e.g. Heckman et al. 2011; Chisholm et al. 2016) and in the high-redshift universe (review by Erb 2015). Galaxies in the local Universe allow for detailed panchromatic mappings of outflow phenomena, from the highest energies (cosmic rays and hot gas) to the longest wavelengths possible (cold and molecular gas) – especially if the outflow occurs ‘edge-on’. The most detailed observational studies that analyse multiple phases of outflows simultaneously are performed on nearby, more evolved, and massive systems. These observations nourish our understanding and provide robust constraints of how mass, energy, and momentum from star formation couple with the different phases of the interstellar-medium (ISM) to drive outflows¹. However, from a cosmological perspective, similar constraints are especially needed for the abundant and, in the early universe, dominating population of low-mass star-forming galaxies. Currently, cosmological models have to implement those processes in the form of phenomenological ‘sub-grid physics’ recipes that are tuned to match observed galaxy population statistics, especially luminosity- and stellar-mass functions (Somerville & Davé 2015).

Understanding the effects of feedback and winds in low-mass systems is especially required to grasp the large-scale physics of the universe during the cosmological phase transition known as the Epoch of Reionisation at z ≳ 6. This is because dense shells in the pre-fragmentation stage are more opaque to hydrogen ionising radiation than perforated bubbles at a later stage (e.g. Fujita et al. 2003). Theoretical modelling requires sustained momentum and energy injection to drive a wind and thereby significant Lyman continuum (LyC) photons out of low-mass galaxies (Kimm & Cen 2014; Paardekooper et al. 2015). Detailed ‘zoom-in’ simulations of individual dwarf star-forming galaxies suggest that LyC escape is highly variable and highly anisotropic (Trebitsch et al. 2017). The angular dependence is modulated by high optical depth sight-lines along cold-flow accretion filaments (Park et al. 2021) and low optical depth sight-lines that follow the direction of the outflow phenomena (Trebitsch et al. 2017). Such anisotropic beamed LyC may have consequences for the mechanisms by which H II zones percolate during reionisation (Paardekooper et al. 2015; Furlanetto & Oh 2016).

While a pan-chromatic approach is required for a detailed understanding of the multiple phases of the outflowing material, narrow band imaging and kinematic analysis of the classical principal emission lines can expose and characterise feedback, outflow, and wind phenomena just from the T ∼ 10⁴ K phase (e.g. Calzetti et al. 2004; Westmoquette et al. 2008; Zastrow et al. 2011; Lee et al. 2016; McQuinn et al. 2019). The increased sensitivity for low surface-brightness line emission offered by modern integral-field spectrographs (IFSs; Bacon & Monnet 2017) allows for unparalleled discoveries in this respect.

Observational evidence for the T ∼ 10⁴ K phase of feedback-driven winds in nearby dwarf galaxies and potential early-universe analogues have been amassed for decades from such imaging and spectral analyses (e.g. Marlowe et al. 1995; McQuinn et al. 2015, 2018; Collins & Read 2022). Combined H I–H II morpho-kinematical investigations of those targets are feasible as well. The studies by van Eymeren et al. (2009a,b, 2010) are noteworthy as they compare 21 cm HI and HII kinematics from Hα in five nearby irregular dwarf galaxies (NGC 4861, NGC 2366, NGC 5408, and IC 4622). In these galaxies, multiple localised velocity offsets between the ionised and the neutral phase are found and are interpreted as being driven by feedback (see also Bomans 2001; Bomans et al. 2007 for reviews and X-ray follow-up). This idea is corroborated by filamentary features seen in Hα that extend towards regions of lower neutral columns, as well as the detection of partially broken super-shells in Hα. However, in the van Eymeren et al. (2009a) studies, feedback was never found to be strong enough to drive outflows into the IGM. A similar conclusion was reached by Westmoquette et al. (2008) for the dwarf starburst NGC 1569, where the hot wind fluid is still presumed to be confined by super shells that extend up to 1 kpc outwards from the super-star clusters in this system. Roychowdhury et al. (2012) compared H I with H II in three compact dwarf galaxies but found only tentative evidence for outflows. More recently, McQuinn et al. (2019) analysed Hα and 21 cm morphologies simultaneously in a sample of 12 nearby dwarfs. They traced potential galactic winds via low surface-brightness Hα emission isophotes (≳3 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻²) that extend beyond regions of their limiting neutral column (N_HI ∼ 10²⁰ cm⁻² or 1 M_⊙ pc⁻²) gas; however, even the most convincing winds found in this study were not deemed strong enough to have material escaping into the IGM.

More compelling early universe analogues for such observational studies are small, low-mass, and low-metallicity galaxies that exhibit high specific star-formation rates – so-called blue compact dwarf galaxies (Gil de Paz et al. 2003, see also the review by Thuan 2008 and the recent, more general review on dwarf galaxies by Henkel et al. 2022). Their shallow gravitational potential and the spatially and temporally concentrated energy and momentum injection from the starbursting population make them prime candidates for powering large-scale outflows (Mac Low & Ferrara 1999). And indeed, the unprecedented sensitivity for mapping low-surface brightness emission lines offered by the Multi Unit Spectroscopic Explorer (MUSE) at ESOs Very Large Telescope (VLT) UT4 (Bacon et al. 2014) unveiled spectacular outflows around the blue compact dwarfs Haro 11 (Menacho et al. 2019), Haro 14 (Cairós et al. 2022), and ESO 0338−IG04 (Bik et al. 2015, 2018), as well as for the nearby HII galaxy Heinze 2−10 (Cresci et al. 2017).

In this paper we report new observational results on H I–H II morphology and kinematics from the interstellar and circumgalactic medium of the observationally well-studied blue compact dwarf galaxy SBS 0335−052E (e.g. Izotov et al. 1990, 2001; Thuan et al. 1997; Papaderos et al. 1998; Johnson et al. 2009; Adamo et al. 2010; Hunt et al. 2014; Kehrig et al. 2018; Wofford et al. 2021, and references therein). Numerous properties render this galaxy a unique analogue of early universe galaxies. With an absolute UV magnitude of M_UV = −16.8, it is significantly fainter than the characteristic absolute UV magnitude of the z ≳ 6 luminosity functions ($M_{\mathrm{UV}}^{\star }<-20.6$; Bouwens et al. 2015). With 12 + log(O/H)∼7.25 (Izotov et al. 2009), it is one of the most metal-poor galaxies in the local Universe. This low-mass galaxy (M_⋆ = 6 × 10⁶ M_⊙; Reines et al. 2008) is also extremely compact; the starburst is taking place within 500 pc, where six super-star clusters with diameters ≲60 pc form stars at ∼0.7 M_⊙ yr⁻¹ or at specific star-formation rates up to 20 M_⊙ yr⁻¹ kpc⁻².

Our previous analysis of VLT/MUSE integral-field spectroscopic data from this galaxy unveiled two prominent ionised filaments in its halo (Herenz et al. 2017, hereafter H17). Here we present an analysis of new MUSE data in concert with new Karl G. Jansky Very Large Array (VLA) B-configuration observations that show that these filaments are likely related to an unprecedented large conical outflow structure that pinches deep into the halo of this compact starburst galaxy (Sect. 4). Prior to this, we summarise observations and data reduction in Sect. 2; also, the morphological and kinematic analyses of neutral and ionised gas are detailed in Sect. 3. In Sect. 5 we summarise and conclude. For distance conversions we assume a flat Λ cold dark matter Universe with Ω_m = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹, and we adopt z = 0.01352 as the cosmological redshift for SBS 0335−052E (Moiseev et al. 2010). This translates to a luminosity distance of 58.42 Mpc and an angular scale of 276 pc/″.

2. Observations and data reduction

2.1. VLT/MUSE

For the present analysis we combine data from two ESO/VLT MUSE programmes: Programme 096.B-0690 (PI: Hayes) with data taken on November 16th and 17th in 2015 (clear to photometric sky, DIMM seeing ≈0.7″–0.9″) and Programme 0104.B-0834 (PI: Herenz) with data taken on October 2nd and 3rd in 2020 (photometric sky, DIMM seeing 1.1″–1.4″). Both programmes were observed in wide-field mode (1′×1′ field of view) with the blue cut-off filter removed (extended mode, wavelength range 465 nm–930 nm) and without adaptive optics. Programme 096.B-0690 centred the field of view on SBS 0335−052E (RA: 03^h37^m44.075^s Dec: −05°02′39.5″). These data lead to the discovery of ionised filaments in the halo of SBS 0335−052E (H17), and were also analysed by Kehrig et al. (2018) and Wofford et al. (2021) to understand the properties of the ionising sources of this galaxy. In order to map the extent of the filament we obtained new MUSE data in programme 0104.B-0834 with a pointing centred north-west of the galaxy (RA: 03^h37^m41.9167^s Dec: −05°01′57.2″). This pointing overlaps partly with the discovery data. The total open shutter times are 5680 s (8 exposures a 710 s) and 5576 s (8 exposures a 697 s) for the central- and offset pointing, respectively.

We reduced the data with version 2.8.3 of the MUSE data processing pipeline (Weilbacher et al. 2020). The reduction of the 2015 data was already detailed in Wofford et al. (2021, Sect. 3.2) and we apply here the same recipes to our new 2020 data. These recipes provide us first with sky-subtracted, flux calibrated, and astrometrically roughly calibrated pixel tables for each exposure. As explained in H17, the spatially extended line emission in the data requires a modification of the default pipeline sky subtraction procedure; significant flux from the extended principal H II emission around SBS 0335−052E would be removed without this modification. The most prominent oversubtracted emission shows up as spikes at the observed wavelengths of Hα and [O III] λ5007 in the sky continuum spectra that the pipeline determines after the best-fit sky emission line spectra have been subtracted (Weilbacher et al. 2020, Sect. 3.9.1). We thus altered the output sky continuum spectra from the initial run of the sky subtraction routine by replacing the values around principal emission lines via linear interpolation from the surrounding spectral bins. The corrected sky continuum spectra and the unaltered sky emission line spectra from the first run were then used as input in a second run of the sky subtraction procedure. This removed the visible oversubtraction artefacts at Hα and [O III] λ5007 in the initial reduction, and to be secure we also repeated the procedure for the lines [O III] λ4363, He IIλ4686, He Iλ5876, [S II] λλ6717,6731, and [O I] λ6300, even if such strong spikes were apparent in the continuum sky spectrum from the initial run of the sky subtraction. Lastly, the individual exposure pixel tables are resampled onto a 3-dimensional 473 × 532 × 3802 grid (datacube); the first two axes sample the spatial domain with 0.2″ × 0.2″ spectral pixels (spaxels) parallel to RA and Dec and the third axis samples wavelength linearly from 4599.96 Å to 9351.21 Å with 1.25 Å bins. The final coordinate transformation between datacube spaxels and (RA, Dec)-coordinates is registered by cross-correlating the Pan-STARRS (Chambers et al. 2016; Magnier et al. 2020) r-band image against a synthetic r-band image created from the datacube.

2.2. VLA

We obtained new HI 21 cm observations with the Karl G. Jansky Very Large Array (VLA) in its B-configuration. The observations (project number 17B-234; PI: Herenz) made use of the L-Band receiver centred on the redshifted 21 cm frequency of the galaxy. The programme was executed in 10 × 3h blocks from Sep. 30th to Oct. 20th, 2017.

All raw VLA data are calibrated following standard prescriptions using the National Radio Astronomy Observatory (NRAO) casa software package (version 5.6.1-8; McMullin et al. 2007; Emonts et al. 2020). Radio-frequency interference is manually identified and excised from the raw data. These data are then calibrated, and once a satisfactory calibration is achieved, the data are continuum subtracted and weighted using CASA’s uvcontsub and statwt functions. “Satisfactory” calibration constitutes an rms noise value in source-free channels of the cube that is within 10% of the theoretically optimal rms value given the duration of the observation, the configuration of the VLA, and the spectral properties of the datacube.

Using tclean and the masking algorithm Auto-Multithresh (AMT; Kepley et al. 2020), we generate an intermediate image cube using casa version 6.3.0.48. From this intermediate cube, we extract the AMT-generated mask. We chose AMT for mask generation as it produced high quality data products with reproducibility and speed. All image cubes are created using Briggs weighting with a robust parameter of 0.5 for a compromise between angular resolution and sensitivity.

We set the AMT noise threshold to 3σ. In doing so, AMT finds and masks all emission associated with our target sources, in addition to noise peaks that rise above this threshold. To address this, we then examine the AMT-generated mask interactively and extract only the regions that include source emission. Through this process of automated mask generation and interactive selection, we end up with a final, robustly generated mask that solely and fully masks the target sources. With this mask in hand, we then initiate a final tclean with the same parameters as used for the mask-generation clean, but now providing the source-only mask generated in the previous step. We clean down to 0.5σ, as this removes the need for residual flux rescaling. The resulting final image cube covers an area of 22.8 arcmin² sampled at 1″/pixel, and the channel width is 10.26 km s⁻¹. The typical rms-noise per channel is 2.5 × 10⁻⁴ Jy beam⁻¹. The beam is 6.1″ × 4.6″ and it is oriented 23.6° west of north. The positional accuracy of the VLA HI data, which is based on milli-arcsecond precision positions of quasars in the VLA calibrator manual², can be assumed to be typically 10% of the synthesised beam size.

We also create another datacube from the B-configuration data, that provides us with a different spatial resolution. This cube is created by applying a 7 kλ UV-taper in the imaging process. This UV-tapering downweights the longest baselines using a Gaussian function in the uv-plane (with the FWHM measured in units of wavelengths), producing a cube at a spatial resolution that is roughly equivalent to VLA C-configuration observations (beam 15.9″ × 14.7″ at 66.4° east of north). The datacube differs from the untapered cube by an adjusted spatial sampling of 3″/pixel and the rms-noise per channel is 4 × 10⁻⁴ Jy beam⁻¹. The 7 kλ-tapered data are less sensitive than the untapered data, but this is compensated by the larger beam that enables the mapping of lower columns over larger areas.

3. Analysis and results

3.1. Detection and morphology of extended H II and H I emission

3.1.1. Ionised gas

To optimally detect extended line emission in our MUSE data, we processed the reduced datacube with routines from the LSDCat software (Herenz & Wisotzki 2017). We started by subtracting a running median in the spectral direction. The band-width of this filter (which effectively enables the removal of low-frequency continua) is 151 px (188.75 Å). We then cross-correlated the continuum-subtracted datacube with a 3D Gaussian template. The result of this procedure is a unit-less datacube, that we dub in LSDCat as “S/N cube” in reference to the signal-to-noise ratio maximising characteristic of matched filtering. However, we here use LSDCat only as a simple linear filter to suppresses high-frequency noise and enhance the emission line signal in the datacube. By experimentation we found that the parameters FWHM_v = 150 km s⁻¹ and σ_G = 2″ optimally augmented the spectrally narrow and spatially extended emission line features in the circumgalactic medium around SBS 0335−052E. We visually examined the resulting S/N cube with the QFitsView software (Ott 2012) in layers with strong nebular line detections in the integrated spectrum (cf. Fig. 3 in Wofford et al. 2021). We here report significant detections of extended filamentary emission in Hα, [O III] λ5007, and Hβ around SBS 0335−052E (Fig. 1). Hβ was not reported in H17, but the increased exposure time in the region where the discovery data overlaps with our new observations renders the detection of this line significant. The filaments remain undetected in other emission lines.

Fig. 1. H II and H I emission as seen with MUSE and VLA around SBS 0335−052E. Top panel: continuum-subtracted Hα narrowband image. Fluxes are encoded via an asinh-scaled cyclic colour map from 0 to 10−18 erg s−1 cm−2 in the low surface-brightness, and from 10−18 erg s−1 cm−2 to 10−16 erg s−1 cm−2 in the central high surface-brightness region. The cyan solid contours demarcate Hα isophotes at SBHα = {1.5, 2.5, 5, 12.5}×10−18 erg s−1 cm−2 arcsec−2 while the white dotted contour demarcates the limiting surface-brightness of SBHα = 7.5 × 10−19 erg s−1 cm−2 arcsec−2. The image has been smoothed with a Gaussian of 0.4″ FWHM. Subdued grey contours indicate HI column densities NHI = {5, 10, 20, 30, 40}×1020 cm−2 from our VLA-B configuration observations; these contours are displayed more prominently in Fig. 4. The dashed white contour demarcates the 2.5σ = 2.5 × 1020 cm−2 detection limit. The VLA-B configuration beam, 6.1″ × 4.6″ oriented 23.6° west of north, is indicated via a grey shaded ellipse in the top left. The displayed field of view is 1′35″ × 1′46″ corresponding to 26.2 kpc × 29.2 kpc in projection. The dashed violet square in the centre indicates the viewport used in Fig. 4. North is up and east is to the left. Bottom panels: signal-to-noise of the three principal emission lines (Hα, Hβ, and [O III] λ5007 in the left, centre, and right panel, respectively) in which extended ionised gas is detected (see text for details). The maps are displayed with a linear cyclic colour map from 0 to 90 and from 90 to 900. White dotted contours mark out regions with SN > 8. The chartreuse dashed rectangles in the Hα panel outline the rectangular regions used for the line ratio analysis in Sect. 3.2.

The Hα narrowband image shown in Fig. 1 was synthesised by summing over the three layers in the median-filter subtracted datacube (6651.21 ± 1.25 Å observed frame). Only these layers contain appreciable Hα emission (see also the line-width analysis in Sect. 3.3.1). The maximum S/N-images in Fig. 1 were created by taking the maximum values from the LSDCat cross-correlated datacube over a few layers around the emission lines of interest. We define ${\mathrm{SN}}_{\mathrm{H}\alpha }^{\mathrm{LSD}}>8$ as a conservative criterion for detected emission line signal. Spaxels in the S/N cube that fulfil this criterion are surrounded by a white dotted contour in Fig. 1.

As can be appreciated from Fig. 1, the new MUSE observations reveal a continuation of the filamentary structure first reported by H17. This Hα emitting complex extends up to ≳15 kpc in projection towards the north-west from the galaxy’s main stellar body. While the pronounced high-SN emission stems from the two filamentary threads, Hα is also detected throughout a contiguous region in between the filaments. The [O III] and Hβ detections, on the other hand, are confined exclusively to the filaments. The northern filament shows quite a strong detection in [O III] (SN ≳ 25), whereas the north-western filament is significantly weaker in [O III] (SN ≳ 10). An analysis of the observed line ratios from the filaments is presented in Sect. 3.2.

To make the discussion of morphological features of the extended ionised structure text more accessible, we display in Fig. 2 Hα surface brightness contours with annotations. The two filaments that protrude towards the north (northern filament; PA ≈ 0°) and north-west (western filament; PA ≈ 315°) are characterised by SB_Hα ≈ 3 × 10⁻¹⁸ erg s⁻¹ cm⁻². The northern filament fragments after ∼19″ (∼5 kpc) into a lacier prolongation. This threadlike extension continues outwards up to ∼58″ (∼16 kpc). The western filament, on the other hand, brightens noticeably after ∼25″ (∼7 kpc) to SB_Hα ≈ 5 × 10⁻¹⁸ erg s⁻¹ cm⁻² (Hα “hot-spot”). After this brightening the western filament fans out into two short branches, one to the west and one to the north-east. The accented filamentary features are surrounded by diffuse emission at SB_Hα ≈ 7 × 10⁻¹⁹ erg s⁻¹ cm⁻². Two shorter and fainter filaments are found towards the north-east of the galaxy; these are labelled as “eastern foreshortened filaments” in Fig. 2. They may trace a similar structure as the northern and the western filaments, but appear shorter due to a different alignment along our sight-line.

Fig. 2. Annotated Hα surface brightness contours (contour levels SBHα = {1.5, 2.5, 5, 12.5, 100, 700}×10−18 erg s−1 cm−2 arcsec−2). In the background the Hα narrowband is shown with a log-stretch to 5 × 10−16 erg s−1 cm−2. The thin-dotted line surrounding the filaments to the north-west demarcates the aperture that is used to estimate the total Hα flux from that structure (Sect. 4).

While Hα in the northern circumgalactic halo shows complex morphological features, the southern halo displays a rather smooth and slowly decaying light profile. At our sensitivity we detect diffuse Hα emission up to 25″ away from the main stellar body of the system. To provide a quantitative assessment of this emission we extract a radial surface-brightness profile, SB_Hα(r), under the assumption of circular symmetry. Therefore, we first determine the photometric centre $(\overline{x},\overline{\mathit{y}})=(\sum _{x\mathit{y}}{I}_{x\mathit{y}}x/\sum _{x\mathit{y}}{I}_{x\mathit{y}},\sum _{x\mathit{y}}{I}_{x\mathit{y}}\mathit{y}/\sum _{x\mathit{y}}{I}_{x\mathit{y}})$, where I_xy is the Hα flux value of the pixel at coordinates (x, y). We find $(\overline{x},\overline{\mathit{y}})=(153.7,162.7)$ which is located at (RA, Dec) = (3^h37^m43.99^s, −5° 02′39.2″). We then extract SB_Hα(r) in semicircular annuli of 0.6″ width south of $\overline{\mathit{y}}$ and centred on $(\overline{x},\overline{\mathit{y}})$. The result of this procedure is shown in Fig. 3. The measured Hα SB profile of SBS 0335−052E is characterised by a central high-SB region (“core”) which drops rapidly from 2 × 10⁻¹⁴ to 2 × 10⁻¹⁶ erg s⁻¹ cm⁻² arcsec⁻² within the first 5″ and then kinks into a more slowly decaying (“halo”) profile. Further out, at ∼13″ and ∼10⁻¹⁷ erg s⁻¹ cm⁻² arcsec⁻², another kink towards an even flatter profile is found.

Fig. 3. Hα surface brightness profile in the south, i.e. not affected by the filaments in the north. The points show the observed profile in circular annuli and error-bars are 2σ. The black line indicates the double exponential profile fitted to the outer (r > 5″) profile, while dot-dashed and dashed lines indicate the individual components with scale-lengths r1 = 1.62″ (447 pc) and r2 = 8.94″ (2.47 kpc). The grey line indicates the light profile of the central (r ≤ 5″) residual component that is obtained after subtracting the double exponential.

As can be seen from Fig. 3, the outer profile can be described adequately by a sum of two exponentials with scale lengths r₁ = 1.62″ (447 pc) and r₂ = 8.94″ (2.47 kpc). We caution that given the extreme intensity contrast – four orders of magnitude between core and outer halo region – diffuse scattered light may have brightened the observed profile at the faintest isophotes. This observational effect of scattered excess light was studied extensively by Sandin (2014, 2015). These papers caution against the physical interpretation of faint halo scale lengths around bright compact galaxies in the absence of an accurate extended point-spread function model for the observations. In particular the here measured flattening at larger radii was put forward as a tale-tale sign of diffuse scattered light (see especially Sect. 5 and Appendix E in Sandin 2015).

The “core” high-SB region (SB_Hα ≳ 2 × 10⁻¹⁶ erg s⁻¹ cm⁻²), denoted as “central circumgalactic zone” in Fig. 2, is characterised by an intricate morphology. This region is detected in a plethora of other emission lines in the MUSE IFS data (see figures in Kehrig et al. 2018). To visualise the morphology of this region adequately we display in the left panel of Fig. 4 a zoomed-in view of the Hα narrow band image from Fig. 1; we now use a different cyclic asinh-stretch that is appropriate to trace the brighter morphological features relevant at this magnification. To better visualise the size difference between circum-galactic ionised gas structures and main stellar body of the galaxy, we inset the Hubble Space Telescope (HST) F550M image into the brightest central (SB_Hα > 10⁻¹⁵ erg s⁻¹ cm⁻²) region. This band is devoid of strong emission lines and thereby traces only the morpholgy of the stellar continuum. We show the HST FR656N image in the right panel of Fig. 4. This resolves morphological features in the central region in greater detail and nicely demonstrates the gain in depth of our MUSE data with respect to HST.

Fig. 4. Left panel: zoomed view (22″ × 22″/6 kpc × 6 kpc) of the Hα narrow band image from the top panel of Fig. 1 with NHI contours from the VLA B-configuration 21 cm observations (NHI = {5, 10, 20, 30, 40}×1020 cm−2, rose lines); these contours are also shown in Fig. 1 as subdued white lines. The cyclic colour map encodes Hα flux in a asinh-scale from 0 to 10−17 erg s−1 cm−2 to 10−15 erg s−1 cm−2. In the brightest central region the archival HST F550M image is inset, with an arbitrary asinh-scaling to highlight the compact super-star clusters that comprise the main stellar body of this system. The violet dashed square indicates the region that is displayed in the right panel. Right panel: inner high-SBHα region as seen in the archival HST FR656N image (13″ × 13″/3.5 kpc × 3.5 kpc) displayed here with a cyclic asinh-scale from 0 to 5 × 10−15 erg s−1 cm−2 to 1013 erg s−1 cm−2. Here the violet dashed square indicates the region that is displayed in Fig. 5.

Figure 4 reveals that the high-SB Hα region exhibits multiple nested arcs or loops, as well as outward pointing filaments or spurs. Generally loops seen in Hα at such large scales trace the limb-brightened edges of so-called super-bubbles or super-giant shells (e.g. Chu 1995; Bomans 2001; Bomans et al. 2007). Disconnected loops are believed to point towards regions where the shells are broken up and where the pressurised interior can vent out. The interaction of the hot wind (> 10⁵ − 10⁷ K) with the colder material in the surroundings and warm gas (∼10⁴ K) entrained in these bubble-blow outs, as well as the merging of super-bubbles, are believed to manifest in outward pointing filamentary spurs (e.g. Cooper et al. 2008; Tanner et al. 2016). Here we see that the two brightest inner arcs towards the north-west resemble closed loops, but variations in surface brightness point at inhomogeneities that indicate cracks in the shells caused by Rayleigh–Taylor instabilities. Further outwards towards the north-west, only barely detected in the FR656N image but clearly visible in the MUSE data, we see a flaring of spur-like features. These structures are very extended, up to ∼10″ or 2.7 kpc from the main star-forming regions. These features may be interpreted as outflowing warm gas that is entrained in the hot wind that vents through the cracks in the shell.

3.1.2. Neutral gas

For our mapping of the signals from the neutral phase we processed the VLA datacubes with the software SoFiA (Source Finding Application; Serra et al. 2015). We used SoFiA’s smooth+clip source detection algorithm that is most commonly employed for searches of extragalactic HI emission. Here the radio datacube is smoothed with a set of top-hat kernels of varying sizes and the final detected source mask is built from the union of the binary masks that result from thresholding each individually smoothed cube. We used the recommended parameters for detection of extragalactic HI signal, i.e. a set of 12 top-hat kernels of zero, three, seven, and 15 channels width and spatial dimensions of 0 × 0 (i.e. no spatial smoothing), 3 × 3, and 6 × 6 pixels (1 pixel = 1 arcsec²), and a detection threshold of 5. The SoFia runs provide us with moment-0 maps, S^ν [Jy km s⁻¹ beam⁻¹], that are of interest here and moment-1 maps that are used for the kinematic analysis in Sect. 3.3.2. The moment-0 maps are converted into maps of neutral column, N_HI [cm⁻²], via

$$\begin{array}{c}\hfill N_{\mathrm{H}}=1.1\times {10}^{24}{(1+z)}^2{S}^{\nu }/(ab),\end{array}$$(1)

where a [″] and b [″] are the FWHM of the major- and minor axis of the synthesised beam, respectively, and z is the redshift of the galaxy (Eq. (78) in Meyer et al. 2017). The N_HI map of SBS 0335−052E from this procedure is displayed in Fig. 6 and the top panel of Fig. 7 shows the moment-0 map encompassing both galaxies SBS 0335−052E & W. The individual layers that contribute to the spatially resolved H I signal are assembled in Appendix B in the form of channel maps.

Translating the detection threshold of the smooth-and-clip algorithm into an N_HI sensitivity limit needs to account for the different smoothing kernels that contributed in building up the 3D source mask. We proceeded empirically by inspecting the final binary-mask cube from SoFiA to determine that the faintest pixels in the outskirts of the moment-0 map are mostly single channel detections. As stated in Sect. 2.2, and verified by the noise cube produced by SoFia, the typical rms-noise per channel in the B-configuration data is 2.5 × 10⁻⁴ Jy beam⁻¹. For a 10 km s⁻¹ wide-channel this translates via Eq. (1) into a 1-σ column-density limit of σ_NHI = 1 × 10²⁰ cm⁻². We finally used the contouring algorithm of ds9 (Joye & Mandel 2003) without smoothing on the N_HI map (Fig. 6) to produce contours at N_HI = {2.5 (2.5σ),5, 10, 20, 30, 40}×10²⁰ cm⁻². The closing 2.5σ contour associated with SBS 0335−052E is shown as a white-dashed line only in Fig. 1; all other contours are displayed in Figs. 1, 4, and 5.

Fig. 5. HST archival F550M image with super-star clusters labelled according to Thompson et al. (2009); NHI contours as in Fig. 4. The viewport, indicated as a violet dashed square in Fig. 4, is 7″ × 7″ or 1.88 kpc × 1.88 kpc.

Fig. 6. 21 cm moment-0 map of SBS 0335−052E with the colour-bar encoding the neutral column NHI in units of 1020 cm−2 using an asinh-scaling. The magenta ellipse indicates the dimension and orientation of the VLA B-configuration beam. This map has been used to draw the NHI contours in Figs. 1, 4, and 5.

From Fig. 4 it becomes apparent that the peak of N_HI is slightly offset to the west from the Hα peak and from the star-clusters. The offset with respect to the star-clusters is made more apparent in Fig. 5, where we overlay the N_HI contours onto the emission-line free HST/ACS F550M image. Such offsets are observed frequently in dwarf starbursts and it is hypothesised that that feedback from starbursting SSCs will mechanically disrupt and ionise H I in their vicinity (e.g. van Eymeren et al. 2010; Cannon et al. 2016; Teich et al. 2016; Jaiswal & Omar 2020).

At lower column densities the resolved H I structures appear to relate spatially with Hα morphological features. Towards the north we find an extended tail at N_HI = 10²¹ cm⁻² that bends towards the east at a slightly lower columns (Figs. 1 and 4). The northern tail is almost co-spatial with a northern Hα spur before it curves around the Hα spur towards the east. At even lower neutral columns, close to the detection limit of our observations (2.5σ = 2.5 × 10²⁰ cm⁻²; white dashed contour in Fig. 1), we trace an extended tail that aligns with the foreshortened filaments towards the east. Towards the north-west and more prominently towards the south-east, two “ears” emerge at the lowest column density, both of which appear also co-spatial with spurs pointing in the same direction.

Given the observed spatial correlations between HI and Hα we suspect that the low-column HI emission stems from cold gas that is interacting with the ionised phase. The entangled interactions between star-forming H II regions and H I phase are observationally well studied on small (≲10² pc) scales in more nearby systems (e.g. Egorov et al. 2014; Cannon et al. 2016; Egorov et al. 2018; and review by Veilleux et al. 2020). These processes are complex, observationally challenging to tackle, and far from being quantitatively fully understood. Here, in the extreme environment of the starburst, we may now have resolved such interactions on kiloparsec scales.

To study the interaction between the ionised and neutral phase on larger scales, we use the VLA B-configuration data tapered at 7 kλ (Sect. 2.2). The N_HI moment-0 maps and the 1-σ detection threshold, 2.2 × 10¹⁹ cm⁻², were computed in the same way as described above for the untapered dataset. The tapered data reveal H I within and around the pair SBS 0335−052E & W. We juxtapose the moment-0 maps for the two different beams in Fig. 7. We compare the large-scale HI contours from tapered data with our MUSE Hα image in Fig. 8. There we use Pan-STARRS imaging (Chambers et al. 2016; Magnier et al. 2020) outside the MUSE FoV. In Pan-STARRS the fainter companion SBS 0335−052W is prominently visible as a blue cometary object close to the position of the western N_HI peak.

Fig. 7. Moment-0 maps from VLA B-configutration 21 cm observations. The colour mapping encodes NHI in 1020 cm−2 and the beam for each dataset is indicated as a black ellipse. Top panel: VLA-B configuration data. Bottom panel: VLA-B configuration observations, tapered at 7 kλ; this map has been used to draw the NHI contours in Fig. 8.

The juxtaposition of the moment-0 maps in Fig. 7 illustrates how different beams are sensitive to different scales of the neutral gas distribution in and around the galaxies. The untapered data resolves the morphology of higher-column H I within the galaxies and in their immediate vicinity. The measured peak N_HI = 4.1 × 10²¹ cm⁻² for SBS 0335−052E is only slightly lower than the inferred $N_{\mathrm{H}}^{\mathrm{COS}}$ from Lyα absorption with HST/COS spectroscopy: $N_{\mathrm{HI}}^{\mathrm{COS}}=(5.0\pm 0.5)\times {10}^{21}$ cm⁻² (James et al. 2014). However, this does not indicate that the gas seen by COS in absorption is spatially resolved with the VLA. Especially the spatial offset between the N_HI peak and the star-clusters discussed above (cf. Fig. 5) precludes such a conclusion. This rather indicates substantial substructure below our resolution limit in the neutral phase. The larger beam of the 7 kλ-tapered data smooths out almost all resolved features in the untapered data (Fig. 7, central panel). This is evidenced by the peak column (9.5 × 10²⁰ cm⁻²) being only slightly higher than the average column (9 × 10²⁰ cm⁻²) in the untapered map. However, the tapered data reveals the diffuse lower-column gas that is spread over larger areas in the outskirts of the system. Here this gas forms a bridge between the galaxy pair and also possible tidal features towards the east and west. The N_HI morphology from this reduction appears broadly consistent with the significant detected H I in the VLA C- and D-configuration data presented in Pustilnik et al. (2001). We find a similar level of agreement between our maps and the GMRT 21-cm observations presented by Ekta et al. (2009). As described both by Ekta et al. (2009) and Pustilnik et al. (2001) in detail, the morphological features (and also the kinematics) of the H I halo gas on larger scales strongly indicate an early stage of a merger between both galaxies. The onset of this merger is hypothesised as having triggered the current star-formation episodes in both galaxies.

Comparing the morphology of the H I envelope and the extended Hα emission, it is evident that the filaments extend into regions that are characterised by lower neutral column densities or, at greater radial distance, are even devoid of detected H I; This is especially so for the eastern filament. However, the western filament appears to overlap more with the low-column gas of the bridge. Especially in Fig. 8 it appears as if the extended Hα structure starts to interact with the extended H I halo of the western galaxy. Interestingly, this boundary coincides with the fanning out of the western filament (Fig. 2). The elongation of the H I halo surrounding the eastern galaxy is oriented almost perpendicular to the direction of the filaments.

Fig. 8. Hα narrowband from MUSE (log-stretch from 0 to 1.25 × 10−15 erg s−1 cm−2) inset into a grz colour composite (Pan-STARRS) with NHI contours from VLA B-configuration data tapered at 7 kλ; contours delinate NHI = {0.5(2.5σ),1, 2, 5, 8}×1020 cm−2. The beam is 15.9″ × 14.7″ oriented at 66.4° and is indicated as a hatched ellipse.

3.2. Emission line ratios along the filaments

Our detections of Hβ and [O III] λ5007, in concert with Hα, enable us to make limited inferences on the physical state of the emitting plasma in the filaments. We show in Fig. 9 the continuum-subtracted narrowband images in [O III] λ5007 and Hβ. The images were synthesised by summing over the datacube layers at 4927.46 ± 1.25 Å and 5074.96 ± 1.25 Å observed frame for Hβ and [O III], respectively. These narrow windows are adequate, since the line emission in the filaments is barely resolved (see kinematic analysis in Sect. 3.3.1).

Fig. 9. Top panel: continuum subtracted [O III] λ5007 narrowband; contours indicate SB[OIII] = {5, 12.5, 25}×10−19 erg s−1 cm−2 arcsec−2. The cyclic colour map encodes flux densities from 0 to 5 × 10−19 erg s−1 cm−2 (colour bar) to 5 × 10−17 erg s−1 cm−2. The image has been smoothed with a Gaussian of 0.4″ FWHM. Bottom panel: continuum subtracted Hβ narrowband. The cyclic colour map encodes flux from 0 to 3 × 10−19 erg s−1 cm−2 to 3 × 10−17 erg s−1 cm−2. The image has been smoothed with a Gaussian of 0.95″ FWHM to enhance the visibility of the low-SB Hβ emission.

We analyse the emission line ratios Hβ/Hα (Sect. 3.2.1) and [O III]/Hα (Sect. 3.2.2) within two rectangular regions that cover each filament of an area where all three lines are detected; these slits are indicated in the bottom left panel of Fig. 1. Both slits cover 6″ across each filament, and the length of the eastern and western slits are 27.6″ and 19.2″, respectively. The flux distribution of Hα, Hβ, and [O III] along (Δl) and across (Δw) the filaments is shown in the top three panels of Fig. 10. These flux distributions were extracted directly from the synthesised narrowband images. Corresponding variance maps in all lines were created by summing the respective layers in the variance datacube. We measured the individual emission line fluxes by summing the emission across the filaments in bins of 1.2″ along the direction of the slit. We then calculate the ratios for each bin after correcting the individual measurements for galactic foreground emission using the Cardelli et al. (1989) extinction curve for a given A_V = 0.123 from the Galactic dust map of Schlafly & Finkbeiner (2011).

Fig. 10. Analysis of emission line ratios along the eastern (left) and western (right) filament. In the top three panels the emission in Hα, Hβ, and [O III] λ5007 is shown, now aligned such that the x-axis runs along the direction of the filament (Δl) and that the y-axis runs across the filament (Δw); the correpsonding areas have also been outlined in the bottom left panel of Fig. 1. The two bottom panels show the Hβ/Hα and [O III]/Hα ratio in bins of 1.2″ along Δl. These ratios are obtained after integrating the emission in each line along Δw. The 1-σ errors on the ratios have been computed with Eq. (2). For Hβ/Hα the Case-A recombination expectations from Storey & Hummer (1995), are indicated as dashed lines for three different temperatures (T = {1, 2, 3}×104 K in green, orange, and red); cf. Fig. 11. For Hβ/Hα the corresponding conventional Balmer decrement values are measured on the right ordinate.

3.2.1. Inverse Balmer decrement Hβ/Hα

Traditionally the Balmer decrement, Hα/Hβ, is used to measure interstellar extinction. We here use the inverse of the Balmer decrement since the statistics for the ratios of measurements for the low surface-brightness emission close to the detection limit are only well defined if the denominator is chosen to have the smaller error-bar. Then the standard 1-σ error on the ratio, Δ(Hβ/Hα), can be calculated via

$$\begin{array}{c}\hfill \mathrm{\Delta }(\mathrm{H}\beta /\mathrm{H}\alpha )=\frac{\sqrt{\mathrm{Var}(\mathrm{H}\alpha )+\mathrm{Var}(\mathrm{H}\beta ){(\mathrm{H}\beta /\mathrm{H}\alpha )}^2}}{\mathrm{H}\beta },\end{array}$$(2)

where Var(Hα) and Var(Hβ) are the variances on the extracted fluxes Hα and Hβ, respectively. Equation (2) is only valid as long as $\sqrt{\mathrm{Var}(\mathrm{H}\beta )}\lesssim 0.25\mathrm{H}\beta$ (e.g. Dunlap & Silver 1986), which is the case here.

To aid the interpretation of this non-standard ratio, we plot in Fig. 11 its behaviour as a function of temperature for a recombining plasma under Case-A and Case-B recombination scenarios. Both in Figs. 10 and 11 we also include ticks on the right ordinate that allow reading off the conventional Balmer decrement.

Fig. 11. Inverse Balmer decrement values for Case-A (dashed line) and Case-B (solid line) recombination as a function of temperature as calculated by Storey & Hummer (1995). The right ordinate measures the conventional Balmer decrement values.

From an integrated spectrum covering the star-forming sites in the galaxy (depicted in Fig. 3 in Wofford et al. 2021) we find an inverse Balmer decrement of 0.34 (or Hα/Hβ = 2.92). As can be seen in Fig. 10, both the eastern and the western filament are predominantly characterised by higher Hβ/Hα ratios (lower Balmer decrements), that is Hβ/Hα ≳ 0.35 (Hα/Hβ ≲ 2.86). Interesting features in the radial ratio profiles are the dips before and after the Hα “hot-spot” in the western filament, for which the higher ratios appear closer to the galaxy; in contrast, the eastern filament shows the highest values at its end. The mean (median) of Hβ/Hα along the eastern and western filament are 0.4 (0.38) and 0.41 (0.39), respectively.

Some of the here observed (high) low (inverse) Balmer decrements in the filaments are not compatible with standard Case-B or even Case-A recombination values (Fig. 11). The ansatz of the Case-B recombination scenario is infinite optical depth in all Lyman series lines (τ_Lyn = ∞), but zero optical depth in all other transitions, whereas for Case-A also τ_Lyn = 0. Case-B delivers realistic values for Hydrogen line ratios in the interstellar medium, since the absorption cross section for Lyman series photons is very high. Both in Case-B and in Case-A Hβ/Hα increases with increasing temperature, but the increase is more rapid for Case-A. However, even for T ≈ 3 × 10⁴ K, the highest temperature provided in the calculations by Storey & Hummer (1995), the resulting Hβ/Hα = 0.385 for Case-A (red dashed line in Fig. 10) and Hβ/Hα = 0.370 for Case-B are well below the measured ratios, 0.4 ≲ Hβ/Hα ≲ 0.5, at the beginning (end) of the western (eastern) filament.

Since the low-surface brightness filaments are observed in the low-density outskirts of the galaxy it is qualitatively conceivable that the filaments originate in gas where neither the τ_Lyn = ∞ nor the τ_Lyn = 0 approximation can be justified. Radiative transfer effects thus may alter the population levels, and neither the Case-A or Case-B approximation are reasonable anymore. Calculations for this intermediate τ_Lyn → 0 regime were performed by Capriotti (1966) and Cox & Mathews (1969); cf. Chap. 4.5 of Osterbrock & Ferland (2006). Especially at low optical depths in Lyα (τ_Lyα < 10³) Balmer decrements compatible with the range of the here measured ratios in the filaments, 2 < Hα/Hβ < 2.5 (0.4 ≤ Hβ/Hα ≤ 0.5), are possible (Fig. 4.3 in Osterbrock & Ferland 2006). Deviations from equilibrium occur in media that are ionised of fast-radiative shocks, but model calculations in low-metallicity environments favour higher Balmer decrements ≈3 (Alarie & Morisset 2019).

3.2.2. [O III]/Hα

The line ratio [O III]/Hβ is typically used to map the excitation of the plasma. Given the limited coverage of detected Hβ emission we here substitute Hβ with Hα in the denominator. The standard error on [O III]/Hα, Δ([O III]/Hα), can then be calculated in an analogous way to Eq. (2). In principle, [O III]/Hα must be understood as a lower bound to the excitation, since extinction lowers the ratio. This needs to be kept in mind, since the unusual low Balmer decrements discovered in the previous section preclude us from making definitive claims on the presence of dust within the filaments.

From the integrated spectrum covering the central star-forming sites (see Fig. 3 in Wofford et al. 2021) we measure [O III]/Hα = 1.13. For both filaments the ratio is significantly lower, with an mean (median) of 0.63 (0.65) and 0.26 (0.23) for the eastern and western filament, respectively. Such low [O III]/Hα ratios have been reported recently in the outskirts of the blue compact dwarf galaxy Haro 14 (Cairós et al. 2022). From Fig. 9 it could be already anticipated that the western filament is significantly dimmer in [O III] (see also Fig. 1 in H17). The ratio in the eastern filament is relatively flat, and only at the end (Δl ≳ 17″) of it does the [O III] flux decrease more rapidly than Hα. For the western filament the ratio appears to show an overall decrease with increasing distance from the galaxy, albeit with significant fluctuations. There appears to be no strong correlation between trends seen in Hβ/Hα and [O III]/Hα.

Interpreting the [O III]/Hα ratio is complicated, as it is regulated by multiple parameters that characterise the physical conditions of the emitting plasma. At low-densities (n ≲ 1000 cm⁻³) the emissivity of [O III] λ5007 depends only on T_e (Luridiana et al. 2015; Morisset et al. 2020), but the ratio of O²⁺ to H⁺ ions is influenced both by the oxygen abundance and the ionisation mechanism. More specifically, if O²⁺/H⁺ is held fixed then the important parameters regulating [O III]/Hα are the hardness of the ionising radiation field and the ionisation parameter (ratio between density of ionising photons to density of H atoms) for photo-ionisation, whereas in fast radiative shocks the shock-velocity, density, and the strength of the magnetic field are influential (e.g. Alarie & Morisset 2019). We recall that the western filament appears to interact with the low neutral column gas that belongs to the tidal bridge that connects the western with the eastern halo, whereas the eastern filament appears to protrude more freely out of the neutral envelope (Sect. 3.1.2 and Fig. 8). Thus, a plausible interaction with the more pristine halo gas in the western filament might lower the oxygen abundance and thereby also attenuate the [O III] line. On the other hand, this interaction could also lead to a higher neutral fraction within the western filament that in turn lowers the ionisation parameter and thereby also attenuates [O III] line.

3.3. Kinematic analysis

3.3.1. H II kinematics

In Fig. 12 we display the line-of-sight velocity field (v_los, top panel) and the velocity dispersion (σ_v, bottom panel) of the line emitting plasma. We created these maps from fits of Gaussian profiles to Hα, Hβ, and [O III] λλ4959,5007 in a Voronoi tessellated³ continuum-subtracted datacube; in the outer filaments only Hα is contributing to the kinematic signal. By experimenting we found that a S/N threshold of 8.5 and a maximum bin-size of 4 square arcseconds provided useful results, that is to say those parameters preserve the morphological appearance of the filaments while also providing a large number of high-SN bins throughout the low-SB regions; all bins in the high-SB region consist of single spaxels. Despite the morphological complexity, the emission lines are all well modelled by a single Gaussian. We show in Appendix A the line profiles and model along the main filaments. Only for a few spaxels in the central regions a weaker secondary component could be added to better fit the observed skew in some profiles. As zero-point velocity we set v₀ = c ⋅ z = 4053 km s⁻¹ and σ_v is corrected for the instrumental line spread function using the prescription provided in the MUSE user manual.

Fig. 12. Line-of-sight velocity map (top) and velocity dispersion map (bottom) of the ionised gas; subdued contours mark various Hα surface brightness levels as in Fig. 1. As zero-point velocity v = c ⋅ z = 4053 km s−1 (Moiseev et al. 2010) was adopted. The kinematic major and minor axis after Moiseev et al. (2010) are indicated with a dash-dotted and dotted line, respectively. The viewport of this figure is 1′35″ × 1′46″ as in Fig. 1. The blue dashed square indicates the 22″ × 22″ viewport centred on the central circumgalactic medium (Fig. 2) shown in the left panel of Fig. 4.

Low S/N σ_v values obtained from Gaussian fits are known to be biased high, since increasing the noise relative to the line flux artificially broadens the fits. Therefore we show in the bottom panel of Fig. 12 only σ_v’s where SB_Hα > 1.5 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻². Moreover, we only include bins in the v_LOS map where the measurement uncertainty Δv < 15 km s⁻¹. In the low-SB filaments typical values for the 1-σ uncertainties on velocity and velocity dispersion are Δv ∼ 10 km s⁻¹ and Δσ_v ∼ 20 km s⁻¹, respectively. Further, it is known for MUSE data, that velocity fields derived from narrow emission lines exhibit a characteristic striping pattern on small scales (see Sect. 3.4 in Weilbacher et al. 2015); this is due to the LSF being undersampled by the MUSE spectrograph. For the display in Fig. 12 we removed the most prominent effects of this striping pattern by slightly degrading the spatial resolution of the velocity fields. This is achieved by smoothing the computed v_LOS map with a circular top-hat kernel of 0.4″ (2 pixels) radius.

The H IIv_LOS map in the central high-SB region appears quite chaotic at first glance. Nevertheless, Moiseev et al. (2010), using SAO/SCORPIO Fabry-Perot spectroscopy of Hα in the central region and using VLA C- and D-configuration 21 cm data from Pustilnik et al. (2001), showed that the velocity field can be characterised by a disk-like north-west to south-east gradient that exhibits a distinct perturbation in the centre. Moiseev et al. (2010) modelled this central disturbance, which is shaped somewhat like a hook, as the effect of an expanding super-bubble disturbing a disk-like velocity field. The best-fitting disk was found to be aligned at a position angle of 53° east of north. We plot this kinematic axis and the orthogonal axis forced through the photometric centre (Sect. 3.1.1) in Fig. 12. It can be appreciated how the orthogonal axis appears almost symmetric in between the two filaments.

We do not observe velocity gradients along the direction of the extending filaments. As already noted in H17, the filaments seemingly inherit the velocity that is found at their base where they connect the high-surface brightness region. The northern filament appears to inherit the blue-shifts from the central perturbation, whereas the western filament appears slightly more redshifted. A conspicuous kinematic feature is the blue-shift (≈ − 35 km s⁻¹) at the western sub-branch at the tip of the western-filament. As analysed in Sect. 3.1.2, this is where the 21 cm maps indicate an overlap between Hα emitting gas and the H I bridge that connects to the western galaxy. Excluding this sub-branch, the velocity difference between the average velocities in the filament is ∼3 − 5 km s⁻¹.

The overall shear of the complete velocity field is v_shear = 1/2 × (v₉₅ − v₅) = 28.5 km s⁻¹, with v₉₅ and v₅ being the upper and lower fifth-percentile of the velocity map. By limiting the calculation of v_shear only to the high-SB region (SB_Hα > 12.5 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻²) we find $v_{\mathrm{shear}}^{>12.5}=22.9$ km s⁻¹. This can be compared to the disk model by Moiseev et al. (2010), that is characterised by a maximal rotation velocity of v_max = 28.2 km s⁻¹ and an inclination of i = 37° (see also Moiseev et al. 2015), that is the projected maximum amplitude of this disk is v_max ⋅ sin i = 17 km s${}^{-1}<v_{\mathrm{shear}}^{>12.5}$. Moiseev et al. (2010) attributed this discrepancy between disk-like and observed motions in the high-SB region to non-circular motions due to the expanding shell near the centre. For the large-scale H II velocity field the discrepancy is even more pronounced.

That neither a large velocity offset between the filaments, nor a strong velocity gradient along the filaments, nor abrupt discontinuities with respect to the central velocity field are observed may be caused by projection effects (i.e. because of material moving mostly perpendicular to our sightline). However, it may also indicate the absence of strong flows of ionised gas along the filaments. The eastern “foreshortened” filaments (cf. Fig. 2) are not characterised by an abrupt break in line-of-sight velocities. Ostensibly they just continue the velocity gradient towards positive velocities along the north-eastern direction. However, the $v_{\mathrm{shear}}>v_{\mathrm{shear}}^{>12.5}$ behaviour is caused by the redshifts of the foreshortened filaments and by the blue-shift in the western sub-branch of the western filament. Hence, The velocities observed in those extended regions are not compatible with an idealised disk and thus the motions of the ionised filamentary halo gas are not predominantly driven by the gravitational potential of the galaxy.

The velocity dispersion map in Fig. 12 (bottom panel) shows that the filaments are characterised by relatively narrow emission, with average values of ≈30 km s⁻¹. This is similar to the galaxy’s luminosity weighted average ${\sigma }_0=\sum _{x,\mathit{y}}{\mathrm{NB}}_{x,\mathit{y}}^{\mathrm{H}\alpha }{\sigma }_{x,\mathit{y}}/\sum _{x,\mathit{y}}{\mathrm{NB}}_{x,\mathit{y}}^{\mathrm{H}\alpha }=29$ km s⁻¹; here σ_x, y denote the dispersion measures in each pixel of the dispersion map and ${\mathrm{NB}}_{x,\mathit{y}}^{\mathrm{H}\alpha }$ are the flux values recorded in each pixel of the Hα narrow band image (Fig. 1). This value is in exact agreement with σ₀ derived from high-spectral resolution (R ∼ 86 000) Fabry-Perot spectroscopy (Moiseev et al. 2010, 2015). We note, however, that in the filaments σ_v/Δσ_v ≲ 2, that is there the Hα line is barely resolved with MUSE.

The ratios of large-scale ordered motions, quantified by v_max or v_shear, to localised unordered motions, quantified by σ₀, are a useful quantity to characterise the kinematical state of a galaxy (see review by Glazebrook 2013). We confirm v_max/σ₀ ≈ 1 for the disk model of Moiseev et al. (2010), and we note v_shear/σ₀ ∼ v_max/σ₀. The v/σ-ratio for SBS 0335−052E is lower that what is observed for typical disk-galaxies (v_max/σ₀ ≳ 4). A large fraction of star-forming galaxies at high-z show also v_max/σ₀ ≲ 1 (e.g. Turner et al. 2017; Wisnioski et al. 2019). Star-formation-driven feedback is deemed to be an important driver of those so-called dispersion dominated kinematics.

The zone of narrow lines at the photometric centre is surrounded by zones of elevated velocity dispersions (45−55 km s⁻¹, towards the north and south). These zones appear as double-peaked Hα profiles in high spectral resolution (R ∼ 10 000) VLT/GIRAFFE ARGUS data (Izotov et al. 2006, their Fig. 7a), but at our lower spectral resolution we here only observe a broadening. These double peaked profiles were interpreted as expanding shells. Moreover, Izotov et al. (2006) conjectured the launch of an outflow perpendicular to the disk from the intricate line profiles seen in the small ARGUS field of view (see especially Figs. 7 and 8). This is consistent with our interpretation of the filaments being large-scale effects of this outflow (see discussion in Sect. 4).

3.3.2. H I kinematics

In Fig. 13 we display the intensity weighted line-of-sight velocity field of the 21 cm signal from SBS 0335−052E & W, both for the untapered and the 7 kλ-tapered data products. The velocity maps in km s⁻¹ were computed from the moment-1 maps produced by the source finding software SoFiA (Serra et al. 2015). The relevant input parameters for our SoFiA runs were stated in Sect. 3.1.2.

Fig. 13. Moment-1 map from VLA-B configuration 21 cm data showing the intensity weighted line-of-sight velocity field of H I (top: untapered data; bottom: 7 kλ tapered data). The purple contours mark various Hα surface brightness levels as in Fig. 1. For the zero-point velocity the optical recession velocity from the eastern galaxy, v0 = c ⋅ z = 4053 km s−1 (Moiseev et al. 2010), was used.

The velocity map from the 7 kλ-tapered datacube appears quantitatively consistent with the maps from GMRT data (Ekta et al. 2009) and with the maps from VLA C- and D-configuration data (Pustilnik et al. 2001). The H I envelope encompassing both systems may be characterised by an overall shear from east to west, whose gradient is, however, not perfectly smooth. Both Ekta et al. (2009) and Pustilnik et al. (2001) argued for the envelope being comprised of two disk-like systems. Judging from our maps, but also from the maps presented in the aforementioned publications (Fig. 6 in Pustilnik et al. 2001 and Fig. 1 in Ekta et al. 2009), this reading appears overly simplistic. Focusing here on the eastern galaxy, another gradient from redshifts in the south-east to blue-shifts the north-west is perceivable. The alignment of this gradient along the directions of the filaments is suggestive of causal relationships between ionised and neutral phase. At the higher spatial resolution of the untapered datacube the velocity fields exhibit an even higher complexity, and simple disk-like gradients appear indiscernible. In the top-left panel of Fig. 14 we show a magnified view of the velocity field of the eastern galaxy; an overall north-west to south-east shear may be envisioned, but this velocity gradient is significantly warped.

Fig. 14. Comparison between H I and H II kinematics on small scales. The H I velocity field (top left panel) is a zoomed-in view of the velocity field shown in the top panel Fig. 13. The H II velocity field (top right panel) is the same as in Fig. 12, but regions outside of detected H I have been subdued. In the bottom left panel we show the H II velocity field convolved with a 2D Gaussian (6.1″ × 4.6″, position angle −23.6°) to mimic the VLA B-configuration observations beam smearing (see text). The bottom right panel shows the smoothed Hα velocity field resampled onto the native grid of the VLA B-configuration cube.

As already noted in Sect. 3.1.2, the faintest H I signals from the outskirts are often only single channel detections. This can also be appreciated from the channel maps in Appendix B. While our data does not allow for robust higher-order measurements, these channel maps supplement the view of high kinematic complexity in the H I on kpc scales in the system. From H I morphology and kinematics alone it is not possible ascribe tidal- or feedback effects as the source of this complexity. We thus proceed by comparing H I and H II kinematics quantitatively in the next section.

3.3.3. Comparison between H I and H II kinematics

We compare the ionised gas velocity field derived from the MUSE data (Sect. 3.3.1) to the H I velocity field from the VLA observations (Sect. 3.3.2). For this comparison we need to account for the lower spatial resolution and the coarser spatial sampling of the H I data products. To this aim we follow the method used by van Eymeren et al. (2009a,b) and convolve the H IIv_los map (Fig. 12) with a 2D Gaussian, whose position angle, major-, and minor-axis are matched to the beam of the respective 21 cm data. We next re-bin this convolved velocity field to the native resolution of the VLA datacubes. The complete process is illustrated in Fig. 14, with the top right panel displaying the H II velocity field as observed, with the bottom left panel displaying it after convolution, and with the bottom right panel displaying this convolved field after re-binning. For the 7 kλ-tapered dataset we omit the intermediate step and show only the resulting resampled velocity field (bottom panel of Fig. 15).

Fig. 15. Comparison between H I and H II kinematics on larger scales; we here display a zoomed-in view of the 7 kλ tapered 21 cm data shown in the bottom panel of Fig. 13. The top panel shows the first moment from the H I data, the middle panel shows the ionised gas velocity field, with regions outside of H I detections being subdued, and the bottom panel shows the H II velocity field convolved with a 2D Gaussian (15.9″ × 14.7″, position angle 66.4°) to mimic the beam smearing and resampled onto the native grid of the H I observations.

Qualitatively, both the velocity fields from the untapered and the 7 kλ-tapered H I data show a high level of congruence with the degraded and resampled H II velocity fields (Figs. 14 and 15). Our method ignored the effect of intensity variations on the resulting resolution-degraded and re-binned H II velocity fields. If we instead convolve each MUSE datacube layer with the beam-imitating 2D Gaussian prior to the kinematic fitting, then the resulting degraded velocity fields show almost no spatial variation, as they mostly trace H II velocities stemming from the brightest regions. Hence, the unknown intrinsic 21-cm signal from the neutral phase does not mimic the extreme intensity variations seen in the H II emitting plasma.

The degraded H II velocity fields appear more smoothly varying compared to the H I measurements, which is partly related to the low-SN of the H I data. We analyse the quantitative differences between H I velocity fields, v_HI, and the resolution degraded and re-sampled H II velocity fields, $v_{\mathrm{HII}}^{\mathrm{rep}}$, in Fig. 16. There we display $\mathrm{\Delta }v=v_{\mathrm{HI}}-v_{\mathrm{HII}}^{\mathrm{rep}}$ maps both for the untapered and the 7 kλ-tapered datasets. These Δv-maps exhibit mostly an amplitude of absolute differences ≲10 km s⁻¹ and thus confirm the qualitative impression of overall congruence already anticipated from Figs. 14 and 15. The similarity between both velocity fields manifests in small mean (median) absolute differences of 0.6 km s⁻¹ (0.2 km s⁻¹) and 1.1 km s⁻¹ (0.6 km s⁻¹) for the untapered and 7 kλ-tapered velocity fields, respectively. There are, nevertheless, some notable differences. In the untapered data we find Δv ≳ 15 km s⁻¹ near the centre. This $v_{\mathrm{HI}}>v_{\mathrm{HII}}^{\mathrm{rep}}$ zone is in close vicinity to the most prominent “hook-like” disturbance of the H II map. A causal connection between H I−H II velocity offsets and the expansion thus suggests itself. In the 7 kλ-tapered data we find a prominent rim $v_{\mathrm{HI}}>v_{\mathrm{HII}}^{\mathrm{rep}}$ along the south-west, parallel to the minor axis. The origin of this feature is not obvious. Both high-Δv zones in both datasets are related to the most prominent complexities of the H I velocity fields that inhibited the identification of simple disk-like velocity gradients in the first place (Sect. 3.3.2).

Fig. 16. Maps showing $\mathrm{\Delta }v=v_{\mathrm{HI}}-v_{\mathrm{HII}}^{\mathrm{rep}}$, where vHI are the H I line-of-sight velocities and $v_{\mathrm{HII}}^{\mathrm{rep}}$ are the resolution matched and re-projected H II line-of-sight velocities. The left panel shows the comparison with the untapered VLA-B configuration data, i.e. the subtraction of the velocity field in the lower-right panel of Fig. 14 from the upper-left panel of Fig. 14, while the right panel shows the comparison for the 7 kλ-tapered data, i.e. the subtraction of the velocity field in the lower panel of Fig. 15 from the upper panel of Fig. 15.

4. Discussion

The presented observational evidence in the previous section leads us to explain the observed features in the circumgalactic halo of SBS 0335−052E as limb-brightened edges of a conical outflow phenomenon. We briefly summarise the main clues towards the suggested interpretation before making an attempt to quantitatively illustrate this scenario with a didactic toy model.

VLT/MUSE reveals two faint ionised filaments of delicate morphology that extent up to 15 kpc to the north-west of the starburst (Fig. 1). The orientation of a perceived symmetry axis that separates the rays of the filaments is perpendicular to the shearing direction of the large-scale velocity field (Fig. 12). This large-scale shear lines up with the major axis of the disk-model derived from SAO/SCORPIO Fabry-Perot observations in the central high-SB_Hα region (Moiseev et al. 2010, 2015); therefore, the symmetry axis of the bifurcating filaments is nearly parallel with the minor axis of this model (Fig. 12). The ionised filaments protrude away from the elongation axis of the H I envelope (Fig. 8). This envelope also exhibits overall disk-like shearing motions (Fig. 15). These observational facts are reminiscent to prototypical conical outflows in disk-galaxies that launch such cones from their central starburst outwards from the disk plane (Veilleux et al. 2005; Bland-Hawthorn et al. 2007; Nelson et al. 2019).

We explore the possibility of this scenario by setting up a heavily idealised astrophysical structure with parameters suggested by the observations. Our toy model geometry is a simple cone of opening angle θ and height h. The apex of this cone is fixed at the centre of a disk at position angle ϑ and inclination i. Thus the observed opening angle, θ_P, is a projection of the true cone opening θ = 2 × arctan[sin(i)×tan(θ_P/2)]. We assume that the Hα filaments are confining this cone. We measure θ_P = 34° as the angle between the filaments from the Hα NB (Fig. 1) and the disk-model of Moiseev et al. (2010) which fits i = 43°, hence θ = 27°. We note, however, that conical winds can be tilted and asymmetric with respect to the centre; our value for θ thus represents an approximate working hypothesis. We also adopt ϑ = 52° from Moiseev et al. (2010). We fix l_p = 10 kpc for the projected length of the filaments from Fig. 1, where we ignore for now the threadlike extension. The height of the inclined cone is then h = l_p ⋅ cos(θ/2)/sin(i) = 16.2 kpc.

Our estimate of the height h may be translated into a velocity requirement that the hypothesised wind from a stellar-population of age t_* must have to blow out such a structure. The oldest star-clusters in the north (SSC 5 & SSC 6; Fig. 5) have ages t_* ≲ 15 Myr (Reines et al. 2008; Adamo et al. 2010). Accounting for the delay after which supernovae start to inject kinetic energy, t_SN ≈ 5 Myr, we have v_wind = h/(t_* − t_SN) = 10 kpc/10 Myr × cos(θ/2)/sin(i)≈1000 km s⁻¹ × cos(θ/2)/sin(i) = 1620 km s⁻¹. This rough limit, which neglects variations in energy input and flow-speed, is broadly consistent with predictions for velocities in the tenuous hot phase (T > 10⁶ − 10⁷ K, n ∼ 10⁻³ − 10⁻⁴ cm⁻³) of galactic winds. This “wind-fluid” is believed to carry most of the momentum and energy of an outflow, and at larger radii even most of the mass (Schneider et al. 2020). If the volume of the cone, V_hot = π/3 ⋅ tan²(θ/2)⋅h³ = 256 kpc³, would be filled uniformly with this fluid, then its mass would comprise M_hot ∼ 6 × 10⁶ − 6 × 10⁵ M_⊙. This requires loading rates, ${\mathrm{\Lambda }}_{\mathrm{hot}}=(M_{\mathrm{hot}})/({\dot{M}}_{\mathrm{SFR}}·\mathrm{\Delta }t)$, of order unity if all the mass is asserted to be loaded into the wind via star-formation over the last Δt = 10 Myr given the determined Ṁ_SFR ≃ 1 M_⊙ yr⁻¹ in SBS 0335−052E (Reines et al. 2008; Thompson et al. 2009). This rough limit for Λ_hot is higher than the expected average mass-return from stellar winds and supernovae in stellar populations, Λ = 0.1 (Leitherer et al. 1999), that is also used as the injected mass fraction for the hot phase in computer simulations of wind phenomena (e.g. Schneider & Robertson 2018; Schneider et al. 2020). However, a significant fraction of the fluid will consist of shock-heated ambient ISM, hence the required mass loaded into the interior of our cone does not appear unphysical.

The hot fluid, which does not emit in optical recombination lines, is assumed to be surrounded by the warm ionised medium that fills the wall of the cone and leads to the observed structure in Hα. Denoting the thickness of the walls as t and writing Δh = 2 ⋅ t ⋅ cos(θ/2) the volume inhabited by the diffuse H II gas is V_HII = π/3 ⋅ tan²(θ/2)⋅((h+Δh)³−h³−Δh³). Figuratively, V_HII is the volume between two identical cones pushed into each other (∝(h + Δh)³ − h³), but ignoring the volume between the apexes of the inner and the outer cone (∝Δh³). From the Hα NB we estimate t = 1.5 kpc and hence V_HII = 142 kpc³.

In the spirit of a toy-model we assume the plasma filling V_HII is in recombination equilibrium, fully aware that deviations from equilibrium are expected in shocked gas (e.g. Allen et al. 2008; Morisset et al. 2020) or due to radiative transfer effects (Sect. 3.2.1). An ionised hydrogen plasma of density n (≡n_e = n_p, i.e. no electrons from other species) and temperature T filling a volume V emits Hα emission at luminosity

$$\begin{array}{c}\hfill L_{\mathrm{H}\alpha }={ϵ}_{\mathrm{H}\alpha }(T)·n^2·V\end{array}$$(3)

in recombination equilibrium, where ϵ(T) denotes the Hα equilibrium emissivity that follows from atomic data. For T = {1.0, 1.5, 2.0}×10⁴ K in case-A recombination ${ϵ}_{\mathrm{H}\alpha }^{\mathrm{A}}(T)=\{2.3,1.5,1.1\}\times {10}^{-25}$ erg s⁻¹ cm³ (Storey & Hummer 1995); the small density dependence on ϵ can be neglected for our purposes. Using a hand-drawn aperture that encompasses the filamentary structure while not covering inner high-SB spurs and shell-like structures (thin dotted line in Fig. 2), we estimate F_Hα ≈ 2.6 × 10⁻¹⁵ erg s⁻¹ cm⁻² as the flux from the supposed cone. From Eq. (3) then n_A = {3, 4, 5}×10⁻² cm⁻³ follows for T = {1, 1.5, 2}×10⁴ K in case-A. We recall from Sect. 3.2.1 that Case-A recombination appears more compatible with the observed Balmer decrement in the filaments than case-B recombination, noting here that ${ϵ}_{\mathrm{H}\alpha }^{\mathrm{B}}(T)<{ϵ}_{\mathrm{H}\alpha }^{\mathrm{A}}(T)$ and that ${ϵ}_{\mathrm{H}\alpha }^{\mathrm{B}}/{ϵ}_{\mathrm{H}\alpha }^{\mathrm{B}}$ is only mildly temperature dependent. Therefore we can approximate ${ϵ}_{\mathrm{H}\alpha }^{\mathrm{B}}\approx 0.65{ϵ}_{\mathrm{H}\alpha }^{\mathrm{A}}$ and hence n_B ≈ 1.24 ⋅ n_A would be required for case-B.

The required densities in the walls of the cone are only a factor of two to three lower than the canonical n = 0.1 cm⁻³ for diffuse warm ionised gas in the circumgalactic medium. They are consistent with the expectations for the T ∼ 10⁴ K phase of galactic winds. The correspondence between the derived parameters in our idealised setup and the physical characteristics expected in outflows appears encouraging and supportive of the idea that we indeed observe a 15 kpc outflow cone emanating from a tiny starburst galaxy. We remark that the mass of ionised hydrogen in the toy-model, M_HII = 1.5 − 1.7 × 10⁸ M_⊙, requires a mass-loading factor Λ_HII ≳ 10 if this mass has been loaded into the wind over the last 10 Myr. Neglecting the cold phase, which is expected to provide only a minor contribution to the total mass of outflowing material, we arrive at a total mass loading factor of Λ = (M_hot + M_HII)/(Ṁ_SFR · Δ>t) ≳ 10. Given the low-stellar mass of SBS 0335−052E, M_⋆ ≈ 6 × 10⁶ M_⊙ (Reines et al. 2008), the required Λ is in line with theoretical expectations (Pandya et al. 2021) and observed trends for the other dwarf starbursts (Collins & Read 2022, their Fig. 4). Interestingly, the ratios Λ_hot/Λ ≈ 0.1 and Λ_HII/Λ ≈ 1 are in line with the expectations of the wind models by Pandya et al. (2021). We caution that our results have to be taken cum grano salis, as we assumed that the diffuse Hα emission in the cone is caused by pure Case-A recombination. As explained above, slightly higher densities would be required for the Case-B scenario, thus the mass in the warm phase would be slightly lower. However, we have shown in Sect. 3.2.1 that the Balmer-decrement in the filaments is suggestive of neither Case-A or Case-B being appropriate. Ultimately photo-ionisation models are required to obtain more stringent constraints on the observed mass in the warm phase.

We now explore the observable Hα morphology of our toy-model cone-wall structure numerically. To this aim we compute the emission on a cubic 3D lattice with 500³ volume cells that cover a 25³ kpc³ volume, that is each cell has a volume of V_cell = 50³ pc³. Then all grid cells C are set to emit Hα according to Eq. (3) with V_cell when their centre-coordinates are within the walls of the cone. Formally, this is the set of cells C = {x′,y′,z′:f(x′,y′,z′|h, θ)≤0 ⊕ f(x′,y′,z′|h + Δh, θ)≤0}, where ⊕ denotes the exclusive disjunction and f(x′,y′,z′|h, θ) = x^′2 + y^′2 − z^′2tan²(θ/2) is the implicit parametric equation for the surface of a single cone and x′, y′, and z′∈[0, h] are the coordinates in a rotated frame that is tilted along ϑ and inclined along i using standard intrinsic Euler rotations. We set up the model grid with cell-coordinates x, y, z such that the y-axis is parallel to the sight-line. We project along this axis by summation after dividing by the surface area of a cell ($6V_{\mathrm{cell}}^{2/3}$). This provides a map of surface luminosities (erg s⁻¹ kpc⁻²). This map is transformed into observable fluxes in erg s⁻¹ cm⁻² on a grid of 0.2″ MUSE-like spaxels for the given distance to galaxy. Prior to this re-binning the received fluxes are convolved with a Gaussian of 1.2″ FWHM, which was the average DIMM seeing in our observations (Sect. 2.1). Finally, we add noise to the re-binned image consistent with the noise determined in the outskirts of the Hα NB (σ_Hα = 1.8 × 10⁻¹⁹ erg s⁻² cm⁻²). We show the result of this exercise for the observationally constrained parameters from above (θ = 27°, i = 37°, ϑ = 52°, h = 16.2 kpc, and t = 1.5 kpc) and (n, T) = (0.05 cm⁻³, 2 × 10⁴ K) in Fig. 17. To provide visual intuition regarding the projected morphology for different geometrical parameters we study the effects of varying θ, ϑ, and i in Appendix C; we note that the projection of the model is the same if the cone is inclined towards us or away from us. Furthermore we examine the effects of varying T and n on the overall observability in Appendix C. We also verified the consistency between numerical computation and the analytical results presented beforehand.

Fig. 17. Toy model to explain the NW filaments as limb-brightened edges of an outflow cone. We assume a cone filled with a hot wind fluid that is surrounded by walls filled with ionised hydrogen in recombination equilibrium (T = 2 × 104 K, n = 0.05 cm−3). Details on the setup and the model parameters that are derived from the observations are given in the text. The left panel shows a projection of the resulting Hα emission (as surface-luminosity in erg s−1 kpc−1) from this model along the Y-axis of the model box; this axis is chosen to be parallel to the line of sight. The right panel shows how the Hα emission from this model would appear in our MUSE data under the assumption that the model is at the distance of SBS 03352−052E. Here we overlay some Hα contours from Fig. 1 (dotted white lines) to indicate the degree of congruence between the simulated observations of this simple toy model and reality.

We overlay in Fig. 17 the two faintest Hα SB isocontours from Fig. 1 to visualise morphological similarities and differences between calculation and observation. These contours were aligned by eye to maximise congruity. The resemblance between the MUSE mock data of the cone-wall model and observed reality is readily apparent. An interesting feature of the model is the brightening along the walls with increasing distance from the apex. This radial brightening is caused by more wall material being intersected perpendicular to the sightline closer to the base of the cone. In reality we observe a brightening with increasing distance only in the western filament (Hα “hot-spot” in Fig. 2); however, compared to the smooth radial gradient in the model this feature emerges more abruptly. In fact, the observed filaments brighten towards the central high-SB zone. Density or temperature variations, as well as variations of the ionisation mechanisms, may be invoked to explain the differences between model and observations and to produce localised features such as the Hα hot-spot.

The observations reveal further morphological sub-structure that is not reproducible with our simple purely geometrical setup. Most prominently, the cone-wall structure cannot explain the threadlike extension and the fanning out of the western filament (Fig. 2). Parts of these filamentary sub-structures coincide spatially with the region where only one wall of the cone would be intersected by the sight-lines. Zones of higher density or lower temperature in the wall section or lumps of gas that have been advected by the hot fluid could lead to brighter Hα patches. Alternatively, these features may be interpreted as the terminating cloud-wind interface where dense gas has been compressed and pushed up by the hot wind fluid and now cools radiatively.

State-of-the-art numerical models of starburst-driven winds (Schneider & Robertson 2018; Schneider et al. 2020) produce small-scale density and temperature variations in wind structures that all have been neglected in our didactic scenario. Sub-structure within the hot wind fluid arises due to interplay between the different gas-phases; cold and warm material swept up from the ISM may survive within the hot flow whereas density variations in the hot phase might lead to condensation and radiative cooling. Substructures in the HII emitting phase produce turbulent mixing layers between the fast moving tenuous fluid and the slower-moving denser warm-phase. The emergence of a conical structure appears, however, a common feature of wind phenomena. This “beaming” is believed to be a purely hydro-dynamical collimation effect due to the intrinsic density distribution in the gaseous disk surrounding the energy and momentum injecting star-formation sites (Nelson et al. 2019). In this respect it appears interesting to note that Micheva et al. (2019) advocate to interpret θ as the opening angle of a Mach cone, that is the surfaces of the cone are created by a shock-front from the overlapping sound waves in a supersonic flow. The opening angle is then related to the Mach number ℳ, that expresses the velocity of the flow as multiples of the speed of sound via ℳ = 1/sin(θ/2). For our cone we have M ≃ 4, that is to say in this interpretation the wind fluid is required to be super-sonic. Intriguingly, M ≃ 4 flows are again consistent with sophisticated model expectations for the tenuous hot wind fluid.

Obviously our toy-model is static and neglects kinematics. A kinematic model requires many assumptions and is beyond the scope of our current analysis. Still, some important kinematical features need to be discussed qualitatively, especially since the absence of a strong velocity gradient and the narrow line widths appear incompatible with the idea that the H II phase is being accelerated in an outflow. Projection affects surely may influence the observed line-of-sight velocities and velocity dispersions given that we look at the cone from the side. However, under the assumption of isotropic random motions the combination of ordered motions and inclination can only broaden the line. This effect is known to as “mixing term” in observational studies of velocity dispersions in disk galaxies (see, e.g. Sect. 2.1 in Bouché et al. 2015). In this scenario then the narrowness of the lines would indeed reflect an intrinsic narrow velocity dispersion. A projection effect is thus only expected if the distribution of velocities is larger parallel to the walls of the cone. Similarly, the lack of an observed velocity gradient along the walls of the cone may be an projection effect. Nevertheless, the large-scale wind simulations by Schneider et al. (2020) show that the H II emitting gas at larger radii is moving at fairly constant velocities. Interestingly, we do not observe strong velocity differences between H I and H II at the base of the cone in the 7 kλ-tapered data (Fig. 15 and right panel of Fig. 16). We thus speculate, that the walls of the cone are H I halo gas that is advected into the flow of the hot fluid. During this advection there might only be a short window of time before the dense warm medium, which is likely in the form of clumps, is being transformed into a more dilute wind fluid. If this is the case, then we expect that here observed emission of warm phase inherits mostly the velocity dispersion of the neutral gas. In this scenario the narrow lines appear consistent with the detection of H I emission in only 2−3 channels (see Fig. B.2), which corresponds to ∼20 − 30 km s⁻¹.

5. Summary and conclusions

We have presented new MUSE and VLA B-configuration 21 cm observations of the extremely low-metallicity blue compact starburst galaxy SBS 0335−052E. Our observations reveal the continuation of two ionised filaments in Hα emission towards the north-west of the galaxy’s main stellar body (Sect. 3.1.1 and Fig. 1). The onset of this structure was observed with MUSE observations centred on the galaxy and described in H17. The morphology of the newly detected prolongation can be described as thread-like lacy fragments. These somewhat detached portions of extended Hα emission are of lower SB (SB_Hα ≈ 1.5 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻²) than the inner straight portions of the filaments (SB_Hα ≈ 3 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻²). The rays of the inner filaments are aligned at a projected angle of θ_P = 43° relative to each other. Diffuse ionised gas with a smoothly declining brightness profile is found up to ∼25″ towards the south of the galaxy.

The inner parts of the filaments are detected in Hα, [O III] λ5007, and Hβ. The flux ratios between those three lines are different compared to typical ISM ratios in the vicinity of star formation. The inferred average Balmer decrement is Hα/Hβ = 2.5 for both filaments, whereas [O III]/Hα = 0.65 and [O III]/Hα = 0.26 for the eastern and western filament, respectively. The unusual low Balmer decrement, which does not agree with standard recombination cascade expectations, is deemed to be caused by radiative transfer effects in the low-density halo gas that alter the populations of the contributing energy levels compared to the classical calculations (Sect. 3.2.1; Osterbrock & Ferland 2006).

Our VLA B-configuration observations enable us to resolve H I at unprecedented spatial resolution within and in the vicinity of the starburst. The peak of the 21 cm signal is spatially offset by ∼1″ to the west from the main super-star-cluster complexes. At lower N_HI (5 − 10 × 10²⁰ cm⁻²), we observe morphological correspondences between neutral and ionised phases. In particular, we observe a bent tail that lines up in the direction of the eastern foreshortened filaments. Our observations are thus mapping the entanglement and interactions between the neutral and ionised phases in feedback-driven gas.

Tapering the VLA B-configuration data at 7 kλ in the UV plane in the imaging process enabled us to map the neutral halo around the interacting galaxy pair SBS 0335−052E and SBS 0335−052W. Our results are consistent with the previous analysis of 21 cm data from VLA C- and D-configuration (Pustilnik et al. 2001) and GMRT (Ekta et al. 2009) observations. The Hα emitting filaments extend outwards from the east to west elongated neutral halo around SBS 0335−052E. However, there is significant overlap between halo gas and the western filament. Especially the coincidence of the Hα hot spot and the fanning-out of the western filament with the extended tidal halo gas from the western galaxy appears intriguing.

Our kinematical analysis of the ionised gas shows that the line emission from the filaments is very narrow, σ_v ≈ 20 − 30 km s⁻¹, and in fact barely spectrally resolved in the MUSE data. We do not map velocity gradients along the filaments. Zones nearer to the star clusters appear in the MUSE data at high σ_v (45 − 50 km s⁻¹), but they have previously been resolved into double-peaked profiles with VLT/GIRAFFE (Izotov et al. 2006). The seemingly chaotic appearance of the H II line-of-sight velocity map near the centre has also been previously described as a kinematic disturbance of an expanding shell superimposed onto a disk-like velocity field (Moiseev et al. 2010). The position angle of the Moiseev et al. (2010) model velocity field, ϑ = 52°, is parallel to the shearing velocities mapped at lower SB levels and in H I. Interestingly, this gradient is aligned perpendicular to the symmetry axis in between the bifurcating filaments.

The line-of-sight kinematics of the neutral phase exhibit a significant degree of complexity, both in the untapered and in the 7 kλ-tapered data products. The complexity on large scales is likely due to tidal effects from the encounter of the two galaxies. However, an interaction with the large-scale wind may also cause deviations from simple disk-like motions in and around SBS 0335−052E. On smaller scales, the velocity gradient appears significantly warped. We find that H I and H II kinematics are broadly in agreement on large and small scales. Nevertheless, there appear to be zones where the bulk motions of the neutral phase are redshifted by ∼15 − 20 km s⁻¹ with respect to the velocities of the ionised phase. These velocity differences are deemed to be caused by different responses of the neutral and ionised phases to injected energy and momentum from stellar winds and supernova explosions.

The observational evidence assembled in our analysis merits an explanation as to the unprecedented detection of more than 10 kpc long Hα-emitting filaments as limb-brightened edges of a conical outflow; this new hypothesis supersedes the initial attempt at explaining the morphology in H17. A toy model of a cone structure, with geometrical parameters read off the MUSE Hα narrow band, is commensurable with the expected physical reality of a starburst-driven wind. Our setup is described by an inner cone that is filled with a hot wind fluid and a surrounding wall filled with ionised plasma. In order to be compatible with the dimensions of the cone, the hot wind fluid would be required to flow at radial velocities ∼1600 km s⁻¹. This hot wind fluid could be loaded easily with mass-loading factors Λ_Hot ≲ 0.1. While there is some degeneracy between density and temperature in order to match the flux of the ionised phase, we generally require very low densities of n ≲ 0.05 cm⁻³. If the ionised phase is purely gas expelled from the starburst, then mass-loading factors Λ ≈ Λ_HII ≳ 10 would be required. A mock MUSE observation of the toy model clearly shows the limb-brightening effect at surface brightness levels similar to the observations. Obviously, our didactic setup fails to reproduce the morphological substructure observed in reality, especially the thread-like lacy prolongation at the largest projected distances.

Given the small size of SBS 0335−052E’s main stellar body, the dimensions of its outflow structure revealed in the ionised phase appear extreme. Similar-sized ∼10 − 20 kpc outflows are known from more massive and larger galaxies. The finding reported here was only possible because of the unprecedented sensitivity of MUSE for low-SB line emission. Currently it is unknown whether such structures are common around other metal-poor compact starbursts. One unique characteristic of the SBS 0335−052 system is its large neutral gas reservoir, and we may speculate here that the size of the ionised structure is connected to the size of the halo, since a significant fraction of the Hα-emitting gas may have been ionised in situ.

In our analysis and discussion, we made no attempt to constrain the actual ionisation mechanism. Photo-ionisation from the starburst and the X-ray emitting wind fluid, as well as radiative shocks between wind fluid and neutral halo gas are expected to contribute. While the three emission lines in this study do not provide discriminating power with respect to these mechanisms, they can be used in connection with photo-ionisation or shock models to design deeper observations that may constrain the ionisation mechanisms. We believe that this is an intriguing avenue for future investigations.

Acknowledgments

We thank the anonymous referee for their valuable contribution to the manuscript. We express our gratitude to the ESO Office for Science in Santiago for funding research internships for H. Salas and C. Moya-Sierralta to contribute to this project. J.M. Cannon and J. Inoue are supported by NSF grant AST-2009894. J.M. Cannon and B. Koenigs acknowledge support from Macalester College. M. Hayes is fellow of the Knut & Alice Wallenberg Foundation. C. Moya-Sierralta acknowledges support from the Becas-ANID scholarship #21211528. We acknowledge the use of the following software packages for the analysis and exploration of the data presented in this paper: Astropy (Astropy Collaboration 2018), SciPy (Virtanen et al. 2020), NumPy (Harris et al. 2020), Matplotlib (Hunter 2007), GNUastro (Akhlaghi 2018), LSDCat (Herenz & Wisotzki 2017), SoFiA (Serra et al. 2015), cmasher (van der Velden 2020), QFitsview (Ott 2012), fips (Kornilov & Malanchev 2019), pyneb (Luridiana et al. 2015; Morisset et al. 2020), and ds9 (Joye & Mandel 2003). Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 096.B-0690 & 0104.B-0834, and based on observations made with the Karl G. Jansky Very Large Array (VLA), project number 17B-234. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Based on observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA) and the Canadian Astronomy Data Centre (CADC/NRC/CSA).

References

Appendix A: Emission line profiles along the filaments

In Sect. 3.3.1 we stated that the emission lines in the outer low-SB regions are adequately modelled by a single Gaussian profile. To back this claim we here visualise the line profile in the main filaments. To this aim we extract 2D spectra along the two pseudo-slits shown in the bottom-left panel of Figure 1. From these extractions we show in Figure A.1, the Hα profiles alongside the fitted Gaussian models and the residuals. It can be appreciated from this Figure in that the single Gaussian is an adequate model and that the line is extremely narrow. To further explore the simple morphology and narrowness of the line profiles we encourage the interested reader to obtain the fully reduced datacube that is released with this publication.

Fig. A.1. Hα emission line profiles extracted along the pseudo slits shown in the bottom left panel of Fig. 1 (top panels) alongside the fitted Gaussian models (centre panels) and the residuals (data – model, bottom panels). The left panels show the results for the eastern filament, whereas the right panels show the result for the western filament and the horizontal extend is the same as in Fig. 10. The x-axis is the cross dispersion axis in units of MUSE spaxels (1 spaxel = 0.2″) and the y-axis is the dispersion axis (1 pixel = 1.25 Å). The colour bar on the left corresponds to the two top panels, whereas the colour bar on the right corresponds to the bottom panel.

Appendix B: 21 cm channel maps

We supplement our analysis with channel maps from the VLA B-configuration datacubes. Figure B.1 shows seven channel maps from -24 km s⁻¹ to 37 km s⁻¹ zoomed onto SBS 0335-052E. Figure B.2 shows 7 channel maps from -65.0 km s⁻¹ to 37.5 km s⁻¹ containing both SBS 0335-052E and SBS 0335-052W. Velocity offsets in both figures are defined with respect to the optical recession velocity v₀ = c ⋅ z = 4053 km s⁻¹.

Fig. B.1. Channel maps from the untapered VLA-B configuration datacube. Each channel is shown with a linear stretch from 0 to 6σ, where σ = 2.5 × 10−4 Jy beam−1 is determined from fitting the negative distribution of pixels in each channel. Contours are drawn at 2σ (dashed), 4σ, and 6σ. The bottom right panel shows the same viewport (46″/12.7 kpc along each axis), but with the Hα NB image and the total NHI contours from Fig. 1.

Fig. B.2. Channel maps from the 7 kλ-tapered datacube, similar to Fig. B.1; here σ = 4 × 10−4 Jy beam−1. The bottom right panel shows the viewport (210″ × 114″/58.1 kpc × 31.5 kpc) over the Pan-STARRS false-colour image with the NHI contours from Fig. 8.

Appendix C: Cone model exploration

Sect. 4 discusses the ionised cone-wall structure for parameters that are derived from the observations. However, it appears instructive to visualise projections of this structures for different parameters. In Fig. C.1 we fix the opening angle θ = 27 from Sect. 4, but vary ϑ and i. In Fig. C.2 we leave i = 37° and ϑ = 52° as in Sect. 4, but vary the opening angle. While for most of the parameters the cone geometry is apparent, projections at low inclination or large opening angles appear more disk-like.

Fig. C.1. Surface-luminosity projections of the cone-wall structure from Sect. 4 for various position angles ϑ and inclinations i as indicated in the panels.

Fig. C.2. Surface-luminosity projections of the cone-wall structure from Sect. 4 for various opening angles θ as indicated in the panels.

From Eq. (3) it is apparent that there is a degeneracy in observable flux between T and n. To obtain a feel for the effects of those parameters on the observable morphology, we explore variations of T at fixed n in Fig. C.3 while in Fig. C.4 we fix T but vary n. As in Sect. 4 we use case-A recombination emissivities for Hydrogen computed by Storey & Hummer (1995).

Fig. C.3. Simulated MUSE Hα observation of the cone-wall structure from Sect. 4 at the distance of SBS 0335-052E for fixed n = 0.05 cm−2 but varying T as indicated in the panels.

Fig. C.4. Simulated MUSE Hα observation of the cone-wall structure from Sect. 4 at the distance of SBS 0335-052E for fixed T = 20000 K but varying n as indicated in the panels.

All Figures

Fig. 1. H II and H I emission as seen with MUSE and VLA around SBS 0335−052E. Top panel: continuum-subtracted Hα narrowband image. Fluxes are encoded via an asinh-scaled cyclic colour map from 0 to 10−18 erg s−1 cm−2 in the low surface-brightness, and from 10−18 erg s−1 cm−2 to 10−16 erg s−1 cm−2 in the central high surface-brightness region. The cyan solid contours demarcate Hα isophotes at SBHα = {1.5, 2.5, 5, 12.5}×10−18 erg s−1 cm−2 arcsec−2 while the white dotted contour demarcates the limiting surface-brightness of SBHα = 7.5 × 10−19 erg s−1 cm−2 arcsec−2. The image has been smoothed with a Gaussian of 0.4″ FWHM. Subdued grey contours indicate HI column densities NHI = {5, 10, 20, 30, 40}×1020 cm−2 from our VLA-B configuration observations; these contours are displayed more prominently in Fig. 4. The dashed white contour demarcates the 2.5σ = 2.5 × 1020 cm−2 detection limit. The VLA-B configuration beam, 6.1″ × 4.6″ oriented 23.6° west of north, is indicated via a grey shaded ellipse in the top left. The displayed field of view is 1′35″ × 1′46″ corresponding to 26.2 kpc × 29.2 kpc in projection. The dashed violet square in the centre indicates the viewport used in Fig. 4. North is up and east is to the left. Bottom panels: signal-to-noise of the three principal emission lines (Hα, Hβ, and [O III] λ5007 in the left, centre, and right panel, respectively) in which extended ionised gas is detected (see text for details). The maps are displayed with a linear cyclic colour map from 0 to 90 and from 90 to 900. White dotted contours mark out regions with SN > 8. The chartreuse dashed rectangles in the Hα panel outline the rectangular regions used for the line ratio analysis in Sect. 3.2.

In the text

	Fig. 2. Annotated Hα surface brightness contours (contour levels SBHα = {1.5, 2.5, 5, 12.5, 100, 700}×10−18 erg s−1 cm−2 arcsec−2). In the background the Hα narrowband is shown with a log-stretch to 5 × 10−16 erg s−1 cm−2. The thin-dotted line surrounding the filaments to the north-west demarcates the aperture that is used to estimate the total Hα flux from that structure (Sect. 4).
In the text

Fig. 3. Hα surface brightness profile in the south, i.e. not affected by the filaments in the north. The points show the observed profile in circular annuli and error-bars are 2σ. The black line indicates the double exponential profile fitted to the outer (r > 5″) profile, while dot-dashed and dashed lines indicate the individual components with scale-lengths r1 = 1.62″ (447 pc) and r2 = 8.94″ (2.47 kpc). The grey line indicates the light profile of the central (r ≤ 5″) residual component that is obtained after subtracting the double exponential.

In the text

Fig. 4. Left panel: zoomed view (22″ × 22″/6 kpc × 6 kpc) of the Hα narrow band image from the top panel of Fig. 1 with NHI contours from the VLA B-configuration 21 cm observations (NHI = {5, 10, 20, 30, 40}×1020 cm−2, rose lines); these contours are also shown in Fig. 1 as subdued white lines. The cyclic colour map encodes Hα flux in a asinh-scale from 0 to 10−17 erg s−1 cm−2 to 10−15 erg s−1 cm−2. In the brightest central region the archival HST F550M image is inset, with an arbitrary asinh-scaling to highlight the compact super-star clusters that comprise the main stellar body of this system. The violet dashed square indicates the region that is displayed in the right panel. Right panel: inner high-SBHα region as seen in the archival HST FR656N image (13″ × 13″/3.5 kpc × 3.5 kpc) displayed here with a cyclic asinh-scale from 0 to 5 × 10−15 erg s−1 cm−2 to 1013 erg s−1 cm−2. Here the violet dashed square indicates the region that is displayed in Fig. 5.

In the text

	Fig. 5. HST archival F550M image with super-star clusters labelled according to Thompson et al. (2009); NHI contours as in Fig. 4. The viewport, indicated as a violet dashed square in Fig. 4, is 7″ × 7″ or 1.88 kpc × 1.88 kpc.
In the text

	Fig. 6. 21 cm moment-0 map of SBS 0335−052E with the colour-bar encoding the neutral column NHI in units of 1020 cm−2 using an asinh-scaling. The magenta ellipse indicates the dimension and orientation of the VLA B-configuration beam. This map has been used to draw the NHI contours in Figs. 1, 4, and 5.
In the text

	Fig. 7. Moment-0 maps from VLA B-configutration 21 cm observations. The colour mapping encodes NHI in 1020 cm−2 and the beam for each dataset is indicated as a black ellipse. Top panel: VLA-B configuration data. Bottom panel: VLA-B configuration observations, tapered at 7 kλ; this map has been used to draw the NHI contours in Fig. 8.
In the text

	Fig. 8. Hα narrowband from MUSE (log-stretch from 0 to 1.25 × 10−15 erg s−1 cm−2) inset into a grz colour composite (Pan-STARRS) with NHI contours from VLA B-configuration data tapered at 7 kλ; contours delinate NHI = {0.5(2.5σ),1, 2, 5, 8}×1020 cm−2. The beam is 15.9″ × 14.7″ oriented at 66.4° and is indicated as a hatched ellipse.
In the text

Fig. 9. Top panel: continuum subtracted [O III] λ5007 narrowband; contours indicate SB[OIII] = {5, 12.5, 25}×10−19 erg s−1 cm−2 arcsec−2. The cyclic colour map encodes flux densities from 0 to 5 × 10−19 erg s−1 cm−2 (colour bar) to 5 × 10−17 erg s−1 cm−2. The image has been smoothed with a Gaussian of 0.4″ FWHM. Bottom panel: continuum subtracted Hβ narrowband. The cyclic colour map encodes flux from 0 to 3 × 10−19 erg s−1 cm−2 to 3 × 10−17 erg s−1 cm−2. The image has been smoothed with a Gaussian of 0.95″ FWHM to enhance the visibility of the low-SB Hβ emission.

In the text

Fig. 10. Analysis of emission line ratios along the eastern (left) and western (right) filament. In the top three panels the emission in Hα, Hβ, and [O III] λ5007 is shown, now aligned such that the x-axis runs along the direction of the filament (Δl) and that the y-axis runs across the filament (Δw); the correpsonding areas have also been outlined in the bottom left panel of Fig. 1. The two bottom panels show the Hβ/Hα and [O III]/Hα ratio in bins of 1.2″ along Δl. These ratios are obtained after integrating the emission in each line along Δw. The 1-σ errors on the ratios have been computed with Eq. (2). For Hβ/Hα the Case-A recombination expectations from Storey & Hummer (1995), are indicated as dashed lines for three different temperatures (T = {1, 2, 3}×104 K in green, orange, and red); cf. Fig. 11. For Hβ/Hα the corresponding conventional Balmer decrement values are measured on the right ordinate.

In the text

	Fig. 11. Inverse Balmer decrement values for Case-A (dashed line) and Case-B (solid line) recombination as a function of temperature as calculated by Storey & Hummer (1995). The right ordinate measures the conventional Balmer decrement values.
In the text

Fig. 12. Line-of-sight velocity map (top) and velocity dispersion map (bottom) of the ionised gas; subdued contours mark various Hα surface brightness levels as in Fig. 1. As zero-point velocity v = c ⋅ z = 4053 km s−1 (Moiseev et al. 2010) was adopted. The kinematic major and minor axis after Moiseev et al. (2010) are indicated with a dash-dotted and dotted line, respectively. The viewport of this figure is 1′35″ × 1′46″ as in Fig. 1. The blue dashed square indicates the 22″ × 22″ viewport centred on the central circumgalactic medium (Fig. 2) shown in the left panel of Fig. 4.

In the text

	Fig. 13. Moment-1 map from VLA-B configuration 21 cm data showing the intensity weighted line-of-sight velocity field of H I (top: untapered data; bottom: 7 kλ tapered data). The purple contours mark various Hα surface brightness levels as in Fig. 1. For the zero-point velocity the optical recession velocity from the eastern galaxy, v0 = c ⋅ z = 4053 km s−1 (Moiseev et al. 2010), was used.
In the text

Fig. 14. Comparison between H I and H II kinematics on small scales. The H I velocity field (top left panel) is a zoomed-in view of the velocity field shown in the top panel Fig. 13. The H II velocity field (top right panel) is the same as in Fig. 12, but regions outside of detected H I have been subdued. In the bottom left panel we show the H II velocity field convolved with a 2D Gaussian (6.1″ × 4.6″, position angle −23.6°) to mimic the VLA B-configuration observations beam smearing (see text). The bottom right panel shows the smoothed Hα velocity field resampled onto the native grid of the VLA B-configuration cube.

In the text

Fig. 15. Comparison between H I and H II kinematics on larger scales; we here display a zoomed-in view of the 7 kλ tapered 21 cm data shown in the bottom panel of Fig. 13. The top panel shows the first moment from the H I data, the middle panel shows the ionised gas velocity field, with regions outside of H I detections being subdued, and the bottom panel shows the H II velocity field convolved with a 2D Gaussian (15.9″ × 14.7″, position angle 66.4°) to mimic the beam smearing and resampled onto the native grid of the H I observations.

In the text

Fig. 16. Maps showing $\mathrm{\Delta }v=v_{\mathrm{HI}}-v_{\mathrm{HII}}^{\mathrm{rep}}$, where vHI are the H I line-of-sight velocities and $v_{\mathrm{HII}}^{\mathrm{rep}}$ are the resolution matched and re-projected H II line-of-sight velocities. The left panel shows the comparison with the untapered VLA-B configuration data, i.e. the subtraction of the velocity field in the lower-right panel of Fig. 14 from the upper-left panel of Fig. 14, while the right panel shows the comparison for the 7 kλ-tapered data, i.e. the subtraction of the velocity field in the lower panel of Fig. 15 from the upper panel of Fig. 15.

In the text

Fig. 17. Toy model to explain the NW filaments as limb-brightened edges of an outflow cone. We assume a cone filled with a hot wind fluid that is surrounded by walls filled with ionised hydrogen in recombination equilibrium (T = 2 × 104 K, n = 0.05 cm−3). Details on the setup and the model parameters that are derived from the observations are given in the text. The left panel shows a projection of the resulting Hα emission (as surface-luminosity in erg s−1 kpc−1) from this model along the Y-axis of the model box; this axis is chosen to be parallel to the line of sight. The right panel shows how the Hα emission from this model would appear in our MUSE data under the assumption that the model is at the distance of SBS 03352−052E. Here we overlay some Hα contours from Fig. 1 (dotted white lines) to indicate the degree of congruence between the simulated observations of this simple toy model and reality.

In the text

Fig. A.1. Hα emission line profiles extracted along the pseudo slits shown in the bottom left panel of Fig. 1 (top panels) alongside the fitted Gaussian models (centre panels) and the residuals (data – model, bottom panels). The left panels show the results for the eastern filament, whereas the right panels show the result for the western filament and the horizontal extend is the same as in Fig. 10. The x-axis is the cross dispersion axis in units of MUSE spaxels (1 spaxel = 0.2″) and the y-axis is the dispersion axis (1 pixel = 1.25 Å). The colour bar on the left corresponds to the two top panels, whereas the colour bar on the right corresponds to the bottom panel.

In the text

Fig. B.1. Channel maps from the untapered VLA-B configuration datacube. Each channel is shown with a linear stretch from 0 to 6σ, where σ = 2.5 × 10−4 Jy beam−1 is determined from fitting the negative distribution of pixels in each channel. Contours are drawn at 2σ (dashed), 4σ, and 6σ. The bottom right panel shows the same viewport (46″/12.7 kpc along each axis), but with the Hα NB image and the total NHI contours from Fig. 1.

In the text

	Fig. B.2. Channel maps from the 7 kλ-tapered datacube, similar to Fig. B.1; here σ = 4 × 10−4 Jy beam−1. The bottom right panel shows the viewport (210″ × 114″/58.1 kpc × 31.5 kpc) over the Pan-STARRS false-colour image with the NHI contours from Fig. 8.
In the text

	Fig. C.1. Surface-luminosity projections of the cone-wall structure from Sect. 4 for various position angles ϑ and inclinations i as indicated in the panels.
In the text

	Fig. C.2. Surface-luminosity projections of the cone-wall structure from Sect. 4 for various opening angles θ as indicated in the panels.
In the text

	Fig. C.3. Simulated MUSE Hα observation of the cone-wall structure from Sect. 4 at the distance of SBS 0335-052E for fixed n = 0.05 cm−2 but varying T as indicated in the panels.
In the text

	Fig. C.4. Simulated MUSE Hα observation of the cone-wall structure from Sect. 4 at the distance of SBS 0335-052E for fixed T = 20000 K but varying n as indicated in the panels.
In the text