© The Authors 2024

1. Introduction

The study of galaxy kinematics provides an invaluable tool for understanding the processes that shape the evolution of galaxies. Specifically, analyzing the motion of ionized gas within these celestial structures uncovers patterns indicative of both their current dynamical state and historical phenomena, such as possible rotations and feedback processes.

This study is focused on I Zwicky 18 (hereafter IZw18), which has been identified as a blue compact dwarf (BCD) galaxy (e.g., Petrosian et al. 1997; Tikhonov 2007; Ramos et al. 2011) located approximately 18 Mpc away (Lebouteiller et al. 2017). IZw18 has attracted attention due to the extremely low metallicity of its ionized gas (12+log(O/H) = 7.11, approximately 3% solar; e.g., Kehrig et al. 2016), which places IZw18 among the most metal-poor star-forming systems in the local Universe. Such a low metallicity suggests the presence of chemically unevolved stellar populations and makes IZw18 a prime candidate for investigating galaxy formation and evolution processes that occurred in the primordial Universe.

IZw18 has a highly irregular morphology, a characteristic common among BCDs, which further enhances its scientific interest. While its morphology could be the result of recent starburst activity or interactions with other galaxies, the underlying mechanisms remain a subject of active investigation (e.g., Legrand et al. 1997; Legrand 2000; Kehrig et al. 2015).

The kinematics of IZw18 has been a subject of diverse investigations that have enriched our understanding, yet have left many questions open. Early works (Petrosian et al. 1995, 1997) noted velocity dispersion values (∼20 km s⁻¹) higher than expected from a self-gravitational model (Terlevich et al. 2015). Furthermore, potential rotational motion in the neutral hydrogen component has been detected, leading to speculative theories about the existence of a dark matter halo around the galaxy (van Zee et al. 1998). These are far from the only intriguing features observed in IZw18.

A study employing long-slit, echelle spectroscopy, specifically exploring the spatial changes in the Hα profile, scrutinized the ionized gas within IZw18 in great detail (Martin 1996). This investigation unveiled evidence of a supergiant shell, manifesting as a bipolar outflow in the southwest (SW) and northeast (NE) extensions of the galaxy, expanding at speeds between 30 km s⁻¹ and 60 km s⁻¹. However, a subsequent study by Petrosian et al. (1997) presented a contrasting view. They suggest that the perceived bipolar outflow might not be a singular, coherent structure but rather a conglomerate of clumpy HII regions. These regions, presumably powered by star formation, challenge the earlier interpretation of a supergiant shell, indicating a more fragmented and fractal star-forming environment within the galaxy.

In IZw18, disparities in gas kinematics between Hα and HeII spectral lines are evident. Notably, the HeII lines show a redshift of about 30 km s⁻¹ compared to the Hα in regions with peak HeII emission (Rickards Vaught et al. 2021). This discrepancy, potentially linked to supernova remnants or young stellar clusters, suggests a complex and diverse kinematic structure within the galaxy.

IZw18 has been studied as part of a collective analysis of BCDs as it exemplifies this class of galaxy, characterized by low metallicities (12+log(O/H) between 7.1 and 8.3), irregular and clumpy morphologies, high star formation rates (SFRs), compact structures, and significant gas richness (e.g., Thuan & Martin 1981; Pilyugin 1992; Izotov & Thuan 1999). These features, indicative of a young stellar population and ongoing star formation processes, offer insights into the early stages of galactic evolution in low-mass environments. In addition, these galaxies often exhibit chaotic kinematics in HI, overlaying a disk rotation (Lelli et al. 2012). In a broader examination of extremely metal-poor galaxies (XMPs; i.e., Z < Z_⊙/10), recent studies (e.g., Olmo-García et al. 2017; Isobe et al. 2023; Xu et al. 2024) have observed typical features, such as rotation in ionized gas, shell outflows, chemical inhomogeneity (although this has not been seen in IZw18; Kehrig et al. 2016), and a high loading factor (over 10), all of which are indicative of substantial matter ejection relative to the SFR. These findings underscore that XMPs are systems influenced by both inflow and outflow mechanisms. The extremely low metallicities in such galaxies could be attributed to the accretion of pristine gas driven by these inflows. Outflows are likely a result of intense stellar feedback, amplified by the low gravitational pull in these low-mass galaxies (e.g., Amorín et al. 2012; Bosch et al. 2019; Hogarth et al. 2020). This is evidenced by the broad components in the emission line profiles (Olmo-García et al. 2017).

The observational data presented here offer a unique opportunity to deepen our understanding of the ionized gas kinematics of IZw18. The data were obtained using the Multi-Espectrógrafo en GTC de Alta Resolución para Astronomía (MEGARA; Gil de Paz et al. 2016), an advanced integral field spectrograph capable of capturing detailed spectral data across a 2D field of view (FoV) with high spectral resolution.

By analyzing the spectral data of the main body (MB) and two regions in the halo of IZw18, this study seeks to generate a comprehensive kinematic portrait of this intriguing galaxy. This investigation constitutes a step forward in our understanding of the complex kinematics in dwarf galaxies, paving the way for future detailed studies on similar low-mass starburst systems in the Universe.

This paper is organized as follows. In Sect. 2 we detail the observations and data reduction techniques employed. Section 3 is devoted to the 2D kinematic analysis of the galaxy. The integrated spectra, highlighting spectral features due to an enhanced signal-to-noise ratio (S/N), are presented in Sect. 4. In Sect. 3.3 we discuss the luminosity-weighted velocity and velocity dispersion analysis. Finally, Sect. 6 summarizes our main conclusions.

2. Observations

The observational data for this study were acquired using MEGARA, situated at the Folded Cass F focus of the 10.4 m Gran Telescopio Canarias (GTC) telescope at La Palma Observatory (Gil de Paz et al. 2016). The integral field unit of MEGARA, also known as the Large Compact Bundle (LCB), offers a FoV of 12.5 × 11.3 arcseconds, consisting of 567 hexagonal spaxels, each measuring 0.62 arcseconds. What sets MEGARA apart as a uniquely suitable instrument for this research is its capability to provide 2D spectroscopic data with high spectral resolution. This enables unprecedented detailed 2D kinematic studies of the ionized gas in IZw18, a level of analysis that has not been conducted before. Consequently, the configuration and capabilities of MEGARA make it an invaluable tool for in-depth examinations of individual regions within nearby galaxies.

Our observations focused on three distinct regions within the galaxy IZw18: the MB and two external regions corresponding to the NE and SW zones of the halo of galaxy. These three regions are represented in Fig. 1.

Fig. 1. High contrast composite image of IZw18 using HST ACS R, I, and V bands, shown in the red, green, and blue channels, respectively. Isocontours delineate the central knots (adapted from Papaderos & Östlin 2012). The three LCB pointings are shown. North is up, and east is to the left here and in all subsequent figures.

The data reduction process was performed using the MEGARA data reduction package (Pascual et al. 2019), version 0.23.2. This software pipeline allowed for sky and bias subtraction, flat-field correction, spectra tracing and extraction, correction of fiber and pixel transmission, and wavelength calibration. Regarding the sky-subtraction, we utilized the sky fibers available in the MEGARA instrument, which are specifically designed for this purpose. MEGARA uses a set of eight dedicated sky fibers that are strategically placed at the edges of the full FoV, which spans 3.5 arcmin. These sky fibers allow for an accurate measurement of the sky background at a considerable distance from the central object of interest. As the primary objective was a kinematic analysis, no flux calibration was applied. Thus, throughout this article, all luminosity values are given in arbitrary units.

As for the specific settings of the MEGARA instrument during the observations, the LR-R volume-phased holographic was employed. This covered a spectral range of 6300Å–7363Å with a spectral resolution of R ∼ 6000. For the Hα line at a wavelength of 6578 Å (redshifted wavelength), this corresponds to a spectral resolution of approximately 1.1 Å. In terms of velocity, this spectral resolution is equivalent to approximately 50 km s⁻¹. In Table 1, we present the log of observations.

3. Main kinematic analysis

In this section we present the main results of this study, which are the maps of luminosity, velocity, and velocity dispersion based on the Gaussian fit of the Hα profile. The methodology followed to construct these maps is also presented.

3.1. Methodology

The data analysis for this study was executed on a total of six data cubes. These were categorized into two sets: the first set comprised three data cubes (one per pointing) with sky subtraction, while the second set correspond to the same data cubes but for which no sky subtraction was performed.

In the initial phase of our analysis we observed that the center of the [OI]λ6300 sky line presented spaxel to spaxel variations (up to 0.7 Å). This led us to hypothesize that the wavelength calibration might not be entirely accurate, at least not to the precision we require for our kinematic analysis. Consequently, we decided to use the [OI]λ6300 sky line as a reference point for calibrating the velocity of the Hα line. For now, to determine the observed central wavelength of each line, we employed a simple Gaussian fit to each spaxel.

Thus, we defined a reference wavelength field (λ_ref), where the center of the [OI]λ6300 sky line in each spaxel was taken as the reference wavelength. This λ_ref was subtracted from the observed Hα wavelength (in all the data cubes) field to compute the intrinsic wavelength field of the Hα line. Importantly, using this methodology, we ensured that any potential inaccuracies in the wavelength calibration were minimized. By referring the position of the Hα line to the sky line of [OI]λ6300, we generated a more accurate wavelength field in Hα. This effectively neutralized any potential systemic shifts introduced by wavelength calibration uncertainties.

At this point, attention was shifted to the cubes with sky subtraction. The Hα line profiles were examined, and Gaussian fits were applied. Initially, we started with a single Gaussian fit to the Hα line profile in each spaxel. However, the observation of broader and more complex profiles in some spaxels led us to explore dual-component fitting. The Hα profile of one of these spaxels exhibits a distinctively broadened base with a wavelength shift (see Fig. 2), indicative of complex underlying kinematic processes.

Fig. 2. Hα profile of a single spaxel. The black line represents the profile. The gray line is a Gaussian with the instrumental FWHM. The vertical line highlights the Hα wavelength.

In this way, for each spaxel, both one- and two-component fits were performed, yielding line luminosity, velocity, and velocity dispersion from each Gaussian. Luminosity was calculated from the Gaussian area. To convert the wavelength field to a velocity field, the relativistic Doppler effect was used:

$$\begin{array}{c}\hfill v=c\left(\frac{{(\lambda /{\lambda }_0)}^2-1}{{(\lambda /{\lambda }_0)}^2+1}\right),\end{array}$$(1)

where v is the velocity, c is the speed of light, λ is the observed wavelength, and λ₀ is the rest wavelength. This Doppler effect conversion was consistently applied in all subsequent velocity calculations.

In order to obtain the corrected velocity dispersion (σ), we took into account the instrumental and thermal widths (σ_inst, σ_thermal). The correction was performed as follows:

$$\begin{array}{c}\hfill {\sigma }^2={\sigma }_{\mathrm{obs}}^2-({\sigma }_{\mathrm{inst}}^2+{\sigma }_{\mathrm{thermal}}^2).\end{array}$$(2)

The instrumental width was calculated from the width of the [OI]λ6300 sky line (σ_inst = 0.45 Å). Meanwhile, σ_thermal was calculated by using the Maxwell-Boltzmann distribution:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{thermal}}^2=\frac{K{T}_e}{m},\end{array}$$(3)

where K is the Boltzmann constant, T_e the electronic temperature retrieved from Kehrig et al. (2016) (with a value of 20 600 K) and m is the hydrogen mass.

In addition to instrumental and thermal widths, we could consider the effects of fine structure and natural broadening (σ_fs and σ_natural). Fine structure, leads to a subdivision of emission lines that can be interpreted as a broadening. The σ_fs is only well characterized for atoms with one electron and is usually around 2–3 km s⁻¹ (e.g., García-Díaz et al. 2008), so for our spectral resolution we could ignore this effect. Natural broadening, governed by the Heisenberg uncertainty principle (ΔE ⋅ Δt ≥ ℏ/2) in our case is minimal. Lifetimes of excited states of the H atom are about 10 nanoseconds, leading to a natural linewidth broadening of only a few m/s, which is significantly smaller than other broadening mechanisms encountered in astrophysical observations.

In our one-component fit analysis, certain spaxels are excluded to maintain data integrity. This exclusion is based on three criteria: spaxels where the full width at half maximum (FWHM) of the Hα line is narrower than the instrumental FWHM, discrepancies between the Hα and Hβ line profiles, and a S/N below 3, where signal represents the peak flux of the Halpha line and Noise denotes the standard deviation of the continuum beneath the line. Those spaxels were tacked away from the representation of the luminosity, velocity and the velocity dispersion maps. In the two-component fitting process, we applied two extra conditions. One crucial condition was using the Akaike information criterion (AIC). If any Hα profile (at a particular spaxel) was better fitted by the single-component fit as indicated by a lower AIC value, it was excluded from the dual-component analysis. Additionally, any Gaussian with an amplitude less than 2.5 times the standard deviation of the continuum was disregarded in the two-component fit.

Due to the high precision required for the two-component fits, all the data cubes were spatially smoothed prior to the fitting process using a Gaussian kernel with a FWHM of 1.4″. This specific choice of FWHM was arrived after exploring several kernel sizes, aiming to optimize the tradeoff between enhancing coherence among adjacent spaxels and minimizing the loss of spatial resolution. The smoothing process, therefore, not only improved the coherence between adjacent spaxels but did so while carefully considering the inevitable compromise in spatial resolution. Ebeling et al. (2006) show the high utility of using Gaussian kernels for smoothing astronomical images.

Finally, it was necessary to distinguish between the primary and secondary components of each dual-component fit. This distinction was made using two criteria: the residual area after subtracting the two Gaussians of the dual fit from the single-component Gaussian, and the absolute velocity difference between each of the dual component velocities and the single component velocity.

For the velocity difference criterion, we defined the velocity of the single-component fit as v_single, and the velocities of the two components of the dual fit as v_1, dual and v_2, dual. Then, we computed the absolute velocity differences, |v_single − v_1, dual| and |v_single − v_2, dual|. The component from the dual fit with the smallest absolute velocity difference was identified as the primary component, thus ensuring that the primary component velocity was as similar as possible to the single component velocity.

For the residual area criterion, we compared the residual areas, where G_single(λ) represents the Gaussian fit for the single-component model, and G_1, dual(λ) and G_2, dual(λ) are the Gaussian fits for the first and second components of the dual fit, respectively. These residual areas are defined as A_residual1 = ∫|G_single(λ)−G_1, dual(λ)|dλ and A_residual2 = ∫|G_single(λ)−G_2, dual(λ)|dλ. The component with the lower residual area, more closely aligning with the single-component fit, was designated as the primary component.

As a final remark, every calculation of the velocity is relative to the velocity of the center of mass of the ionized gas emitting Hα. Since the mass of gas is proportional to the luminosity, we can set the velocity of the center of mass as the zero velocity.

3.2. Luminosity, velocity, and velocity dispersion maps

The result of this methodical process was a comprehensive suite of nine maps including each of the three pointings (MB, NE halo, and SW halo). Each one displaying luminosity, velocity, or velocity dispersion of the Hα line for the three different components (the two components of the dual fit and the one of the single fit).

The maps are designed to offer a holistic yet detailed view of the different regions of the galaxy. Each map features a background layer sourced from deep Hubble Space Telescope (HST) observations (same as used in Fig. 1), which cover a wide FoV to capture the galaxy in its entirety. This panoramic backdrop serves to contextualize the more targeted data overlays from the MEGARA/GTC observations. Specifically, three distinct overlays are superimposed on the HST background in each map, each corresponding to one of the observed regions of the galaxy: the MB, the SW halo, and the NE halo.

3.2.1. Single component

We began our exploration with the maps derived from a single-component Gaussian fit. Initially, we examined the luminosity map, depicted in Fig. 3, which delineates the distribution of Hα luminosity across the different regions of the galaxy IZw18. The Hα luminosity distribution serves as a bedrock for comprehending the structure of the ionized gas within the galaxy. The luminosity map of the IZw18 MB distinguishes three, well-known main features: two star-forming knots, one situated in the southeast and the other in the northwest, and an arc-like structure (“plume”) that has one end rooted in the vicinity of the NW knot (e.g., Davidson et al. 1989; Dufour & Hester 1990; Dufour et al. 1996; Vílchez & Iglesias-Páramo 1998; Papaderos & Östlin 2012; Kehrig et al. 2016). These features underscore the known morphology of IZw18, presenting a lucid portrayal of the distribution and nature of ionized gas, which is particularly insightful for our kinematic study. When focusing on the two halos, it is apparent from the luminosity map that their luminosity is significantly diminished compared to that of the MB (around 1 or 2 orders of magnitude). Moreover, these maps also trace the outer shells of IZw18, as we can see by comparing them with the background HST image, thereby providing a nuanced understanding of the ionized gas distribution and its overarching structure.

Fig. 3. Hα luminosity maps for the MB and halos obtained from the single-component fit. The background layer, presented in a cyclical gray scale, represents the wide-field deep image from the HST using the f814w and f606w filters of the WFC2 instrument. This background is used in all the images of IZw18 in which the MB and both halos are represented.

Following the examination of the luminosity map, the focus transitions to the velocity map illustrated in Fig. 4. This map emerges as an indispensable tool for unraveling the kinematic behavior of the ionized gas within IZw18, offering insights into the internal motion within both the MB and the halos, and unearthing patterns and flows that are pivotal to understand the dynamical state of the galaxy.

Fig. 4. Hα velocity maps for the MB and halos obtained from the single-component fit. The contours correspond to the Hα luminosity in log scale obtained from the same fit component. The same contours are represented in all Hα velocity and velocity dispersion maps.

Within the MB, the velocity map distinctly unveils a rotation pattern between the two central knots, with a velocity difference of ∼40 km s⁻¹. This harmonize with findings from several studies that have explored such kinematics in Hα and HI using the 21 cm line (Petrosian et al. 1997; van Zee et al. 1998; Lelli et al. 2012; Rickards Vaught et al. 2021; Isobe et al. 2023; Xu et al. 2024). Notably, this rotational feature seems to dissipate toward the NE and SW peaks of the MB map, where a shift toward redder velocities (i.e., positive velocities or gas moving away from us) is observed, indicative of outward or recessionary motion.

Transitioning to the halos, the velocity maps presented herein are pioneering, unveiling for the first time the 2D kinematics of ionized gas in these outer realms of IZw18 and tracing the movement of the shells enveloping the galaxy. Previous studies (e.g., Rickards Vaught et al. 2021) studied the kinematics of the Hα line up to ∼0.75 kpc from the center of the galaxy. In our study we observe these two regions of the halo up to ∼3 kpc away from the center. Particularly in the NE halo, a relatively uniform map is observed, with velocities leaning toward the blue, albeit with some velocity gradients present. Conversely, the SW halo exhibits a more complex kinematic structure, with gas manifesting a large-scale turbulent motion, possibly hinting at underlying dynamical processes or interactions occurring within this region.

Wrapping up our exploration of single-component maps, the velocity dispersion maps depicted in Fig. 5 merits a close examination. This map elucidates the spectrum of velocities inherent within each region, serving as a discerning gauge of areas of dynamic stability or potential turbulence.

Fig. 5. Hα velocity dispersion maps for the MB and halos obtained from the single-component fit.

A foundational aspect to consider is the theoretical velocity dispersion IZw18 would exhibit if the gas in the galaxy adheres to the virial theorem, as delineated in Terlevich et al. (2015). We retrieved the Hβ integrated luminosity from the Hβ flux (1.59 × 10⁻¹³ erg s⁻¹ cm⁻²) in Kehrig et al. (2016) and the distance (18.2 Mpc) in Aloisi et al. (2007). From this Hβ luminosity (6.3 × 10³⁹ erg s⁻¹) we conclude, in accordance with Terlevich et al. (2015), that if the observed velocity dispersion surpasses the threshold of ∼20 km s⁻¹, it signifies the presence of another process (beyond the natural random motion induced by gravitational interaction) that escalate the kinetic energy of the gas, thereby augmenting the velocity dispersion.

In scrutinizing the map, it is notable that a majority of the regions encompassing the MB and the NE halo exhibit velocity dispersions < 25 km s⁻¹. This suggests that, broadly, these velocity dispersions could predominantly be attributed to the virial motion of the gas. However, within the MB, particularly in the top left and bottom right corners (see again Fig. 5), velocity dispersions surge to up to 40 km s⁻¹. This escalation hints at broadening or added complexity in the Hα profile in these regions. Regarding the halos, the NE halo presents a pretty mild velocity dispersion that does not go above 35 km s⁻¹, indicating that there is not much turbulent motion in this region. The SW halo nevertheless presents the higher value going up to 55 km s⁻¹ in most regions where the velocity is zero (see the regions with white velocity according to the colorbar in Fig. 4 and the zones with brown – white velocity dispersion in Fig. 5).

The corresponding one-component fit of the spectrum shown in Fig. 2 is represented in Fig. 6. As we can see, this simple fit does not capture all the complexities inherent in the observed line profile, which exhibits broadening in the base. Regions where the single-component fit shows high velocity dispersions are often indicative of a more complex or broad Hα line profile.

Fig. 6. One-component fit of a Hα profile of a single spaxel in the SW of the MB of the galaxy (same spaxel as in Fig. 2). The profile is depicted in black, and the one-component fit in yellow.

3.2.2. Double component

In contrast, a double-component Gaussian fit, applied to the same spectrum, provides a more nuanced and accurate representation (see Fig. 7). The small, broad Gaussian component successfully fits the broadening in the base of the profile, thereby resolving the complexity in the profile. This dual-component approach proves particularly valuable in regions where the single-component fit is inadequate, as it disentangles the overlapping velocity structures and offers a clearer understanding of the internal dynamics of IZw18.

Fig. 7. Double-component fit of a Hα profile of a single spaxel in the SW of the MB of the galaxy (same spaxel as in Fig. 2). The profile is depicted in black, and the double-component fit in yellow. The green lines represent each of the Gaussians corresponding to the double-component fit.

Expanding our discussion to encompass the double-component Gaussian fits, we simultaneously examined both the principal and secondary components. In doing so, we aimed to provide a more comprehensive understanding of the intricate dynamics at play within the different regions of the galaxy.

We initiate our discussion with the Hα luminosity maps, depicted in Figs. 8 and 9 for the principal and secondary components, respectively. It is noteworthy that the sum of the luminosities from the principal and secondary components in each spaxel closely approximates the luminosity captured by the single-component fits. This indicates a high degree of correspondence between the two approaches, yet with the double-component fit offering a more nuanced view. The principal component is typically more luminous compared to the secondary component (around four times more in the MB and two times more in the NE halo, though about the same in the SW halo), fitting well with the two knots in the MB of the galaxy. Meanwhile, the secondary component appears to trace the plume of the galaxy.

Fig. 8. Hα luminosity maps for the MB and halos corresponding to the principal double-component fit.

Fig. 9. Hα luminosity maps for the MB and halos corresponding to the secondary double-component fit.

Transitioning to the Hα velocity, the maps presented in Figs. 10 and 11 elucidate the internal kinematics of the ionized gas in IZw18. Remarkably, the velocity map of the principal component bears a close resemblance to that of the single-component fit, both exhibiting patterns indicative of rotation. This similarity is especially pronounced in regions where the single-component velocity dispersion is low (< 25 km s⁻¹), suggesting that the Hα profile is less complex and that the single-component fit is sufficiently accurate in those areas. However, this is not the case for the SW halo, where the velocity of the principal component divides this region into two distinct zones, north and south, with markedly different velocities. The transition between these zones is abrupt, indicating a significant variation in the kinematic properties of the SW halo.

Fig. 10. Hα velocity maps for the MB and halos corresponding to the principal double-component fit.

Fig. 11. Hα velocity maps for the MB and halos corresponding to the secondary double-component fit.

The secondary component presents a significantly different picture. Its velocity range is more extensive than that of the principal component (principal component goes from −40 km s⁻¹ to 55 km s⁻¹, secondary instead goes from −55 km s⁻¹ to 80 km s⁻¹), possibly pointing to higher velocities in subtler areas of the galaxy. The shape and distribution of velocities in this secondary component are also more challenging to interpret. While the principal component largely aligns with previously studied HI gas rotations (e.g., Lelli et al. 2012), the secondary component could either be tracing different gases in different regions or potentially be indicative of large-scale turbulent motions within the galaxy.

Lastly, we explored the Hα velocity dispersion maps for both components, as shown in Figs. 12 and 13. The velocity dispersion in the principal component closely resembles that of the single-component fit, particularly in regions where both are low (< 25 km s⁻¹). Notably, in areas where the single-component dispersion is high (> 25 km s⁻¹), the principal component often shows a lower dispersion value, hinting at the limitations of the single-component approach in capturing the full dynamical complexity.

Fig. 12. Hα velocity dispersion maps for the MB and halos corresponding to the principal double-component fit.

Fig. 13. Hα velocity dispersion maps for the MB and halos corresponding to the secondary double-component fit.

The secondary component, on the other hand, presents a markedly different scenario. Its velocity dispersion values are up to more than five times greater than those retrieved from either the single or principal components, underscoring the presence of high turbulence and or kinematic complexity in more nuanced regions of the galaxy. Most strikingly, the eastern (and with less relevance also the western) part of the MB exhibits velocity dispersions reaching up to 225 km s⁻¹, potentially tracing underlying outflows in this particular zone of IZw18. A 3D representation of all these kinematic features of the observed regions of IZw18 can be seen in Appendix A.

3.3. Luminosity-weighted velocity

In this section we delve into the analysis of the kinematics of IZw18 through luminosity-weighted histograms for velocity and velocity dispersion, based on the maps discussed in Sect. 3.2. The core of this approach lies in weighting the velocity and velocity dispersion values by the luminosity associated with each spaxel, which is derived from the Hα emission line. This luminosity is directly proportional to the number of ionized hydrogen atoms, and thus to the mass of ionized gas, providing a mass-weighted perspective of the kinematic properties. By employing this method, we aimed to emphasize the contributions from regions with higher mass concentrations, offering a more nuanced understanding of the dynamic behavior of the galaxy.

To construct the velocity histograms, we started by establishing the range of velocities, defined by the maximum and minimum values observed in the velocity map. This range was then divided into uniformly sized bins. For each bin, we aggregated the luminosities of all spaxels whose velocities fall within the range of the bin, effectively creating a sum of luminosities that represent the mass of ionized gas moving at those velocities. Normalizing these sums by the total luminosity of all spaxels yields the luminosity-weighted velocity histogram. A peak in this histogram indicates a predominant velocity at which a significant mass of ionized gas is moving, providing insight into the dominant kinematic flows of the galaxy.

Similarly, the velocity dispersion histograms are constructed using the same luminosity-weighting approach, allowing us to assess the spread of velocities around the mean motion in a mass-weighted manner. This method enriches our understanding of the internal motions within the ionized gas, revealing the complexities of its dynamics and the forces at play.

Starting with the MB of IZw18 in Fig. 14, we present the constructed histograms, elucidating the distributions of velocity and velocity dispersion, each weighted by luminosity. The top left histogram in Fig. 14 delineates the luminosity-weighted velocity histogram, highlighting the predominant velocities within the ionized gas. A noteworthy observation from this histogram is the resemblance in shape between the single-component fit and the principal double-component fit, with both showcasing peaks at around −5 km s⁻¹, a velocity that aligns with the brightest part of the MB of IZw18, the northern knot. These components also present a secondary peak at around 30 km s⁻¹, which corresponds with the southern knot. However, the secondary double-component fit displays a more dispersed pattern in velocity, extending the breadth of the observed velocities.

Fig. 14. Histograms of the luminosity-weighted velocity and velocity dispersion. Top row: MB. Middle row: NE halo. Bottom row: SW halo. Left column: Luminosity-weighted velocity. Right column: Luminosity-weighted velocity dispersion. The black color represents calculations made with the single-component fit, and the blue and red colors represent the principal and secondary double-component fit, respectively. Vertical solid thick lines represent the median of each distribution. The width of each cell on the x axis is 2 km s−1

Conversely, the top right histogram in Fig. 14 illustrates the luminosity-weighted velocity dispersion histogram, shedding light on the prevailing velocity dispersions within the galaxy. Similar to the velocity histogram, the single-component fit in the velocity dispersion histogram bears a resemblance to the principal double-component fit in shape. Nevertheless, there are stark differences in the values they reach. The principal double-component fit achieves lower velocity dispersions, with its median situated around 15 km s⁻¹, whereas the median of the single-component fit is approximately 20 km s⁻¹. The secondary double-component fit, on the other hand, reveals a markedly higher velocity dispersion, with a median of 50 km s⁻¹, that is much more spread that the other components.

The middle panels in Fig. 14 show the histograms of the NE halo of IZw18, offering insights into the luminosity-weighted distributions of velocity and velocity dispersion. The middle left histogram in Fig. 14 distinctly outlines the luminosity-weighted velocity histogram, emphasizing the prevailing velocities within the ionized gas of this region. On the other hand, the middle right histogram in Fig. 14 portrays the luminosity-weighted velocity dispersion histogram, revealing the dominant velocity dispersions within the NE halo of the galaxy.

Upon close examination of these histograms, we observe stark kinematic distinctions between the principal and secondary components of the double-component fit. The principal component is distinctly shifted to bluer, or lower, velocities, exhibiting a median at around −30 km s⁻¹ and presenting a relatively low velocity dispersion with a median around 13 km s⁻¹. In contrast, the secondary component is skewed toward redder velocities, with a median at around 3 km s⁻¹, and manifests a higher velocity dispersion, characterized by a median around 32 km s⁻¹.

This pronounced separation in the kinematics of the different components is corroborated by the velocity and velocity dispersion maps of the NE halo (see Figs. 10–13), which vividly depict the disparate kinematic patterns. Interestingly, the single-component fit occupies a middle ground between these two configurations, serving as a poignant illustration of the limitations inherent in single-component fits in capturing the full spectrum of kinematic information. Indeed, even in the regions of the galaxy that are relatively kinematicly quiescent (i.e., those presenting lower velocity differences and reduced velocity dispersion), significant kinematic diversity is observable, underscoring the necessity of employing nuanced analytical approaches to unravel the multifaceted kinematic landscape of the galaxy.

Similarly, for the SW halo of IZw18, our focus shifts to the bottom histograms in Fig. 14. These figures unveil the developed histograms, demonstrating the distributions of velocity and velocity dispersion, each weighted by luminosity, within this segment of the galaxy. Bottom left histogram in Fig. 14 offers an in-depth perspective on the luminosity-weighted velocity histogram, highlighting the prominent velocities within the ionized gas of the SW halo. Conversely, bottom right histogram in Fig. 14 reveals insights into the luminosity-weighted velocity dispersion histogram, illustrating the predominant velocity dispersions in this region.

A salient observation here is the comparable luminosities between the principal and secondary components of the double-component fit, rendering the distinction between a principal and a secondary component somewhat ambiguous in this instance. In terms of velocity, the principal component is bifurcated into two regions, around −15 km s⁻¹ and 45 km s⁻¹, with the latter leaning toward the redder, or higher, velocities. The secondary component has its median also around −15 km s⁻¹ but extends to bluer, or lower velocities, distinguishing it as the bluer component this time. The single-component fit situates itself between these two components, without reaching the extreme velocities exhibited by them.

In relation to velocity dispersion, the observations here parallel those in other regions of the galaxy, with the principal component exhibiting lower and the secondary component displaying higher and more varied velocity dispersion. However, it is noteworthy that the shape of the single component in this region aligns more with the secondary component than with the principal one.

An additional noteworthy point is the heightened kinematic excitement in the SW halo compared to the NE halo. This added layer of kinetic complexity likely intensifies the intricacy of the system, making the interpretations more challenging. This nuance underscores the importance of meticulous analysis in decoding the diverse and intricate kinematic patterns inherent in different regions of the galaxy.

This comprehensive approach, encompassing different regions of the galaxy, furnishes a holistic understanding of the disparate kinematic behaviors and patterns pervading throughout the various segments of IZw18.

4. Integrated spectra

We also took advantage of our IFS data to produce the 1D spectra of selected galaxy regions. These integrated spectra allow us to study each galactic region (MB knots, NE halo, and SW halo) as a whole, while at the same time improving the S/N.

While the integrated spectra (for each region) could be constructed using all available spaxels, a more discerning approach is employed to optimize the quality of the spectra obtained. Rather than indiscriminately incorporating every spaxel, we strategically added one spaxel at a time, starting with the most luminous ones, to strike a balance between enhancing signal strength and minimizing noise interference. This method allowed us to identify a point where the S/N is optimized. We halted the addition of spaxels at this point to avoid diluting the quality of the spectra with excessive noise from less luminous spaxels. In essence, this meticulous approach ensures the resultant integrated spectra are both robust and representative of the true characteristics of the selected galactic regions. As a way to demonstrate this approach in Fig. 15 is represented the S/N versus the number of spaxels in the MB of IZw18. The signal (S) is the Hα luminosity, meanwhile the noise (N) is the standard deviation of the whole continuum times the full width at zero intensity (FWZI). This allows us to compare the area of the profile with an area defined by the noise, which makes S/N to be dimensionless. We can see that the maximum S/N is reached using the 85 brightest spaxels in this case. Furthermore, in Fig. 16 we present the optimal spectrum of the IZw18 MB along with the integrated spectra retrieved from taking all spaxels into account. We can clearly see that the optimized spectrum shows a lower level of noise than the other one.

Fig. 15. S/N vs. number of spaxels in the MB of IZw18.

In the integrated spectra of the MB of IZw18, we observe a rich variety of spectral lines beyond the Hα emission. Notably, the spectrum reveals the presence of the [NII]λλ6548,6584 doublet and the [SII]λλ6716,6731 doublet, as well as [OI]λ6300, [SIII]λ6312, recombination lines from HeI, and [ArIII]λ7135. These additional spectral features are evident in Fig. 16. In contrast, the spectral analysis of both halos of IZw18 presents a markedly different scenario. Here, we observe only the Hα line, with no discernible presence of other spectral lines. This is attributed to the faintness of the emission in these regions.

Fig. 16. Integrated spectrum of the MB. The spectrum that optimizes the S/N is shown in black, and the integrated spectrum where all spaxels were used in blue.

To provide a visual representation of our selective approach in constructing integrated spectra, in Fig. 17 we present the Hα maps of the observed galactic regions showing only the spaxels that were used to create the integrated spectra. Each panel within the figure represents only the spaxels that were utilized to make the optimal integrated spectra. As we can see in this image, the MB is divided into two disconnected regions, the northern knot and the southern knot. In this way, we created one integrated spectrum for each knot in order to distinguish their kinematic behavior. We can compare this image with the one that shows all spaxels (see Fig. 3). Noticeable, only a small area comprising the brightest spaxels of each region is used in order to construct the integrated spectra.

Fig. 17. Luminosity map of the single-component fit. The spaxels represented here are those we used to form the integrated spectra in each region.

To delve deeper into the kinematic properties of each galactic region, we present a detailed view of the Hα line of the integrated spectra in Figs. 18–21. In each region, a multi-Gaussian fit is performed on the Hα line to extract the kinematic properties that represent these regions as a whole, utilizing the optimal S/N provided by our data. The AIC is employed to determine the optimal number of Gaussians for the fit, ensuring a balance between the goodness of fit and the complexity of the model. This approach allows us to accurately model the Hα profile of each region and gain insights into the underlying kinematic structures with enhanced reliability. To quantify the uncertainties associated with the Gaussians in the fit, we employed the bootstrapping method (Efron & Tibshirani 1985). Table 2 shows the luminosity percentage, velocity, velocity dispersion, and FWZI (along with their corresponding errors) of the components that fit each region.

Fig. 18. Zoomed-in view of the Hα line in the integrated spectrum of the northern knot of the MB. The observed spectrum is shown in black and the multi-Gaussian fit in yellow. In green we show the continuum plus the individual Gaussians. The Gaussian with a high standard deviation is shown in red (σ > 10 Å). This same color scheme is used in the representations of integrated spectra shown in the subsequent figures.

Fig. 19. Zoomed-in view of the Hα line in the integrated spectrum of the southern knot of the MB.

Fig. 20. Zoomed-in view of the Hα line in the integrated spectrum of the NE halo.

Fig. 21. Zoomed-in view of the Hα line in the integrated spectrum of the SW halo.

In examining the Hα profile of the north (south) knot of the MB (see Figs. 18 and 19), the AIC suggests a fit with four (three) Gaussian components for the Hα line, alongside two additional Gaussians for the [NII] lines. The four (three) Gaussians comprise two (one) principal components, which account for ∼90% (in both cases) of the total luminosity and exhibit low velocity dispersion (from 6 km s⁻¹ to 21 km s⁻¹), and two additional components, contributing around 7% and 4% to the luminosity (in both cases), characterized by velocity dispersions of 65 km s⁻¹ (56 km s⁻¹) and an exceedingly broad velocity dispersion of 730 km s⁻¹ (840 km s⁻¹), respectively (see Table 2). Furthermore, the FWZI for the northern and southern knots is measured to be 3000 km s⁻¹ and 2600 km s⁻¹, respectively, underscoring the extensive breadth of the Hα profile, particularly in the context of the very broad components.

The emergence of these broad and very broad components introduces a complex layer to our understanding of the kinematic structure of the galaxy. Notably, the very broad components not only are found shifted to the red in comparison to the narrow components by at least ∼40 km s⁻¹, but they also exhibit a velocity spatial gradient that is almost three times that of the narrow components. Specifically, the velocity difference between the southern and northern knot for the very broad components is approximately 110 km s⁻¹, which, after adjusting for the 40 km s⁻¹ attributed to the rotation of the galaxy, results in a net difference of 70 km s⁻¹. This finding indicates that the very broad components are not only redshifted relative to the narrow ones but also suggest that the dynamics of the very broad components are significantly distinct, highlighting that a fraction of the ionized gas in the MB of IZw18 possesses remarkably high kinetic energy, potentially due to dynamic processes within the galaxy.

The likely physical mechanisms that could impart such a significant input of kinetic energy into the gas include supernovae explosions and/or black hole accretion. Supernovae, with their colossal energy output, can create shock waves that propagate through the interstellar medium, thereby inducing turbulent motion in the gas. Similarly, if there is a black hole present in IZw18, its accretion activity and associated energetic outputs (such as jets or outflows) could also serve as a potent source of kinetic energy, stirring the surrounding gas into a frenzied, disordered state. Even broader components (FWHM > 10⁴ km s⁻¹) are seen in a similar (very metal-poor dwarf) galaxy called SBS 0335-052E in Hatano et al. (2023), where they argued about the possibility of an active galactic nuclei being the source behind this broad emission. The discernment of this very broad component not only challenges our previous kinematic analyses but also beckons a deeper exploration into the internal dynamics and potential energy sources within IZw18, offering a rich avenue for further study and discussion.

The absence of discernible indicators of these broad components in earlier stages and the fact that so many components are found in the integrated spectra of the MB of IZw18 indicate that the dual-component fit in the previous 2D-kinematic analysis could be conflating these (up to four) components, especially in the secondary one. This revelation necessitates a reevaluation of the kinematic structure and the potential interplay between these components. Currently, our 2D-kinematic study extends up to the dual-component fit, emphasizing the two predominant (brightest) components discerned in the fit of the integrated spectra. The newfound complexity and multiplicity of components in the integrated spectral data invite further, more granular exploration into the kinematic properties of IZw18.

In the context of the NE halo of IZw18, the integrated spectrum discernibly align with a two-component fit (see Fig. 20). Our comprehensive 2D-kinematic examination of this region corroborates this dual-component structure, revealing distinct kinematics in both the maps representations (see Figs. 10–13). The two components manifest as a brighter, blue-shifted element with lower velocity dispersion and a dimmer, red-shifted counterpart exhibiting higher velocity dispersion. This dual-structure is consistently evident in the integrated spectrum fit (see Fig. 20), where both components are distinctly recognizable. The NE halo stands out as a region of relative kinematic simplicity, featuring two discernible components that can be readily attributed to distinct ionized gases exhibiting differing kinematics. This scenario exemplifies a case where our mathematical approach of Gaussian fitting converges with direct physical interpretations, offering insights into the movements of the ionized gas within this region.

Observations of the SW halo of IZw18 reveal a consistent two-component fit within the integrated spectra for this region (see Fig. 21), mirroring the dual-structure found in the NE halo. However, the components here exhibit an inverse configuration compared to those in the NE halo. Specifically, the brighter component is red-shifted and exhibits a wider spread, contrasting with the narrower, blueshifted dimmer component. The components in this region tend to have higher velocity dispersion, indicative of a broader kinetic range, reflecting the inherent kinematic complexity of the SW halo already seen in Sect. 3.2.

5. Very broad component

Based on the previous discovery of the very broad component present in the knots of the MB and the strikingly different velocities of this very broad component with respect to the narrow ones, we took a similar approach to the one presented in Sect. 3.2. This approach enabled us to generate luminosity, velocity, and velocity dispersion maps specifically for this very broad component.

Nevertheless, as we have noted before (Fig. 2), the very broad component cannot be seen in individual spaxels. One way to solve this is to perform a harder Gaussian spatial smoothing on the spaxels (though at the cost of losing further spatial resolution). Doing so, we were able to reach a point at which we have enough signal to distinguish this dim component from the continuum. This time we used a Gaussian kernel with a FWHM of 2.5″, which makes the area of the FWHM of the Gaussian kernel three times larger than the previous one used.

Furthermore, to properly fit the very broad component, we had to get rid of the more luminous narrow Hα emission and the [NII] lines. For this, we took away the points in the spectra located less than 5 Å (2.5 Å) from the center of the Hα line ([NII] lines).

We then made a one-component Gaussian fit to the profiles that appeared following this procedure. The fit corresponding to the same spaxel as Fig. 2 is represented in Fig. 22. Here we can see that thanks to the enhanced smoothing of the new kernel the very broad component is clearly visible and can be well fitted.

Fig. 22. One-component fit of the defined broad Hα component. The profile comes from the same spaxel as used in Fig. 2. The profile is depicted in black, and the one-component fit in yellow.

The same criteria for excluding spaxels as the one applied in Sect. 3 is done here. In this case, the criterion that the S/N be below 3 removes most of the spaxels.

To effectively compare these new results with the one-component fit described in Sect. 3, we reapplied the one-component fitting procedure to the entire Hα line, employing the new Gaussian kernel. This approach reaffirms the previously presented results, but with enhanced smoothing.

The maps of the very broad component are represented in Fig. 23 and luminosity ratio (very broad component luminosity / full Hα luminosity) and velocity difference (very broad component velocity – Full Hα velocity) maps are in Fig. 24.

Fig. 23. Three maps of the very broad component, showing the luminosity (top), velocity (middle), and velocity dispersion (bottom). In all maps, the contours of the total Hα luminosity and a circle representing the FWHM of the seeing are overplotted.

Fig. 24. Luminosity ratio and velocity difference. Top: Luminosity ratio (i.e., the very broad component luminosity divided by the full Hα luminosity). Bottom: Velocity difference (i.e., the very broad component velocity minus the full Hα velocity). In both maps the contours of the total Hα luminosity and a circle representing the FWHM of the seeing are overplotted.

The maximum of emission of the very broad component is slightly to the south of the maximum of the total Hα emission (see the top map in Fig. 23). Nevertheless, the maximum luminosity ratio is in the eastern and western part of the MB (see the top map in Fig. 24), where we reach up to 35 % of the emission in the form of this very broad component. These regions are actually the same that presents a velocity dispersion of around 200 km s⁻¹ in the secondary component of the dual-fitting (see Fig. 13). In this way, as we mentioned before, this secondary component fits several components of the Hα line profile seen in the integrated spectra of the knots (see Figs. 18 and 19), making it difficult to interpret. In the southern knot we do not see any secondary maximum in luminosity regarding the very broad component, as it is seen in the luminosity of the whole Hα line (see the top map in Fig. 23). This is reflected in the luminosity ratio map where in this region its value is reach a minimum of 3%. The south part of the northern knot along with the previously mentioned eastern and western regions of the MB present the maximum luminosity ratio, indicating possible zones dominated by outflows. The wide extended spatial region that shows very broad component emission indicates that the reason behind this high kinematic input in the gas should not be local, as we could expect from localized supernovae or active galactic nuclei in the northern knot.

Furthermore, looking at the velocity and velocity dispersion maps (see the middle and bottom maps in Fig. 23), we can see the tremendous kinematic excitement of this gas. The velocity has a range of more than 160 km s⁻¹, which is more than three times the velocity range presented in the one-component fit of the MB. The velocity dispersion does not present much spatial variation ranging from 600 km s⁻¹ to more than 800 km s⁻¹ in the southern knot. The velocity difference between the very broad component and the whole Hα line (see bottom map in Fig. 24) shows that the very broad component is mostly shifted to the red with respect the narrow component. Only in the west of the northern knot does this difference in velocity go to zero, and it even switches to the blue. This high velocity difference between the narrow and the very broad component indicates that these gases maybe occupy a different region along the line of sight (the very broad component gas likely being behind the narrow component gas). A possible explanation of this behavior is the presence of high density gas near the origin of kinematic input. This high density gas could act like a wall that, to some extent, reflects back the momentum of the gas itself, producing the observed change in the velocity of the low density gas in an opposite direction. In this way, in the red velocities of the bottom map in Fig. 24, this theoretical wall would be in the front part of the galaxy with respect the line of sight (closer to us). On the other hand, the blue velocities in this same map indicates that this high density gas is in the back of the galaxy with respect to the line of sight (far away from us). This kinematic interaction of different gases within the interstellar medium is reminiscent of the turbulent mixing layer mechanism, where emission-line broadening occurs at the interface of cold disk gas and warmer outflow gas, as discussed by Hogarth et al. (2020). The observed shifts in velocities between broad and narrow components have been seen in other galaxies like GP J 0820 (Bosch et al. 2019). This might be explained by this reflection of momentum along the line of sight. In fact, this could be a result of the high-density gas acting as a barrier, influencing the dynamics of the gas flow and leading to the velocity shifts seen in the spectral lines.

6. Conclusions

This study presents a comprehensive analysis of the kinematic behavior of the ionized gas in IZw18, a galaxy known for its extremely low metallicity and close proximity. Our research delves into the kinematic complexities of IZw18, exploring its structure and dynamics through an in-depth examination of the Hα line profiles across different regions. By applying both single and double Gaussian component fittings to the spaxel data, we successfully generated luminosity, velocity, and velocity dispersion maps for each Gaussian component in the MB of the galaxy and two distinct regions of the halo (the NE and SW halos). Additionally, our methodology included creating integrated spectra for select regions of the galaxy, focusing on maximizing the S/N. Finally, we extended our analysis by constructing luminosity, velocity, and velocity dispersion maps of the very broad component (discovered in the integrated spectra), employing a larger Gaussian kernel to enhance the signal, thereby allowing for a more detailed observation of this elusive component.

In the MB, our analysis reveals a notable rotational pattern connecting the northern and southern knots. However, the secondary double-component fit uncovers a different and more complex kinematic pattern. This complexity is further highlighted in the integrated spectra of both knots, where up to four Gaussian components were required to adequately fit the Hα profile. Notably, one of these components exhibited a FWHM of almost 2000 km s⁻¹ (i.e., 4% of the total flux), suggesting the presence of a high-energy outflow. Particularly, this very broad component exhibits a significant spatial extent that not only challenges the notion of localized high-energy events, such as supernovae or active galactic nucleus activity, but also questions the likelihood of bipolar outflows. Such outflows are generally expected in scenarios involving localized sources of kinetic energy. The widespread nature of this component instead suggests more global, galaxy-scale processes at work, indicating a complex interplay of dynamics that extends beyond traditional models of localized energetic events.

The halo regions of the galaxy reflect the behavior of the shells in IZw18. The NE halo exhibits a relatively tranquil kinematic state, with two Gaussian components in the integrated spectra showing a low velocity dispersion (≤20 km s⁻¹). These components align with the two components of the double-component fit, presenting bluer velocities than the MB. In contrast, the SW halo is characterized by a more complex kinematic behavior than the NE halo. This is reflected in the presence of higher velocities and velocity dispersions and more complex kinematic patterns in the maps.

This research aligns with and expands upon the existing body of work on the kinematics of dwarf galaxies, particularly those with characteristics similar to IZw18, a BCD galaxy. Consistent with observations of other galaxies of its class, our study confirms that IZw18 exhibits chaotic kinematics along with a discernible disk rotation pattern. This finding suggests that while a primary rotational pattern encompassing the knots of the MB is evident, secondary and more intricate mechanisms are also at play. These mechanisms, possibly stemming from stellar feedback, introduce additional kinetic energy into the interstellar medium, resulting in a more chaotic and unpredictable motion of the ionized gas of the galaxy.

Furthermore, a critical aspect of our findings revolves around the significant velocity differences observed between the narrow and the very broad components of the Hα line. This disparity hints at a distinct spatial segregation of these gases within the line of sight, with the very broad component likely residing behind the narrow component. Intriguingly, this spatial arrangement might be influenced by the presence of high-density gas near the origin of the kinematic input. This dense gas could function to some extent as a wall, reflecting back the momentum of the surrounding gas. Such a reflection would manifest as a change in the velocity of the center of mass of the gas, directed opposite to this “wall” of dense gas. In this way, the redshifted velocities in the bottom map of Fig. 24 suggest that this dense gas is located toward the front of the galaxy, relative to the line of sight (i.e., closer to us). Conversely, the blueshifted velocities imply that the high-density gas lies toward the back of the galaxy (farther from us). These observations provide a novel interpretation of the internal gas dynamics in IZw18, pointing to complex, multilayered kinematic interactions within the galaxy and challenging conventional understandings of gas motion and distribution in such environments.

The study of IZw18, with its notably low metallicity, offers a crucial window into early galaxy formation and evolution in the Universe. The metal-poor environment in IZw18 closely resembles the conditions prevalent in the first galaxies, providing a unique opportunity to observe and infer behaviors and processes that were likely common in the high-redshift Universe. Our proximity to IZw18, coupled with the high spatial resolution of our observations, has allowed us to discern its kinematic complexity in unprecedented detail. This complexity, likely a characteristic of many galaxies, becomes vividly apparent in IZw18 due to our advantageous observational position. However, the intricate nature of these kinematic patterns presents significant interpretation challenges, underscoring the sophisticated dynamical processes at play in such early-stage galactic environments. Through IZw18, we gain insights into the nascent stages of galaxy evolution, observing firsthand the intricate dance of forces shaping galaxies at the dawn of the cosmos.

Acknowledgments

We thank the referee for their constructive comments. Author Antonio Arroyo Polonio acknowledges financial support from the grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033. RA acknowledges support from ANID FONDECYT Regular Grant 1202007. The authors acknowledge the plan PID2021-123417OB-I00.

References

Appendix A: 3D kinematic representation

Fig. A.1. 3D representation of the kinematic data of IZw18. Top row: MB. Middle row: NE halo. Bottom row: SW halo. Left column: Data from the single-component fit. Middle column: Data from the principal double-component fit. Right column: Data from the secondary double-component fit. Each 3D figure presents the following features: The x and y axis correspond to the position of the galaxy in arcseconds. The z axis (vertical axis) represents the normalized luminosity. The colorbar corresponds to the measured velocity, and the size of each point is related to the velocity dispersion (as indicated in the top-left legend in each 3D figure).

All Tables

All Figures

	Fig. 1. High contrast composite image of IZw18 using HST ACS R, I, and V bands, shown in the red, green, and blue channels, respectively. Isocontours delineate the central knots (adapted from Papaderos & Östlin 2012). The three LCB pointings are shown. North is up, and east is to the left here and in all subsequent figures.
In the text

	Fig. 2. Hα profile of a single spaxel. The black line represents the profile. The gray line is a Gaussian with the instrumental FWHM. The vertical line highlights the Hα wavelength.
In the text

	Fig. 3. Hα luminosity maps for the MB and halos obtained from the single-component fit. The background layer, presented in a cyclical gray scale, represents the wide-field deep image from the HST using the f814w and f606w filters of the WFC2 instrument. This background is used in all the images of IZw18 in which the MB and both halos are represented.
In the text

	Fig. 4. Hα velocity maps for the MB and halos obtained from the single-component fit. The contours correspond to the Hα luminosity in log scale obtained from the same fit component. The same contours are represented in all Hα velocity and velocity dispersion maps.
In the text

	Fig. 5. Hα velocity dispersion maps for the MB and halos obtained from the single-component fit.
In the text

	Fig. 6. One-component fit of a Hα profile of a single spaxel in the SW of the MB of the galaxy (same spaxel as in Fig. 2). The profile is depicted in black, and the one-component fit in yellow.
In the text

	Fig. 7. Double-component fit of a Hα profile of a single spaxel in the SW of the MB of the galaxy (same spaxel as in Fig. 2). The profile is depicted in black, and the double-component fit in yellow. The green lines represent each of the Gaussians corresponding to the double-component fit.
In the text

	Fig. 8. Hα luminosity maps for the MB and halos corresponding to the principal double-component fit.
In the text

	Fig. 9. Hα luminosity maps for the MB and halos corresponding to the secondary double-component fit.
In the text

	Fig. 10. Hα velocity maps for the MB and halos corresponding to the principal double-component fit.
In the text

	Fig. 11. Hα velocity maps for the MB and halos corresponding to the secondary double-component fit.
In the text

	Fig. 12. Hα velocity dispersion maps for the MB and halos corresponding to the principal double-component fit.
In the text

	Fig. 13. Hα velocity dispersion maps for the MB and halos corresponding to the secondary double-component fit.
In the text

Fig. 14. Histograms of the luminosity-weighted velocity and velocity dispersion. Top row: MB. Middle row: NE halo. Bottom row: SW halo. Left column: Luminosity-weighted velocity. Right column: Luminosity-weighted velocity dispersion. The black color represents calculations made with the single-component fit, and the blue and red colors represent the principal and secondary double-component fit, respectively. Vertical solid thick lines represent the median of each distribution. The width of each cell on the x axis is 2 km s−1

In the text

	Fig. 15. S/N vs. number of spaxels in the MB of IZw18.
In the text

	Fig. 16. Integrated spectrum of the MB. The spectrum that optimizes the S/N is shown in black, and the integrated spectrum where all spaxels were used in blue.
In the text

	Fig. 17. Luminosity map of the single-component fit. The spaxels represented here are those we used to form the integrated spectra in each region.
In the text

Fig. 18. Zoomed-in view of the Hα line in the integrated spectrum of the northern knot of the MB. The observed spectrum is shown in black and the multi-Gaussian fit in yellow. In green we show the continuum plus the individual Gaussians. The Gaussian with a high standard deviation is shown in red (σ > 10 Å). This same color scheme is used in the representations of integrated spectra shown in the subsequent figures.

In the text

	Fig. 19. Zoomed-in view of the Hα line in the integrated spectrum of the southern knot of the MB.
In the text

	Fig. 20. Zoomed-in view of the Hα line in the integrated spectrum of the NE halo.
In the text

	Fig. 21. Zoomed-in view of the Hα line in the integrated spectrum of the SW halo.
In the text

	Fig. 22. One-component fit of the defined broad Hα component. The profile comes from the same spaxel as used in Fig. 2. The profile is depicted in black, and the one-component fit in yellow.
In the text

	Fig. 23. Three maps of the very broad component, showing the luminosity (top), velocity (middle), and velocity dispersion (bottom). In all maps, the contours of the total Hα luminosity and a circle representing the FWHM of the seeing are overplotted.
In the text

	Fig. 24. Luminosity ratio and velocity difference. Top: Luminosity ratio (i.e., the very broad component luminosity divided by the full Hα luminosity). Bottom: Velocity difference (i.e., the very broad component velocity minus the full Hα velocity). In both maps the contours of the total Hα luminosity and a circle representing the FWHM of the seeing are overplotted.
In the text

Fig. A.1. 3D representation of the kinematic data of IZw18. Top row: MB. Middle row: NE halo. Bottom row: SW halo. Left column: Data from the single-component fit. Middle column: Data from the principal double-component fit. Right column: Data from the secondary double-component fit. Each 3D figure presents the following features: The x and y axis correspond to the position of the galaxy in arcseconds. The z axis (vertical axis) represents the normalized luminosity. The colorbar corresponds to the measured velocity, and the size of each point is related to the velocity dispersion (as indicated in the top-left legend in each 3D figure).

In the text