© The Authors 2023

1. Introduction

While models accounting for the star-formation-driven dynamical interstellar medium (ISM) successfully describe the overall picture of the various ISM components observed in nearby galaxies, details of the processes governing the distribution and physical conditions of the different phases remain an open issue. The warm ionized medium (WIM) or diffuse ionized gas (DIG) observed in the halos of star-forming disk galaxies (for a review see Haffner et al. 2009; Putman et al. 2012) poses questions regarding its origin, ionization, and heating. The existence of an extraplanar layer of low-density ionized gas was first proposed by Hoyle & Ellis (1963) to explain an absorption feature in the Galactic diffuse radio continuum background. The Galactic DIG has since been detected through the dispersion of pulsar signals (Guélin 1974; Reynolds 1991a), in emission of strong optical spectral lines in the direction of high-velocity H I clouds (Weiner & Williams 1996; Tufte et al. 1998), and from absorption lines of intermediate ions in spectra of background sources (Collins et al. 2009). The DIG in external galaxies was originally studied using Hα imaging (Dettmar 1990; Rand et al. 1990) or long-slit spectroscopy (Dettmar & Schulz 1992), and more recently using integral field spectroscopy (Ho et al. 2016; Jones et al. 2017; López-Cobá et al. 2017; Levy et al. 2019; Rautio et al. 2022).

The DIG in the Milky Way occupies the temperature range between 6000 K and 10 000 K and densities between 0.03 cm⁻³ and 0.08 cm⁻³ with a filling factor around f ≈ 0.4–0.2 in the volume 2–3 kpc above the galactic midplane (Reynolds 1991b; Hill et al. 2008; Haffner et al. 2009). For the external galaxy NGC 891 Rand et al. (1990) arrived at a similar density estimate. More commonly, however, the DIG is characterized in terms of the behavior of its optical emission lines rather than its thermodynamical conditions. Generally, the DIG displays elevated levels of the line ratios [S II]λ6717/Hα and [N II]λ6583/Hα and lower [O III]λ5007/Hα values than classical planar H II regions (e.g., Mathis 1986; Reynolds et al. 1999). Additionally, [S II]λ6717/Hα and [N II]λ6583/Hα tend to rise with increasing altitude above the midplane (e.g., Golla et al. 1996; Levy et al. 2019). These varying line ratios are contrasted by a nearly constant [S II]λ6717/[N II]λ6583 ratio (e.g., Haffner et al. 1999; Reynolds et al. 1999). While the abnormal ratios of forbidden lines to Hα relative to H II regions can be explained by photoionization at low ionization parameters, q (Domgorgen & Mathis 1994), the constancy of [S II]λ6717/[N II]λ6583 requires an additional source of heating (Reynolds & Cox 1992; Rand 1998; Reynolds et al. 1999; Collins & Rand 2001). Wood & Mathis (2004) successfully modeled this additional heating source with the spectral hardening of the ionizing photon field from O- and B-type stars as it propagates from the midplane into the halo. These general properties of the optical DIG emission lines are commonly observed, but not ubiquitous. Golla et al. (1996), for example, observed a rising trend for [S II]λ6717/[N II]λ6583 with increasing altitude in NGC 4631 rather than the usual constancy, but at a lower rate than predicted by the Domgorgen & Mathis (1994) models. We note, however, that the interaction between NGC 4631, NGC 4656, and NGC 4627 may introduce complications that are not considered in general models of the DIG.

The necessary energy budget for the production and maintenance of the DIG implies that the radiation field from young massive stars in the disk is the primary ionizing source (Reynolds 1984; Mathis 1986; Domgorgen & Mathis 1994; Haffner et al. 2009). Additional sources that have been proposed to contribute at a significant level include, among others, shocks (Martin 1997), radiation from hot-gas-cool-gas interfaces (Haffner et al. 2009), and hot low-mass evolved stars (HOLMES; Sokolowski & Bland-Hawthorn 1991; Flores-Fajardo et al. 2011).

The extraplanar Hα emission in DIG-halos is found to correlate with tracers of star formation in the disk (Dettmar 1992; Rossa & Dettmar 2003a,b; Jones et al. 2017; Rautio et al. 2022) on local and galactic scales. This coupling to the star formation indicates that the ionized gas originates in the disk and is transported into the lower halo by feedback processes. The necessary gas transport process in this paradigm connects to the study of galactic winds, which are defined as gas flows with a velocity component that exceeds the escape velocity of the galactic potential, and galactic outflows, which remain bound to their host galaxy. Many questions regarding the detailed workings of galactic outflows and winds remain unsolved (see reviews by Veilleux et al. 2005; Rupke 2018; Zhang 2018), but it is clear that they are important mechanisms regulating the evolution of galaxies, the circumgalactic medium, and the intergalactic medium (e.g., Heckman et al. 1990), and consequently also the disk-halo interface. In the context of this paper, the questions about the ionization mechanism(s), the structure of the base of the wind, and the kinematics of the flow are especially noteworthy. The ionization of the gas in diffuse ionized halos and in galactic winds is a longstanding, not fully resolved question, as discussed above; photoionization from OB stars is clearly an important contributor (e.g., Heckman et al. 1990), but shocks can also play a role (e.g., Rich et al. 2010), as well as ionization by photons from the hot gas inside the outflow or wind (e.g., Sarkar et al. 2022). Even non-equilibrium effects can play a role (Gray et al. 2019a,b; Sarkar et al. 2022). To understand the structure of the outflow or wind and its energetics, knowing the main driving mechanism is critical; current models include contributions from overpressure due to the hot gas and its mechanical feedback on the surrounding gas, radiation pressure, and the interaction of cosmic rays with the magnetic field frozen within the ionized gas (e.g., Fielding & Bryan 2022). This short introduction to outflows and winds is focused on the warm, ionized outflow component that we discuss in this paper but it is important to keep in mind that this is only one part of a multiphase phenomenon, with hot gas inside the core of the outflow, cool gas at the compressed shell, clouds inside the flow, and cosmic rays interacting with the magnetic field frozen into the ionized gas.

2. Data

We used archival Integral Field Unit (IFU) data observed with the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) originally intended by Comerón et al. (2019) to investigate the stellar kinematics in the thick disks of star-forming edge-on galaxies. Rautio et al. (2022) used this same data set to analyze the DIG of the majority of galaxies in this sample. In this work we picked IC 1553 from this sample to study the DIG of a single galaxy in more detail.

The MUSE spectrograph has a channel width of 1.25 Å spanning the spectral range of 4750 − 9351 Å, and spaxels of size 0.2″ × 0.2″ over a total 1′×1′ field of view. The data cube was pre-reduced with the MUSE pipeline (Weilbacher et al. 2012) and ZAP (Soto et al. 2016). For a detailed description of the reduction process see Comerón et al. (2019).

We used the fitting code 3-Dimensional Emission Line Analysis (3DELA, Schilling & Bomans, in prep.) to extract the emission line properties. 3DELA is a collection of Python scripts designed for the fast evaluation of IFU data cubes. The treatment of spectra with a small signal-to-noise ratio has been of special importance in detecting features of the gas dynamics far away from the galaxy disk, such as the presence of outflows. The spectral line fitting requires an initial estimate of the systemic velocity, preferably obtained by measuring the wavelength of the shifted Hα line in the spectral cube in a region near the suspected center of mass. The spectral range for the fit of the line under consideration is known therefrom.

Primarily the code has been developed for the study of the galaxies of the Calar Alto Legacy Integral Field Area Survey (CALIFA; Sánchez et al. 2012; García-Benito et al. 2015). Due to the survey’s low spectral resolution of 6 Å in the V500 setup, the spectral lines were only defined by a few data points. For this reason 3DELA fits lines with a single Gaussian component. The current version of 3DELA allows the fitting of emission lines from data cubes from other instruments – like MUSE – and surveys (e.g., Sloan Digital Sky Survey (SDSS; Jönsson et al. 2020) or Mapping Nearby Galaxies at APO (MaNGA, Bundy et al. 2015)) by transposing them to an internal standard format conforming to FITS. This procedure requires some value-substitution strategies due to missing values in the FITS headers. Typical tasks of preparation or post-processing (e.g., estimation of the systemic velocity, dereddening, the plotting of specific results or the calculation of diagnostic diagrams) are performed by discrete scripts using the internal standard format of the data cubes.

The data cube was rebinned to increase the signal-to-noise ratio with a regular grid of 50 × 50 resolution elements over the field of view (i.e., 6 × 6 spaxels per resolution element), resulting in an effective spatial resolution of 1.2″. We chose a regular grid over other commonly used practices like the Voronoi tessellations (e.g., Cappellari & Copin 2003) to simplify the decomposition of a bin’s position into its projected radial distance from the center of IC 1553 and its height above the disk, and to be able to sample this coordinate space in a consistent pattern. In order to have the spaxels aligned with the major and minor axes of the galaxy, the orientation of the grid was determined in the following manner. We extracted a synthetic V-band image from the data cube and separated it into ten bins along a rough estimate of the major axis. The midplane in each bin is then taken to be the position that ensures the V-band flux is equally distributed above and below the plane of the disk. Finally, these positions were fit with a linear function to determine the plane of the disk. The nominal systemic velocity of IC 1533 was ascertained by fitting the Hα profile in the spatially integrated spectrum of the disk of IC 1553 (masking all bins of altitude |z|> 10″) with the busy function (Westmeier et al. 2014a,b). This fit yields a centroid wavelength for each of the approaching and receding halves of the galaxy (see Fig. 1), which can then be averaged to calculate the systemic velocity. In what follows we take the center of IC 1553 to be the position along the midplane of the disk, where the Hα emission line is Doppler-shifted with this systemic velocity. We assume an uncertainty in the determination of the systemic velocity equivalent to half the width of the spectral channels, in other words ±0.625 Å, and obtain a systemic velocity value of v_sys = (2914 ± 29) km s⁻¹. The results of the rebinning and the fits to the disk plane and the center of the galaxy are shown for the synthetic V-band image in Fig. 2. Positive values for the projected radial distance, d, signify a position north of the center along the dashed line in Fig. 2. Positive values for the height above the midplane, z, signify a position east of the disk. We note that the emitting extended object in the southeastern corner of the field of view in Fig. 2 is a background galaxy and not interacting with IC 1553. This background galaxy contaminates the Hα line of IC 1553 (see Appendix A); we therefore masked this area in all diagnostic maps using Hα.

Fig. 1. Visualization of the busy function fit to determine the systemic velocity of IC 1553. The integrated Hα profile was obtained by collapsing the data cube along both spatial axes for bins with altitude |z|≤10″. The central wavelength of the profile is taken as the mean of the centroid wavelengths of the red- and blue-shifted components in the busy function fit. We estimate the uncertainty in this determination as ±0.5 channel width.

Fig. 2. Synthetic V-band images constructed from the data cube at the original resolution (left), and after rebinning with a regular 50 × 50 grid of resolution elements (right). The dashed line and the white circle are the estimates for the plane of the disk and the dynamic center of IC 1553, respectively.

In some quantitative considerations in the following sections we employ data quality cuts on the Hα flux. In these considerations the detection limit is given as a multiple of the background level, σ ≈ 5.13 × 10⁻⁸ ergs s⁻¹ cm⁻² sr⁻¹ ≡ 0.21 Rayleighs (R, 1 R = 10⁶/(4π) photons s⁻¹ cm⁻² sr⁻¹), which we estimate with the standard deviation of non-zero Hα flux values in bins far away from the disk of IC 1553 (25″ ≤ |z|≤30″).

Within this framework we have produced maps for the emission line fluxes commonly used to characterize the warm ionized medium – including Hβ, [O III]λ4959, 5007, [O I]λ6300, Hα, [N II]λ6548, 6583, [S II]λ6716, 6731, and [S III]λ9069 – as well as their observed velocity offsets from the rest-frame system and their full width at half maximum (FWHM).

The emission line fluxes were corrected for dust attenuation using the same methods as Domínguez et al. (2013). Assuming typical gas conditions in star-forming regions for the electron temperature, T_e = 10⁴ K, and density, n_e = 100 cm⁻³, yields a largely constant Balmer decrement (Hα/Hβ)_int = 2.86 for Case B recombination (Osterbrock & Ferland 2006). The extinction in magnitudes at wavelength λ is then given from the measured Balmer decrement as (Domínguez et al. 2013):

(1)

where k(λ) can be sampled from the reddening curve of Calzetti et al. (2000). The observed Hα and Hβ fluxes contain not only the emission from the ionized gas, but also some component from the emission of the stellar population along the line of sight, as do all recombination lines. Typically, the stellar population produces a strong absorption line feature in Hβ and a weaker absorption line feature in Hα. This alters the measured Balmer decrement, and constitutes a form of contamination in the context of Eq. (1). We corrected the observed Balmer decrement for this effect by separating the emission components from the stellar population and the ionized gas using the Penalized PiXel-Fitting code (pPXF; Cappellari & Emsellem 2004; Cappellari 2017) with the E-MILES stellar population synthesis models (Vazdekis et al. 2016). The left panel of Fig. 3 shows a map of the derived gas-phase Balmer decrements.

Fig. 3. Maps of the Balmer emission in IC 1553. Left: map of the Balmer decrement Hα/Hβ obtained from the spectral template fitting with pPXF. Middle: map of the extinction-corrected Hα emission line flux in Rayleighs with linear scaling. Right: the same as the middle panel, but with logarithmic scaling. The contours are arbitrary equidistant levels from the synthetic V-band image meant to give orientation for the eye. The white square in the lower left corner corresponds to the masked background galaxy. The red circles in the middle panel mark the H II regions used for the oxygen abundance measurement listed in Table 1. The boxed region marked A in the right-hand panel is identical to the gray shaded area in Fig. 5 below.

The FWHM was corrected for widening from the line spread function (LSF) via

(2)

The parametrized LSF given by Bacon et al. (2017) for the udf10 field (their Eq. (8)) was used in this correction.

We assume the same distance to IC 1553 as Comerón et al. (2019), namely 36.5 Mpc (Tully et al. 2016). Some of the general properties of IC 1553 are summarized in Table 1.

3. Results

3.1. Hα distribution

The strongest tracer of the ionized component of the warm ISM is the Hα emission line. The middle and right panels of Fig. 3 show maps of the extinction-corrected Hα fluxes on a linear and logarithmic scale under contour lines derived from the V-band continuum image. From the linear scale map it is evident that IC 1553 exhibits a pronounced asymmetry in Hα brightness between its northern and southern sides, with the majority of the flux being produced in the south. This indicates that the southern half is currently undergoing significantly more star formation than the north.

In Fig. 4 we show maps of alternative tracers for star formation as contour lines. The first panel shows the brightness in the W4-band at 22 μm from the ALLWISE Data Release (Cutri et al. 2021) of the Wide-field Infrared Survey Explorer (WISE), which traces star formation via the emission of dust heated by stellar ultraviolet light (Rieke et al. 2009; Kennicutt & Evans 2012; Cluver et al. 2017). The second panel of Fig. 4 shows contours from the first data release of the Rapid ASKAP Continuum Survey (RACS; McConnell et al. 2020) at 887.5 MHz. Radio continuum emission is produced by synchrotron radiation from relativistic electrons, which are accelerated either by an active galactic nucleus (AGN) or by supernovae (Condon 1992; Vaddi et al. 2016). The fact that IC 1553 does not show an AGN signature in optical emission line ratios (see Sect. 3.2 below and Rautio et al. 2022) suggests that here the radio continuum emission is related to star formation. Investigations of global star formation rates and radio continuum luminosities show a linear correlation, but usually use radio data at 1.4 GHz due to the strong far-infrared/radio correlation (Murphy et al. 2011; Kennicutt & Evans 2012). Examinations of local star formation rates and radio luminosities down to 150 MHz even show a superlinear relation (Heesen et al. 2019, 2022), so that we can use the RACS image for a relative comparison of the star formation in the northern and southern halves of IC 1553. The nominal angular resolution is 12″ for WISE W4 images (Wright et al. 2010) and 15″ for the RACS-low band, so both are able to separate the northern and southern halves of the galaxy’s ∼1′-wide disk from each other.

Fig. 4. Synthetic V-band image of IC 1553 overlaid with equidistant contours from the W4-band image from the ALLWISE Data Release at 22 μm (left), from the RACS-low image at 887.5 MHz (middle), and from the Hα image extracted from our data cube smoothed with a Gaussian kernel (right). The nominal angular resolutions are 12.0″ for the WISE W4-band, 15.0″ for RACS-low, and the Hα image was convolved with a 12.0″-wide Gaussian for comparison.

As can be seen in Fig. 4, in both the W4 and RACS-low bands the galaxy appears reasonably symmetrical across the extent of its stellar disk. For comparison, the right-hand panel in Fig. 4 shows the Hα image extracted from the data cube convolved with a 12.0″-wide Gaussian kernel to emulate a similar spatial resolution. The primary difference between these tracers is that Hα is sensitive on a shorter timescale of ≲10 Myr while infrared and radio continuum emission are sensitive to star formation for up to 100 Myr (Kennicutt & Evans 2012). This implies that the northern half of IC 1553 underwent comparable star formation rates to what its southern half currently exhibits, but was quenched sometime within the last 10–100 Myr.

From the right-hand panel of Fig. 3 it can be seen that diffuse Hα emission is detected up to altitudes of ∼20″ above the midplane. The north-south asymmetry is also present in this more diffuse layer of ionized gas, where the emission above the southern star-forming cluster is generally stronger than above more quiescent regions of the disk.

We use this extended emission to characterize the spatial distribution of the DIG. The DIG of edge-on galaxies is usually well fit in the vertical direction by a two-component exponential function (Reynolds 1991a), or a single exponential component, if the plane of the disk is masked sufficiently (e.g., Levy et al. 2019). Rautio et al. (2022) obtained a two-component fit for IC 1553 averaged over its projected radial extent and characterized its radial variation in terms of the total Hα flux. We aim to model the radial DIG distribution by fitting exponential functions at all projected radial positions according to our binning pattern and using all bins with Hα fluxes above 5σ = 1.06 R.

First, we attempted a two-component fit to the measured Hα flux I_Hα along the vertical axis, z, at every radial position, d, of the form

(3)

where the first component characterizes the H II gas contained in the disk and the second component characterizes the extraplanar diffuse ionized gas (eDIG), but we found that the data are not sufficient to constrain the free parameters in this approach. We therefore fit the two components individually by treating bins in the plane of the disk separately from bins in the halo:

(4)

In this approach the height, z_sep, where the two components are separated should be large enough to suppress a bias in the eDIG parameters from the disk emission and small enough for the eDIG component to contain sufficient data points for the fit. In addition to fitting the linear Hα flux according to Eq. (4) we have used logarithmic weighting to give more power to the low-brightness bins at high altitudes. Initially, we separated the disk and eDIG components at z_sep = 5″, but found that the logarithmic weighting resulted in systematically larger eDIG scale heights than those resulting from the linear weighting at all radial positions, which we interpreted as contamination from the disk component. We, therefore, increased z_sep to 8″ in this analysis, at which point both weighting schemes produce consistent eDIG scale heights throughout the northern half of IC 1553 (for a discussion about the southern half, see below). The left panel of Fig. 5 shows an example of the results from this fitting process at d = −15.0″ for both components and both weighting schemes. Appendix C contains analogous plots for the fitting at every projected radial distance.

Fig. 5. Visualization and results of the Hα scale height fitting. Left: visualization of the IHα fit at the example projected radial distance of d = −15.0″. The dotted lines mark the separation height, zsep = 8″, between bins used for fitting the disk H II (red curves) and eDIG (blue curves) components. The solid curves are derived with linear weighting of IHα according to Eq. (4) and the dashed curves result from logarithmic weighting of IHα during the fit. Right: radial distribution of scale heights for the eDIG (blue squares from logarithmic weighting, green circles from linear weighting) and disk H II (red crosses from logarithmic weighting, orange crosses from linear weighting) components obtained from the fitting. The shaded areas around the scale height curves give the 1σ errors obtained from the corresponding estimated parameter covariance. The dashed lines are best fits of Eq. (5) to the eDIG component scale heights within |d|≤20″. The gray shaded area, A, marks the region where the logarithmic weighting produces significantly larger eDIG scale heights than the linear weighting.

The right panel of Fig. 5 presents the radial behavior of all scale heights we derived in this way for IC 1553. For the galaxy’s northern half we obtain eDIG scale heights only out to d = 16.2″ because the Hα flux falls below the 5σ quality cut at z_sep beyond that point. Within the sampled range in the north it can be seen that the Hα scale height declines linearly with the radial distance from the center on average from h_eDIG ≈ 1 kpc to 0.6 kpc, and that both weighting schemes yield consistent results. For the linear weighting scheme this trend is mirrored in the southern half down to d ≲ −20″, leading to a behavior of the form

(5)

with the parameter values c₁ = 0.949 kpc, c₂ = 0.118, and d₀ = −0.002 kpc. The logarithmic weighting, however, produces significantly larger h_eDIG at −20″ ≤ d ≤ −5″ (the shaded gray area in Fig. 5 and boxed area in the last panel of Fig. 3, both labeled as A) than the linear weighting, leading to best-fitting parameters for Eq. (5) of c₁ = 1.087 kpc, c₂ = 0.116, and d₀ = −0.880 kpc. At the southern edge of the disk (d ≤ 20″) h_eDIG appears to increase with the radial distance, but the fit is not entirely reliable here, because of several problematic bins in the Hα map and the smaller sampled z-range resulting from the masking of the background galaxy (compare with the right panel of Fig. 3). We interpret and discuss this result in Sect. 4.1 below.

3.2. Ionization structure of the DIG

Figure 6 shows maps of the emission line ratios commonly used to derive information about the ionization conditions of interstellar gas. These maps are consistent with some of the observed properties of the DIG in other galaxies and the Reynolds-layer of the Milky Way: (1) The line ratios [N II]λ6583/Hα and [S II]λ6717,6731/Hα generally rise with increasing height above the disk (e.g., Golla et al. 1996; Haffner et al. 1999; Collins & Rand 2001; Levy et al. 2019) and (2) [S II]λ6717,6731/Hα and [N II]λ6583/Hα reach elevated levels compared to classical H II regions (e.g., Reynolds 1985, 1988; Golla et al. 1996; Haffner et al. 1999; Levy et al. 2019). For IC 1553, however, the magnitude of these effects varies drastically over the projected radial position. [N II]λ6583/Hα has a clear conical structure of increased values centered slightly north (∼5″) of the galactic center. The same can be seen for [O I]λ6300/Hα, though here the structure is even further off-centered (∼10″) and less pronounced. [S II]λ6717,6731/Hα displays a similar structure of increased value, producing a clear asymmetry between the northern and southern halves of IC 1553. [O III]λ5007/Hβ on the other hand peaks in the plane of the disk coincident with the prominent star-forming regions in the south (compare with Fig. 3), but everywhere else its value is larger in the eDIG than in the planar gas.

Fig. 6. Maps of the emission line ratios commonly used in the classical diagnostic diagrams. First row: [N II]λ6583/Hα (left), [S II]λ6717,6731/Hα (middle), and [S II]λ6717/[N II]λ6583 (right), which is commonly analyzed in studies of the DIG. Second row: [O I]λ6300/Hα (left) and [O III]λ5007/Hβ (middle), as well as the η′ parameter (right) as described in Sect. 3.2 and defined in Appendix B, which quantifies the distance to the Kewley et al. (2001) starburst limit in the [O III]/Hβ vs. [S II]/Hα diagnostic diagram. The black demarcation lines in the η′ map separate bins below the starburst limit (η′≤0.5) from bins above that limit (η′> 0.5). The gray and white contours are at the same equidistant arbitrary levels from the synthetic V-band image as in Fig. 3.

One unusual feature of the DIG in IC 1553 is the map of the [S II]λ6717/[N II]λ6583 ratio. Multiple studies found this ratio to remain largely constant with z (e.g., Golla et al. 1996; Otte & Dettmar 1999; Haffner et al. 1999; Rand 2000). Here, however, this is only observed in the volume coincident with the conus of increased [N II]λ6583/Hα. Outside of this volume, IC 1553 displays a rising trend of [S II]λ6583/[N II]λ6583 with increasing height.

In Fig. 7 we show the classical ionization diagnostic diagrams (Baldwin et al. 1981; Veilleux & Osterbrock 1987) with the limiting surface for starburst galaxies from Kewley et al. (2001) and the best-fitting photoionization model grid from Levesque et al. (2010) for comparison. Figure 7 (as well as Figs. 8 and 9) contains only data points where Hα was detected at a 5σ-level (1.06 R). The population of bins most closely resembles the models for a current instantaneous star formation burst of age 0 Myr with an electron density of n_e = 100 cm⁻³. The grid covers a range of metallicities, Z, between 0.001 and 0.04 where Z_⊙ ≡ 0.02, and ionization parameters, q, between 1 × 10⁷ cm s⁻¹ and 4 × 10⁸ cm s⁻¹. The photoionization models can reproduce the observed population in the log([O III]/Hβ) vs. log([N II]/Hα) diagram quite nicely (left panel of Fig. 7) for roughly solar metallicity (Z ≈ 0.02) and an ionization parameter around 2 − 4 × 10⁷ cm s⁻¹. Roughly half of the bins fall slightly outside of the volume covered by the grid and closer to the Kewley et al. (2001) starburst limit. In principle this could be explained by spectral hardening of the photons during their propagation from the disk to the halo (Wood & Mathis 2004). However, the other two diagnostic diagrams display a clear excess of [S II] and [O I] emission compared to the models, with a significant fraction of the population actually crossing the starburst limit.

Fig. 7. Typical emission line diagnostic diagrams for the analysis of ionization sources. The data points are the population of bins with Hα flux ≥5σ, colored by their η′ value. The black dotted line marks the theoretical limit for starburst galaxies from Kewley et al. (2001). The solid and dashed lines mark the metallicity (Z) and ionization parameter (q) grid of photoionization models described in Levesque et al. (2010) for the best-fitting star formation age of 0 Myr and electron density of 100 cm−3.

Since the behavior of the emission line ratios varies strongly over the projected radial position, we are interested in mapping the bins in the plane of the sky to their distance from the Kewley et al. (2001) starburst limit in the diagnostic diagrams to gauge the spatial distribution of the ionization conditions more clearly. Rautio et al. (2022) used the η-parameter from Erroz-Ferrer et al. (2019) for this purpose, which quantifies the position in the [O III]/Hβ vs. [N II]/Hα diagram. For IC 1553 there is merit in using the [S II]/Hα or [O I]/Hα diagnostic instead, because there the data deviate much more strongly from the photoionization models. We therefore define the parameter η′ analogously to η so that η′< 0.5 marks a position to the left of the Kewley et al. (2001) starburst limit in the [O III]/Hβ vs. [S II]/Hα diagram, and η′> 0.5 to its right. The specific calculation of η′ is described in Appendix B. The bottom-right panel of Fig. 6 shows the map of η′ for IC 1553 and visualizes that most of the volume occupied by the sampled eDIG produces emission line ratios beyond the starburst limit. Only the innermost halo above the star-forming cluster in the southern disk exhibits bins with η′< 0.5.

The Wood & Mathis (2004) radiative transfer models can reproduce the observed [S II]/Hα and [O I]/Hα values in principle, but not without also producing an elevated [N II]/Hα. They can therefore explain the ionization behavior in the eDIG above the center of the galaxy, but not above the outskirts of the disk, where [S II]/[N II] rises with increasing height.

This deficit in [N II] emission compared to [S II] and [O I] in conjunction with the moderate ratio values of [O III]/Hβ observed in the data leads us to consider shocks as an additional source of ionization for the DIG of IC 1553. To this end, we compare the observed emission line ratio to the library of shock models from Allen et al. (2008). The best agreement with our data is achieved under the conditions of solar abundances with dust depletion, as described by Dopita et al. (2005), and a preshock density of n = 1.0 cm⁻³. These models are plotted for a range of shock velocities of v = 100 − 300 km s⁻¹ and magnetic fields of B = 10⁻⁴ − 10 μG in Fig. 8.

Fig. 8. Same diagnostic diagrams as shown in Fig. 7 (dotted black and red lines are the starburst limit and photoionization models described therein), now including the velocity (v) and magnetic field (B) grid (dashed and solid lines) for the best-fitting shock models from Allen et al. (2008), in which the preshock density n = 1 cm−3 and the chemical abundances described in Dopita et al. (2005) were used.

The line ratios for the majority of bins fall somewhere between the parameter space covered by the Levesque et al. (2010) photoionization models, and the curve given by the Allen et al. (2008) shock ionization models for the shock velocity of v = 225 km s⁻¹. This distribution indicates that both photo- and shock ionization are relevant for the ionization state of the eDIG in IC 1553, and that the emission lines in each spatial bin may contain components from both processes with differing relative weights. If we want to characterize properties of the shock, we should consider the edge of the distribution, which is furthest away from the photoionization model grid, because there the contribution from the shocked region to the total emission is maximal. By this reasoning, the data points in Fig. 8 are consistent with a mixture of emission produced by photoionization on the one hand and a shock with the approximate velocity of 225 km s⁻¹ on the other hand. The highest observed ratios for [N II]/Hα and [O I]/Hα cannot be reproduced in detail by the models, but it is conceivable that a more fine-tuned model, specifically in terms of a lower gas density or slightly adapted abundances, might be able to.

If shocks make a significant contribution to the ionization of the DIG we expect to see a positive correlation between observed line widths and shock-sensitive line ratios (Kewley et al. 2019). We therefore plot the corresponding ratio values over the Hα FWHM in Fig. 9 and indeed such a correlation is present for all ratios with a Spearman correlation coefficient r_s of ≈0.64 and p-values ≪0.1. For [N II]/Hα the correlation is only clear for bins with η′≳0.5, in other words for lines of sight in the halo sufficiently far away from the southern star-forming cluster (compare with Fig. 6). Figure 9 also shows that the Hα linewidth correlates with η′, which indicates that shocks play a significant role in shaping the ionization structure of the disk-halo interface of IC 1553.

Fig. 9. Shock-sensitive emission line ratios [N II]λ6583/Hα, [S II]λ6717, 6731/Hα, and [O I]λ6300/Hα vs. the Hα FWHM of all bins with Hα flux ≥5σ, colored by their η′ value.

3.3. DIG kinematics

The top panels of Fig. 10 show maps of the line-of-sight component of the velocity, and the FWHM as a proxy for the velocity dispersion, measured on the Hα emission line. The line-of-sight velocity map is consistent with the ones published by Comerón et al. (2019) and Rautio et al. (2022), displaying observed rotation velocities up to v_c ∼ 140 km s⁻¹ as well as a northward bend in lines of constant v_c. In the velocity map as well as the Hα flux map (right-hand panel in Fig. 3) we also see a slight warp of the star-forming disk at the northernmost edge of the FOV, but not in the continuum contours. The NASA/IPAC Extragalactic Database lists no potential companion in a 1 deg radius, making recent interaction an unlikely explanation for this warp.

Fig. 10. Maps of the emission line kinematics of the DIG in IC 1553. Top left: map of the Hα line-of-sight velocity in km s−1 with contours from the continuum map. Top right: map of the Hα FWHM in km s−1 after correcting for LSF broadening with Eq. (2). Bottom row: FWHM maps measured on [S II]λ6717 (left) and [N II]λ6583 (right), LSF-corrected in the same way. The white contours are from the continuum map as before. The violet contours visualize the alignment with the respective emission line flux by separating bins based on the following flux thresholds relative to Hα: [S II]λ6717,6731/Hα ≷ 1.0 (left) and [N II]λ6583/Hα ≷ 0.60 (right).

The FWHM map shows a fairly homogenous velocity dispersion across the disk of FWHM = 60 − 70 km s⁻¹, with one exception: in the northern half at d = 8″ there appears to be a channel with FWHM ≈ 90 km s⁻¹ in Hα connecting roughly conical volumes in the halo with significantly increased dispersion of FWHM ≥ 120 km s⁻¹. This picture strongly suggests the presence of an outflow, though the spectral resolution of ≳55 km s⁻¹ is not sufficient for a separation of kinematic components.

Outside of the volume of this suspected outflow, the velocity dispersion in the eDIG reaches elevated levels of FWHM ∼ 100 − 110 km s⁻¹ if it is measured on Hα, but only up to ∼80 km s⁻¹ if measured on [S II]λ6717 and [N II]λ6583, which are shown in the bottom panels of Fig. 10. Aside from this difference, the linewidths of [S II]λ6717 and [N II]λ6583 are generally morphologically consistent with Hα. The most noticeable difference is that [S II]λ6717 exhibits a smaller FWHM in the region of the suspected outflow. Among the bins with I_Hα ≥ 5σ the FWHM peaks at the position (d, z) = (9.0″, −13.9″) with the maximum FWHM = 165 km s⁻¹ in [S II] and FWHM = 178 km s⁻¹ in [N II] and Hα. The maximum FWHM approaches the expected shock velocity of ∼225 km s⁻¹ inferred from the ionization diagnostic diagrams (Fig. 8).

The bottom panels in Fig. 10 also show a violet contour line that encompasses bins of high [S II]/Hα and [N II]/Hα values, respectively, to compare the alignment of the line width- and flux-producing mechanisms. The corresponding thresholds for these contour lines are [S II]λ6717, 6731/Hα ≷ 1.0 and [N II]λ6583/Hα ≷ 0.6. While for [S II] the volume of highest relative flux coincides with the volume of highest FWHM, there appears to be a misalignment for [N II], where the volume of highest relative flux is centered a few arcseconds further south. This misalignment indicates that the [N II] emission traces a composite of photo- and shock ionization – both mechanisms produce similar quantities of [N II]/Hα in the models shown in Fig. 8 – while [S II] primarily traces shock ionization in the DIG of IC 1553.

Assuming that an off-centered outflow is responsible for some of the discussed kinematic and ionization properties of IC 1553, we are interested in whether or not the escape velocity is reached; depending on the answer, the flow of matter can enrich primarily either the intergalactic medium or the eDIG. We therefore modeled the galactic gravitational potential with galpy (Bovy 2015) under the condition that it reproduces the rotation curve published by Comerón et al. (2019). We have confirmed that the data processing in this study yields a consistent rotation curve (see Fig. C.1). The galactic potential is modeled with three components:

– A Miyamoto-Nagai (MN) potential component (Miyamoto & Nagai 1975) for the disk with a minor-to-major axis ratio of b/a = 0.1. For comparison, Comerón et al. (2018) measured radial scale lengths for the disk of IC 1553 between 12.8″ and 14.6″ and vertical scale heights of 0.9″ and 2.7″ for the thin and thick disk, respectively.

– A Hernquist (H) potential component (Hernquist 1990) for the bulge.

– A spherically symmetric Navarro-Frenk-White (NFW) potential component (Navarro et al. 1996) for the halo.

Under these conditions, each component requires one parameter describing the scale radius of the potential (a_disk, a_bulge, a_halo) and one parameter for the amplitude given by the mass of the component. We describe the individual amplitudes in terms of relative weights (α, β, γ), so that the superposed potential is given as

(6)

where A gives the absolute mass scale of the amplitude.

We constrain the seven free parameters with a grid search over the parameter grid given in Table 2. The scanned ranges of A and β are informed by Comerón et al. (2018) giving a total mass of the luminous matter in IC 1553 of around 15 − 16 × 10⁹ M_⊙, and a mass ratio between the bulge (or central matter concentration) and disk of around ∼0.05. The other parameter ranges and increments were chosen to facilitate reasonable computation times. The grid search was conducted in two steps. First, we produced potentials for all combinations of the six shape-determining parameters (all except the amplitude A) and normalized the resulting rotation curve so that it reaches a value of v = 1 at the end of the sampled radius interval. The similarity to the rotation curve model of Comerón et al. (2019) sampled to the observed extent and normalized in the same way was then quantified with the R² score function for regression. For the second step, all potential models with scores of R² < 0.99 were discarded. For each remaining model, we obtained the best-fitting amplitude from the physical velocity scale of the rotation curve, again using the R² score and discarding parameter combinations with R² < 0.99, leaving 612 499 combinations for further evaluation. While these models reproduce the observed rotation curve quite well, they suggest an unusually high mass fraction of the disk component of more than 50%. This property and its influence on the escape velocity estimate is further discussed in Sect. 4.3.

Since galpy supports the evaluation of escape velocities only in the midplane in its current version (v1.7), we obtained the escape velocity at the position x₂ = (9.0″, −13.9″), where the highest FWHM is observed in the data, from the relation

(7)

using the known quantities at the corresponding position x₁ = (9.0″, 0″) within the midplane. Assuming the projected radial distance of these positions is subject to projection effects, the actual radial distance will be r = 9.0″/cos φ for a projection angle, φ, between the radial vector and the plane of the sky. v_esc(x₂) was calculated for multiple projection angles up to φ = 75°, at which point the deprojected radial distance approaches the radial extent of the disk of IC 1553. Finally, the resulting v_esc populations were translated into probability density distributions using a Gaussian kernel density estimate; these are shown in Fig. 11.

Fig. 11. Gaussian kernel density estimates of the escape velocity, vesc, at the position of the spaxel with the highest Hα FWHM derived from fitting gravitational potentials to the rotation curve of IC 1553. vesc was modeled for multiple projection angles, φ, between the plane of the sky and the deprojected positional vector from the galactic center to this spaxel.

The models yield lower limits on the escape velocity in the range 190 ≤ v_esc/km s⁻¹ ≤ 225 depending on the projection angle, φ, which is exactly the velocity range between the maximum measured FWHM in the region of the outflow and the shock velocity inferred from the best-fitting shock models in the ionization diagnostic diagrams (Fig. 8). The probability density distribution of v_esc only peaks below 225 km s⁻¹ for high projection angles φ > 45° and extends up to 275 km s⁻¹. Assuming that the highest velocity is reached within the sampled volume and neglecting drag forces from the high-altitude halo gas, it is therefore unclear whether the bulk of the warm ionized phase of the supposed outflow reaches the escape velocity or not, according to our gravitational potential models. This result is generally consistent with other studies concerning kinematics of ionized gas produced by stellar feedback (e.g., Martin 1998; Bomans & Weis 2014).

4. Discussion

4.1. eDIG scale heights

In Sect. 3.1 we derived Hα scale heights as a function of the projected radial position for two single-component exponential models representing the eDIG and the ionized ISM confined to the disk. We fit the model parameters to both the I_H II and log(I_H II) profiles to give more weight to the emission from low- and high-altitude gas respectively.

The resulting eDIG scale heights range from 2.9″ ≡ 0.51 kpc to 7.0″ ≡ 1.24 kpc within the central 20″ of projected radial distance to the dynamical center of IC 1553, which is generally comparable to eDIG scale heights obtained in other studies (e.g., Levy et al. 2019). The southern tip of the disk at d ≤ −20″ might host larger scale heights but the fit in this region is not reliable. Rautio et al. (2022) report the averaged scale height h_eDIG = 1.39 kpc for IC 1553, which falls slightly above the highest values we find. Rautio et al. (2022) obtained this value from narrowband Hα imaging extending to an altitude of 6 kpc, while we sampled the MUSE data out to ∼3.5 kpc. The biggest methodological difference is that we fit single-component exponential functions to the eDIG Hα flux at |z|≥8″ ≡ 1.4 kpc whereas Rautio et al. (2022) use a two-component model for data at |z|≥0.25 kpc. It is therefore possible that this discrepancy is caused by a bias in our results from the tail of the disk component, though we see no evidence for this within the sampled z-range (see Fig. C.2).

The two weighting schemes produce consistent eDIG scale heights in the northern half of IC 1553 but diverge above the prominent star-forming clusters in the south, where the logarithmic weighting results in scale heights up to ∼50% larger than when using the linear weighting. This indicates that the eDIG above the active star formation produces more emission at high altitudes than expected from a single-component exponential behavior. Since Rautio et al. (2022) sample almost twice the z-range, it is possible that their results are dominated by this more extended component that we only detect above the star-forming clusters, which could also cause the discrepancy. This would indicate that the profile of the eDIG halo at larger heights may be dominated by the same process that shapes the eDIG profile above the star-forming region on the southern side of the galaxy at lower heights.

If the emitting ionized gas is optically thin, then the Hα flux in each bin is given in Rayleighs by the line-of-sight integral (e.g., Reynolds 1991b; Dettmar 1992)

(8)

The increased emission at high-altitude therefore corresponds to regions along the line-of-sight that are denser and/or cooler compared to expectations from the innermost disk-halo interface. Since this extended emission coincides with the actively star-forming regions in IC 1553, gas transport from feedback processes may play a role here. Additionally, Rautio et al. (2022) suggest that the galaxy is currently accreting gas from a diffuse source based on the kinematics of the ionized gas, which could also cause overdensities at large z.

Effects from a relatively low gas temperature at higher altitudes are also interesting to consider, given that other studies have found that the observations instead require additional heating mechanisms beyond ionizing radiation from planar OB stars for the eDIG (e.g., Reynolds & Cox 1992; Reynolds et al. 1999; Collins & Rand 2001; Flores-Fajardo et al. 2011). This additional heating is usually invoked to explain the near constancy of [S II]/[N II] and the rising trend of [O III]/Hβ (or equivalently [O III]/Hα assuming an approximately constant Balmer decrement) with increasing altitude. However, neither of these trends is observed in the southern eDIG of IC 1553 (compare with Fig. 6). If gas accretion is indeed relevant here, lower temperatures at larger heights may be indicative of significant mixing of the eDIG with a colder phase of the accreted material.

4.2. Ionization structure

The most intriguing property of IC 1553 related to the ionization state of its eDIG is certainly the [N II]λ6583 emission. The bins with the highest flux relative to Hα are not directly spatially related with the potential sources of ionization and instead form a slightly off-centered conical structure. The usually observed DIG property that [S II]λ6717/[N II]λ6583 remains roughly constant with increasing altitude only holds for this conical volume. A rising trend of this ratio with increasing altitude – like we see outside of this cone – has been observed before (Golla et al. 1996), but never before in tandem with the usual constancy within the same galaxy. Outside of this conical region the gaseous halo of IC 1553 exhibits a deficit of [N II]λ6583 emission with respect to [S II]λ6717 compared to the radiation transfer models of Wood & Mathis (2004). We note that most studies of the DIG were conducted with long-slit spectroscopy, so that in a case like IC 1553 the analysis will yield dramatically different interpretations depending on where exactly the slit was placed. It might therefore be necessary to revisit the target galaxies of these previous studies using IFU data, to extract the full picture of the relevant physical processes.

This behavior can be explained in principle, if the primary ionization mechanism for the eDIG of IC 1553 is actually shock ionization from stellar feedback. Shocks have been invoked before as a significant contributor to the maintenance of the ionization state in the DIG of dwarf galaxies (Martin 1997) and in the warm component of galactic winds (e.g., Sharp & Bland-Hawthorn 2010; Rich et al. 2010). Following the Allen et al. (2008) models, slow shocks naturally produce a lower [N II]λ6583 flux relative to [S II]λ6717 when compared to photoionization (compare with Fig. 8), assuming Dopita et al. (2005) abundances, resulting in the observed deficit. The observed rising trend of this emission line ratio with altitude in the outskirts of the halo can then be expected from a mixing sequence between the predominant shock ionization in the halo and photoionization in the disk. The evident correlations between the diagnostic line ratios and the Hα FWHM (see Fig. 9), which is valid for the majority of observed bins, certainly supports this idea, and the ionization diagnostic diagrams are consistent with this scenario as well. It is unclear, however, how the observed morphology of the [N II]λ6583 emission would be produced, since the highest fluxes relative to Hα occur closer to the galactic center than the origin of the outflow.

Alternatively, if photoionization from planar OB stars is responsible for the maintenance of the eDIG ionization state, one needs to explain the [N II] emission deficit outside of the conical structure and the [N II] morphology. The low [N II] flux can be the result of deficient nitrogen abundances in the eDIG or the ionization conditions favoring other ionization states for nitrogen. Figure 12 shows maps of the abundance diagnostic ratios

(9)

(10)

Fig. 12. Maps of the metallicity diagnostic ratios log([N II]λ6584/{[S II]λ6717,6731}) (N2S2, left) and log({[S II]λ6717,6731 + 3.44 × [S III]λ9069}/Hβ) (S​23, right).

though as [S III]λ9531 lies outside the MUSE spectral range we assume its flux is approximately equal to 2.44 × [S III]λ9069 from atomic physics (Pérez-Montero 2017).

The central conus has an N2S2 value similar to that of the thin disk, leading to a slim, cross-shaped structure with values in the range N2S2 ≈ −0.1 to −0.18, with a few localized knot-like regions reaching −0.05. Using the calibrations of Kewley & Dopita (2002) with an oxygen abundance of log(O/H) ≈ 8.45 (see Table 1) this gives a lower limit on the ionization parameter in these regions of q ≥ 2 × 10⁷ cm s⁻¹. This cross-shape appears encased in a thicker layer of gas with N2S2 ≈ −0.2 to −0.3, which would indicate a lower ionization parameter around 1 − 2 × 10⁷ cm s⁻¹. In this map IC 1553 appears fairly symmetrical; the biggest difference between the northern and southern disk is that this thicker layer is slightly less extended in the north. Outside of this layer, N2S2 decreases further to values between −0.3 and −0.4, corresponding to a smooth further decrease of q down to 0.5 − 1 × 10⁷ cm s⁻¹.

In S​₂₃, however, the eDIG of IC 1553 exhibits a pronounced asymmetry. Outside of the thin disk S​₂₃ appears to have an almost regular gradient along the d-axis from S​₂₃ ≈ 0.2 in the south up to 0.75 in the region of the northern outflow. Again using the Kewley & Dopita (2002) calibration, the ionization parameters inferred from the N2S2 map at the same oxygen abundance can only produce S​₂₃ ≤ 0.35, with higher ratio values requiring much smaller ionization parameters. This is contradicting the inference from the N2S2 map for the central eDIG, where a large q is required to explain the central conus, and produces a discrepancy in the required ionization parameter in the northern eDIG spanning orders of magnitude. The gas within the thin disk exhibits values of S​₂₃ = 0.1 to 0.4, which translates to q ≤ 1.5 × 10⁸ cm s⁻¹ and is consistent with the N2S2 map.

In summary, the N2S2 and S​₂₃ maps together indicate that the state of the ionized gas in the plane of the disk and in the southern eDIG above the prominent star formation is consistent with the Kewley & Dopita (2002) photoionization models, while the diagnostics are contradictory in the central and northern eDIG. Consequently, photoionization from OB stars is unlikely to be the primary ionization mechanism outside of the southern eDIG.

Rautio et al. (2022) discuss HOLMES as a possible additional ionizing source. The corresponding model grids published by Flores-Fajardo et al. (2011) fit slightly better to the observed emission line fluxes in the [O III]/Hβ vs. [N II]/Hα diagram than the Allen et al. (2008) shock models, since they can reproduce the highest observed [N II]/Hα values. The only other model grid presented by Flores-Fajardo et al. (2011) uses the [O II]λ3726,3729 lines, which are outside of the spectral range of MUSE for local galaxies. While the HOLMES models do provide the best fit to the observed line ratio, it is particularly hard to differentiate sources of ionization in the [O III]/Hβ vs. [N II]/Hα diagram alone because the relevant model grids all fall relatively close to each other and partially overlap. Additionally, photoionization from HOLMES is not expected to produce the observed correlation between emission line ratios and the Hα FWHM and the [N II]/Hα and [S II]/[N II] morphologies are even more puzzling than in the shock ionization scenario, since HOLMES are expected to be fairly smoothly distributed throughout the thick disk and halo. Shocks, therefore, appear to be the more promising ionization mechanism in the northern and central eDIG.

4.3. Models of the gravitational potential of IC 1553

As described in Sect. 3.3, we assembled a set of model gravitational potentials, which reproduce the rotation curve of IC 1553 measured by Comerón et al. (2019). The first panel of Fig. 13 shows a histogram of the resulting total masses in these model potentials. On average the models predict a total mass for IC 1553 of M_tot = (70 ± 12)×10⁹ M_⊙. If we assume its luminous mass is equal to 16.36 × 10⁹ M_⊙ (the sum of all component masses reported by Comerón et al. 2018), then this corresponds to a dark matter fraction of (81 ± 3)%.

Fig. 13. Histograms of the mass distribution properties in the model gravitational potentials. Left: distribution of the total mass in the potential models that reproduce the rotation curve from Comerón et al. (2019) for IC 1553, as described in Sect. 3.3 (R2 ≥ 0.99). Right: distribution of the component masses for the same set of gravitational potential models.

However, the potential models move the majority of the mass into the disk component, as can be seen in the right panel of Fig. 13. On average the halo component contains only ∼10% and the disk component ∼80% of the total mass. It is therefore evident that the total mass distribution in our model potentials is significantly more concentrated in the plane of the disk than would be expected. The average mass of the disk component in this set of models is M_disk = (55 ± 8)×10⁹ M_⊙, which is more than three times as much mass as the 16.36 × 10⁹ M_⊙ measured by Comerón et al. (2018). If the total mass estimated from these models can be trusted then the models may overestimate the escape velocities we present in Sect. 3.3 as a result.

It is possible that this underestimation of the halo component is caused by the fact that the radial range covered by the MUSE data does not extend sufficiently to where the rotation curve plateaus, making it hard to constrain the properties of the halo component.

4.4. Evidence for an off-centered outflow

The off-centered galactic outflow we find in IC 1553 is certainly unusual. The energetic source of galactic outflows and winds discussed in the literature mostly falls into one of two categories: nuclear starbursts, typically in more massive disk galaxies, and global starbursts, often in dwarf galaxies and merging galaxies. Additionally, starburst regions at the end of an elongated stellar main body or bar-like structure appear in Magellanic dwarf galaxies (type IBm) like NGC 2366 (van Eymeren et al. 2007, 2009a), NGC 4861 (van Eymeren et al. 2007, 2009b), or Kiso 5639 (Elmegreen et al. 2016). Global starburst-like activity over the whole disk has also been observed in massive galaxies like NGC 3079 (e.g., Cecil et al. 2001; Li & Wang 2013), where the outflow is not confined to the famous starburst bubble in the core (Cecil et al. 2001). Another special case is the dwarf galaxy NGC 1569, where a wind seems to be driven by the entire stellar disk (Martin et al. 2002; Veilleux et al. 2005). In NGC 1569 a clear age gradient is apparent and the two older star clusters carved a large hole into the gas disk and drive a large scale outflow (Heckman et al. 1995), while the youngest cluster is still embedded in its dusty, dense H II region (Clark et al. 2013). In more massive galaxies, localized star-forming regions driving large (∼1 kpc size) bubbles out of the disks are not uncommon, for example NGC 55 (Graham & Lawrie 1982; Bomans & Grebel 1994) and other edge-on galaxies (e.g., Rossa & Dettmar 2003a). The case of IC 1553, however, may be unique in that it appears to have a large-scale outflow, with its base located far from the dynamical center of the host galaxy but not at the end of the disk. Depending on projection effects, though, the base of the outflow may indeed be located near the end of the disk of IC 1553, similar to what has been observed in some dwarf galaxies.

While the spectral resolution of the MUSE data is not sufficient for us to identify multiple kinematic components in the nebular emission lines, the following indirect evidence for the presence of an outflow was found: (1) The Hα FWHM, which is sensitive to the velocity of gas flows in the presence of a shock, shows a roughly biconical structure centered on a region of the disk approximately 1.4 kpc away from the galactic center in projection. (2) The emission line ratios [S II]/Hα and [O I]/Hα outside of the disk are consistent with shock ionization models, rather than photoionization models. [N II]/Hα is also consistent with this scenario, but cannot distinguish between shock- and photoionization. (3) The shock-sensitive emission line ratios [N II]/Hα, [S II]/Hα and [O I]/Hα are correlated with the Hα FWHM.

Rautio et al. (2022) originally suspected that IC 1553 hosts a large biconical outflow related to superbubbles. The observations listed above together make a strong case that this suspicion is indeed correct. This would make IC 1553 the first known host of a localized, off-centered galactic outflow or wind.

4.5. The disk-halo interface

The brightest spot of the eDIG in Hα at any given altitude in the inner halo is the region above the prominent star formation region in the south. As a rough quantification of this observation we selected fiducial intervals in the radial direction representing the southern (−23″ ≤ d ≤ −13″), central (−5″ ≤ d ≤ 5″), and northern (10″ ≤ d ≤ 20″) regions of IC 1553, and compared the median Hα fluxes in these regions at different heights in Fig. 14. We only considered bins in the northwestern direction above the midplane (where z < 0) because of the contamination from the background galaxy in the field of view.

Fig. 14. Ratios of the median Hα fluxes between the southern, central, and northern fiducial regions of IC 1553, as defined in Sect. 4.5, as a function of the height above the midplane in the northwestern direction (where z < 0). The blue data points give the south-to-center flux ratio and the orange data points give the south-to-north ratio. The dashed line marks the ordinate value of 1.

In terms of this estimate, within the inner 10″ of altitude the southern halo is around ∼2.5 times as bright in Hα as the central halo. The north-to-south comparison at low altitude (−z ≲ 5″) is modulated by the warp in the Hα-disk of IC 1553 at its northern edge and therefore does not represent the general trend. At higher altitudes the ratio has a maximum around the value ∼4. In general, though, both ratios tend to steadily decrease to a value of ≲1 at altitudes above ∼17″. We note that at this height the south-to-center ratio is contaminated by the extraplanar H II region mentioned by Rautio et al. (2022).

According to Eq. (8) the higher Hα flux in the southern eDIG corresponds to a decreased electron temperature and/or increased electron density along the line of sight (assuming that scattered light from the disk is negligible). Wood & Mathis (2004) argued that the (electron) temperature of DIG is expected to increase with distance to the star formation because the radiation field hardens along its path. For this effect to be the sole cause of the asymmetric Hα morphology, the average temperature in the northern eDIG (at |z|≤10″) has to be around times as high as in the southern halo, and the temperature in the central halo has to be around as high as in the south. Temperatures of the eDIG are not well constrained in the literature, but it is a common assumption that the bulk of the gas is in the temperature range 6000–10 000 K (e.g., Haffner et al. 2009), which would not allow for the required temperature ratios. However, IC 1553 is certainly unusual in the distribution of its star formation, so it is conceivable that it produces an unusual temperature structure in its eDIG.

Alternatively, the Hα morphology can be caused by increased densities in the southern eDIG relative to the other regions and, since I_Hα scales with an average , the required difference in the density is significantly lower than it is for the temperature. Such differences in the density could be caused by star formation feedback. In this case, since the southern eDIG produces approximately four times as much Hα flux as the northern eDIG up to the height of z ≈ 13″ (see Fig. 14), this would correspond to twice the average eDIG density. Interestingly, if all ionized gas in the eDIG originates in the disk, then this would imply that the gas transport process into the disk-halo interface must be greater in magnitude for the active star formation in the south than it is for the outflow in the north. Given that the emission line ratios indicate that the shock velocity is on the same order as the escape velocity in our models, it is possible that the bulk of the gas in the outflow is instead transported well into the halo or beyond. Additionally, Rautio et al. (2022) concluded that IC 1553 is currently undergoing accretion of gas from a diffuse source, which could contribute to the density of the extraplanar ionized gas.

In practice, it is likely that both the density- and temperature structure produce the relatively enhanced Hα emission in the southern halo synergetically. It is interesting to note that, while the highest fluxes are therefore linked to the current star formation in the plane of the disk, the ionization of the eDIG appears to be driven in large part by the off-centered outflow on the other side of the galaxy. In this sense, we do not find evidence that the disk-halo interface is dominated by a single process. Rather, it appears to be shaped by star formation on different timescales in different aspects.

5. Summary and conclusions

We analyzed the optical emission line spectrum of the extraplanar DIG in the edge-on star-forming galaxy IC 1553 using archival IFU data from the MUSE spectrograph. We derived the Hα scale height h_eDIG(d) as a function of the projected radial distance from its center and find that it decreases roughly linearly with |d|. The scale height values that we measure for IC 1553 are generally lower than those found by Rautio et al. (2022), which sampled the eDIG to higher altitudes. This discrepancy might indicate the presence of a more extended layer of warm ionized gas, which would appear to be connected to the active star formation in the southern half of IC 1553. Alternatively, it may be caused by peculiarities in the density or temperature structure of the eDIG.

Rautio et al. (2022) proposed that IC 1553 hosts a biconical outflow originating from a superbubble in the disk, which we can confirm. Our analysis suggests that the shock produced by this outflow is a significant, if not dominant, source of ionization for the eDIG. The bulk of the gas transport into the disk-halo interface, however, is most likely coupled to the active star formation in the southern disk rather than the outflow in the north. The spatial separation between the base of this outflow and the current star formation leads to very different ionization conditions across different regions in the eDIG.

We compiled a population of gravitational potential models that reproduce the observed rotation curve for IC 1553 from Comerón et al. (2019) and predict an escape velocity in the range 190–275 km s⁻¹, which encompasses the best-fitting shock velocity of the outflow of 225 km s⁻¹. However, while these models have reasonable total masses for IC 1553, they predict a very large mass fraction of the disk component, which may cause the escape velocity to be overestimated.

Based on the results of this study we argue that, in a case like IC 1553, in which the physical conditions of the eDIG change significantly along the galaxy’s major axis, IFU spectroscopy is a necessity for an accurate analysis of the physical processes involved. The majority of studies concerned with the DIG in other galaxies, however, only had access to long-slit spectra, so it might be necessary to revisit the corresponding target galaxies with the current generation of IFU spectrographs if we wish to understand the role of the DIG in the context of galaxy evolution.

Acknowledgments

We thank Sébastien Comerón for his helpful comments during the review process, which improved this work significantly. R.J.D. and D.J.B. acknowledge funding from the German Science Foundation DFG, via the Collaborative Research Center SFB1491 “Cosmic Interacting Matters – From Source to Signal”. P.K. acknowledges financial support by the German Federal Ministry of Education and Research (BMBF) Verbundforschung grant 05A20PC4 (Verbundprojekt D-MeerKAT-II).

References

Appendix A: Contamination from a background galaxy

The extended source at the southeastern edge of the field of view visible in the V-band continuum and Hα kinematic maps (Figs. 2 and 10) was determined to be a background galaxy based on its spectroscopic redshift of z = 0.052 corresponding to a distance of ≳200 Mpc, whereas IC 1553 has a distance of only 36.5 Mpc (Tully et al. 2016). The [OI]λ6300 emission line of this background galaxy is redshifted to the same wavelength as Hα at the systemic velocity of IC 1553, which leads to false detections in Fig. 10 appearing like a satellite or coherent cloud of ionized gas. We have masked the corresponding area in figures showing any diagnostic that uses Hα throughout this paper as a result.

The background galaxy’s spectrum contains typical nebular emission lines produced by star formation and has an estimated maximum rotational velocity of around 90 km/s. A comparison of example spectra from the central regions of IC 1553 and this background galaxy is shown in Fig. A.1.

Fig. A.1. Comparison of example spectra of the central regions of IC 1553 (black) and the background galaxy (red) extracted from the spectral cube.

Appendix B: Calculation of the η′-parameter

The calculation of η′ as described in Sect. 3.2 is largely analogous to the determination of η in Erroz-Ferrer et al. (2019). The largest difference is that Erroz-Ferrer et al. (2019) consider two demarcation lines (the empirical and theoretical curves from Kauffmann et al. 2003; Kewley et al. 2001) resulting in three separate regions in the [O III]/Hβ vs. [N II]/Hα diagram, while we only have the equivalent theoretical limit from Kewley et al. (2001) in the [O III]/Hβ vs. [S II]/Hα diagram.

Following the nomenclature in Erroz-Ferrer et al. (2019), the distance, d₁, between a bin’s position in the ionization diagnostic diagram P₁ = (x₁, y₁) and any point (x, f(x)) on the curve from the Kewley et al. (2001) starburst limit f(x) = 0.72/(x − 0.32)+1.30 is

(B.1)

The minimal distance therefore corresponds to the solution of ∂d₁(x)/∂x = 0, or equivalently to the solution of

(B.2)

Out of the two real roots of this polynomial, we selected the one that produces the smaller distance, d₁(x₀). Following the definition of the η-parameter, we required that η′ = 0.5 if P₁ ∈ f(x), so that

(B.3)

Appendix C: Additional figures

Figure C.1 shows the line-of-sight velocity component for all bins in the midplane of IC 1553 as a function of the projected angular distance to the dynamic center, measured on Hα. Following Comerón et al. (2019), we show rotation curve fits to the expression (Courteau 1997)

(C.1)

Fig. C.1. Line-of-sight velocity profile for all the bins intersecting the midplane of IC 1553, measured on the Hα emission line (gray circles). The dashed lines mark fits of Eq. (C.1) to the data in this work (dashed blue line) and from Comerón et al. (2019; dashed red line). The solid yellow curve marks the rotation curve from the best-fitting model of the gravitational potential obtained, as described in Sect. 3.3.

where v_c is the circular velocity, χ is the normalized axial coordinate, and γ determines the sharpness of the transition from the inner increase in the rotation curve to its flattening further out. The curves shown in Fig. C.1 include the best fit to the measured velocity profile in this work (v_c = 134 km s⁻¹,χ = 13.3″,γ = 2.81), the rotation curve obtained by Comerón et al. (2019) (v_c = 142 km s⁻¹,χ = 13.5″,γ = 2.02), and the rotation curve of the best-fitting gravitational potential model (see Sects. 3.3 and 4.3).

Figure C.2 visualizes the Hα scale height fits presented in Sect. 3.1 for all projected radial positions, analogous to the left panel of Fig. 5.

Fig. C.2. Visualization of the Hα scale height fit analogous to the left panel of Fig. 5, shown for all projected radial positions, d.

Fig. C.2. continued.

All Tables

All Figures

Fig. 1. Visualization of the busy function fit to determine the systemic velocity of IC 1553. The integrated Hα profile was obtained by collapsing the data cube along both spatial axes for bins with altitude |z|≤10″. The central wavelength of the profile is taken as the mean of the centroid wavelengths of the red- and blue-shifted components in the busy function fit. We estimate the uncertainty in this determination as ±0.5 channel width.

In the text

	Fig. 2. Synthetic V-band images constructed from the data cube at the original resolution (left), and after rebinning with a regular 50 × 50 grid of resolution elements (right). The dashed line and the white circle are the estimates for the plane of the disk and the dynamic center of IC 1553, respectively.
In the text

Fig. 3. Maps of the Balmer emission in IC 1553. Left: map of the Balmer decrement Hα/Hβ obtained from the spectral template fitting with pPXF. Middle: map of the extinction-corrected Hα emission line flux in Rayleighs with linear scaling. Right: the same as the middle panel, but with logarithmic scaling. The contours are arbitrary equidistant levels from the synthetic V-band image meant to give orientation for the eye. The white square in the lower left corner corresponds to the masked background galaxy. The red circles in the middle panel mark the H II regions used for the oxygen abundance measurement listed in Table 1. The boxed region marked A in the right-hand panel is identical to the gray shaded area in Fig. 5 below.

In the text

Fig. 4. Synthetic V-band image of IC 1553 overlaid with equidistant contours from the W4-band image from the ALLWISE Data Release at 22 μm (left), from the RACS-low image at 887.5 MHz (middle), and from the Hα image extracted from our data cube smoothed with a Gaussian kernel (right). The nominal angular resolutions are 12.0″ for the WISE W4-band, 15.0″ for RACS-low, and the Hα image was convolved with a 12.0″-wide Gaussian for comparison.

In the text

Fig. 5. Visualization and results of the Hα scale height fitting. Left: visualization of the IHα fit at the example projected radial distance of d = −15.0″. The dotted lines mark the separation height, zsep = 8″, between bins used for fitting the disk H II (red curves) and eDIG (blue curves) components. The solid curves are derived with linear weighting of IHα according to Eq. (4) and the dashed curves result from logarithmic weighting of IHα during the fit. Right: radial distribution of scale heights for the eDIG (blue squares from logarithmic weighting, green circles from linear weighting) and disk H II (red crosses from logarithmic weighting, orange crosses from linear weighting) components obtained from the fitting. The shaded areas around the scale height curves give the 1σ errors obtained from the corresponding estimated parameter covariance. The dashed lines are best fits of Eq. (5) to the eDIG component scale heights within |d|≤20″. The gray shaded area, A, marks the region where the logarithmic weighting produces significantly larger eDIG scale heights than the linear weighting.

In the text

Fig. 6. Maps of the emission line ratios commonly used in the classical diagnostic diagrams. First row: [N II]λ6583/Hα (left), [S II]λ6717,6731/Hα (middle), and [S II]λ6717/[N II]λ6583 (right), which is commonly analyzed in studies of the DIG. Second row: [O I]λ6300/Hα (left) and [O III]λ5007/Hβ (middle), as well as the η′ parameter (right) as described in Sect. 3.2 and defined in Appendix B, which quantifies the distance to the Kewley et al. (2001) starburst limit in the [O III]/Hβ vs. [S II]/Hα diagnostic diagram. The black demarcation lines in the η′ map separate bins below the starburst limit (η′≤0.5) from bins above that limit (η′> 0.5). The gray and white contours are at the same equidistant arbitrary levels from the synthetic V-band image as in Fig. 3.

In the text

Fig. 7. Typical emission line diagnostic diagrams for the analysis of ionization sources. The data points are the population of bins with Hα flux ≥5σ, colored by their η′ value. The black dotted line marks the theoretical limit for starburst galaxies from Kewley et al. (2001). The solid and dashed lines mark the metallicity (Z) and ionization parameter (q) grid of photoionization models described in Levesque et al. (2010) for the best-fitting star formation age of 0 Myr and electron density of 100 cm−3.

In the text

Fig. 8. Same diagnostic diagrams as shown in Fig. 7 (dotted black and red lines are the starburst limit and photoionization models described therein), now including the velocity (v) and magnetic field (B) grid (dashed and solid lines) for the best-fitting shock models from Allen et al. (2008), in which the preshock density n = 1 cm−3 and the chemical abundances described in Dopita et al. (2005) were used.

In the text

	Fig. 9. Shock-sensitive emission line ratios [N II]λ6583/Hα, [S II]λ6717, 6731/Hα, and [O I]λ6300/Hα vs. the Hα FWHM of all bins with Hα flux ≥5σ, colored by their η′ value.
In the text

Fig. 10. Maps of the emission line kinematics of the DIG in IC 1553. Top left: map of the Hα line-of-sight velocity in km s−1 with contours from the continuum map. Top right: map of the Hα FWHM in km s−1 after correcting for LSF broadening with Eq. (2). Bottom row: FWHM maps measured on [S II]λ6717 (left) and [N II]λ6583 (right), LSF-corrected in the same way. The white contours are from the continuum map as before. The violet contours visualize the alignment with the respective emission line flux by separating bins based on the following flux thresholds relative to Hα: [S II]λ6717,6731/Hα ≷ 1.0 (left) and [N II]λ6583/Hα ≷ 0.60 (right).

In the text

	Fig. 11. Gaussian kernel density estimates of the escape velocity, vesc, at the position of the spaxel with the highest Hα FWHM derived from fitting gravitational potentials to the rotation curve of IC 1553. vesc was modeled for multiple projection angles, φ, between the plane of the sky and the deprojected positional vector from the galactic center to this spaxel.
In the text

	Fig. 12. Maps of the metallicity diagnostic ratios log([N II]λ6584/{[S II]λ6717,6731}) (N2S2, left) and log({[S II]λ6717,6731 + 3.44 × [S III]λ9069}/Hβ) (S​23, right).
In the text

	Fig. 13. Histograms of the mass distribution properties in the model gravitational potentials. Left: distribution of the total mass in the potential models that reproduce the rotation curve from Comerón et al. (2019) for IC 1553, as described in Sect. 3.3 (R2 ≥ 0.99). Right: distribution of the component masses for the same set of gravitational potential models.
In the text

	Fig. 14. Ratios of the median Hα fluxes between the southern, central, and northern fiducial regions of IC 1553, as defined in Sect. 4.5, as a function of the height above the midplane in the northwestern direction (where z < 0). The blue data points give the south-to-center flux ratio and the orange data points give the south-to-north ratio. The dashed line marks the ordinate value of 1.
In the text

	Fig. A.1. Comparison of example spectra of the central regions of IC 1553 (black) and the background galaxy (red) extracted from the spectral cube.
In the text

Fig. C.1. Line-of-sight velocity profile for all the bins intersecting the midplane of IC 1553, measured on the Hα emission line (gray circles). The dashed lines mark fits of Eq. (C.1) to the data in this work (dashed blue line) and from Comerón et al. (2019; dashed red line). The solid yellow curve marks the rotation curve from the best-fitting model of the gravitational potential obtained, as described in Sect. 3.3.

In the text

	Fig. C.2. Visualization of the Hα scale height fit analogous to the left panel of Fig. 5, shown for all projected radial positions, d.
In the text

	Fig. C.2. continued.
In the text