Spatially-resolved chemodynamics of the starburst dwarf galaxy CGCG 007-025: Evidence for recent accretion of metal-poor gas

Nearby metal-poor starburst dwarf galaxies present a unique opportunity to probe the physics of high-density star formation with a detail and sensitivity unmatched by any observation of the high-z Universe. Here we present the first results from a chemodynamical study of the nearby, gas-rich starburst dwarf CGCG 007-025. We use VLT/MUSE integral field spectroscopy to characterise the properties of the star-forming (SF) gas, from its metal content to its kinematics. The star formation rate (SFR) surface density presents a clumpy distribution, with the brightest knot hosting a 5 Myr young, Wolf-Rayet (WR) population (revealed by the presence of the characteristic 5808\AA~WR bump). The ionised gas kinematics are dominated by disordered motions. A superposition of a narrow ($\sigma \approx$ 50 km s$^{-1}$), intermediate (150 km s$^{-1}$) and broad (1000 km s$^{-1}$) kinematic components are needed to model the emission line profiles in the brightest SF region, suggesting the presence of energetic outflows from massive stars. The gas-phase metallicity of the galaxy spans 0.6 dex and displays a strong anti-correlation with SFR surface density, dropping to 12+log(O/H) = 7.7 in the central SF knot. The spatially-resolved BPTs indicates the gas is being ionised purely by SF processes. Finally, the anti-correlation between the SFR and the gas metallicity points out to accretion of metal-poor gas as the origin of the recent off-centre starburst, in which the infalling material ignites the SF episode.


INTRODUCTION
Observations (Maiolino & Mannucci 2019) and models of galactic chemical evolution (Ma et al. 2016) indicate that massive galaxies tend to have higher gas-phase metallicities than low mass systems. The so-called mass-metallicity relation (MZR, Tremonti et al. 2004) naturally emerges from the combination of secular and dynamical processes. First, the gas content in massive galaxies is in a more evolved stage and hence are chemically richer (Tinsley & Larson 1978;Edmunds 1990). Second, low mass galaxies, thanks to their lower gravitational potential, are unable to retain most of its metals, which are lost through galaxy winds and outflows (e.g., Chisholm et al. 2015;Tortora et al. 2022;Xu et al. 2022). The MZR presents substantial scatter -specially and more severely in the low mass regime (Curti et al. 2020)-which seems to be associated with the star formation rate (SFR, Mannucci et al. 2010; Lara-López et al. ★ E-mail: macarena.garciavalle@ed.ac.uk 2010), so that for galaxies with the same stellar mass, the less metallic systems will have higher SFR. The interplay between galaxy mass, SFR and metallicity is thought to arise as a result of the complex equilibrium between inflows, outflows, and mergers.
Theory and simulations predict that galaxies with high baryonic mass fractions ( gas > 0.3) have a turbulent interstellar medium (ISM, Hayward & Hopkins 2017), characterised by clumpy structures and whose star formation histories (SFHs) are regulated by inflows and outflow of material (Muratov et al. 2015). Although high gas mass fractions were common among the high-z galaxy population, in the local Universe this almost exclusively occurs in low mass systems (Bradford et al. 2015). Thus, nearby dwarf starburst galaxies are considered the best candidates when looking for highz analogues, sharing stellar masses, metallicities and specific SFRs with the average high-redshift galaxy (Weisz et al. 2011;Morales-Luis et al. 2011;Brammer et al. 2012;Garland et al. 2015;Sánchez Almeida et al. 2016;Izotov et al. 2021).
However, low mass galaxies reside in low mass haloes which are usually highly inefficient in capturing, retaining (Davé 2009;Chris-tensen et al. 2018) and converting baryons into stars (Geha et al. 2006;Kaufmann et al. 2007) due to their shallow potential wells. In addition, the absence of stabilising structures -like bars or spiral armsin dwarf galaxies reduce the effectiveness of the starbursts. Alternatively, mergers and interactions between gas-rich galaxies with comparable masses can induce strong starbursting episodes (Luo et al. 2014;Knapen & Cisternas 2015;Pearson et al. 2019). Simulations show that the tidal torques by interactions canalise large quantities of preexisting metal-poor gas from the outskirts towards the central regions (Muratov et al. 2015). Moreover, the interaction-induced high gas pressures not only translate into higher star formation efficiencies, but also facilitate the formation of star clusters and super star cluster complexes (see Adamo et al. 2020 for a review). Hence, interactions between dwarf galaxies can effectively spark and sustain star formation for longer periods of time (Patton et al. 2013;Moreno et al. 2021). Although dwarf-dwarf interactions must be extremely abundant in the early Universe (due to the hierarchical nature of structure formation), the vast majority of dwarf galaxies in the Local Universe are satellites or isolated galaxies, with less than 5% being part of a dwarf galaxy pair or a dwarf-only group (Stierwalt et al. 2015;Besla et al. 2018).
Cold accretion of gas from the intergalactic medium (IGM) is also predicted to produce star formation in dwarf galaxies, at high and low redshifts (e.g., Dekel & Birnboim 2006;Dekel et al. 2009). Direct detections of gas inflowing into a galaxy are relatively difficult to observe (Sánchez Almeida et al. 2014a). First, its absorption features are very weak, mainly because of its expected low density and low metallicity (e.g., Fumagalli et al. 2011). Second, its covering fraction may be very small (less than 10%, e.g., Fumagalli et al. 2011;Rubin et al. 2012), requiring an alignment with the line of sight (Ho et al. 2019). Various efforts have been made to detect the pristine gas accreting into galaxies by using absorption systems (e.g., Crighton et al. 2013;Zabl et al. 2019). However, metal-poor gas accretion leaves some imprints in the characteristics of the galaxy that can be more easily detected. Bright, off-centre starburst clumps are typically seen as a result of intense star formation episodes induced by gas accretion. These off-centre starburst clumps have been detected vastly in intermediate and high redshift galaxies Amorín et al. 2014Amorín et al. , 2015Guo et al. 2015;Calabrò et al. 2017). As good analogues of these high redshift star-forming galaxies (Papaderos et al. 2008;Papaderos & Östlin 2012), local dwarf galaxies also show various signs of star formation in their outskirts. For example, Sánchez  studied the bright star-forming regions of 7 tadpole galaxies (also known as cometary, Loose & Thuan 1986), a particular type of object formed by an elongated diffuse component with a bright clump in one of the extremities. They found these peripheral star forming regions have a deficit in metallicity, and associated its origin to a recent event of gas accretion (see also Sánchez Almeida et al. 2014b. In this vein, starburst dwarf galaxies (in which we can include blue compact dwarfs, BCDs; extremely metal poor galaxies, XMP; or HII galaxies) are the most promising local analogues for such chemically unevolved dwarf galaxies, since they tend not only to be metal-poor systems but also to have a high fraction of gas (HI mass fraction above 30%, Amorín et al. 2009Amorín et al. , 2016Filho et al. 2013Filho et al. , 2016Trebitsch et al. 2017) as well as strong stellar feedback (Muratov et al. 2015). Their strong bursts of star formation are usually on top of a much older stellar population, dominating the detected optical emission (Loose & Thuan 1986;Bergvall & Östlin 2002;Amorín et al. 2007Amorín et al. , 2009). Kinematical analyses of some BCDs indicate erratic velocity fields that point to dwarf galaxy mergers and/or in-falling of large gas clouds as the sources of their starburst (Östlin et al. 2001;Moiseev et al. 2015;Watts & Bekki 2016).
The development of Integral Field Spectroscopy (IFS) allowed not only to improve these studies but also to compare the properties of the SF regions with their surroundings. Cresci et al. (2010) reported inverted metallicity gradients in high redshift galaxies using near-IR IFS data: the gas-phase metallicity decreases towards the centre of the galaxies, coincident with the star forming regions. Similar findings are presented in Troncoso et al. (2014) for Lyman-break galaxies at ∼ 3.4. This suggest very metal poor gas was accreted into the centre of the galaxies. Using data from the MaNGA survey, Hwang et al. (2019) found some off-centred young star forming regions with an anomalously low-metallicity. They found these regions are more easily found in blue low-mass galaxies, and they concluded their origin was the accretion of metal-poor gas. Sánchez-Menguiano et al. (2019) found using 700 MaNGA galaxies (Bundy et al. 2015) that local enhancements in star formation coincide with drops in gas phase metallicity. Recently, the spatial metallicity variations (ΔO/H ≈ 0.7 dex) seen in Fernández-Arenas et al. (2022) were attributed to the accretion of pristine gas which fell into the system during merging stages.
In this paper we use MUSE IFS observations to study the ionised gas properties of the dwarf galaxy CGCG 007-025. The paper is organised as follows. In Section 2, we introduce the data and spectral fitting routine. In Section 3, the distribution and properties of the ionised hydrogen gas are presented. We obtain the chemical abundances and compute the ionisation maps in Section 4, in addition with the spatially resolved BPT diagrams (BPT, after Baldwin et al. 1981). Section 5 is devoted to the individually detected star forming knots. In Section 6, our main findings are discussed in the framework of the gas accretion scenario. In Section 7, we give a summary and present our conclusions.

Target description
CGCG 007-025, from the Catalogue of Galaxies and Clusters of Galaxies (Zwicky et al. 1968, other cross-IDs are MCG +00-25-010 andSDSS J094401.86-003832.1, De Propris et al. 2007;Ahumada et al. 2020), is known as a low metallicity starburst galaxy with extremely intense emission lines. At a distance of = 23 ± 5 Mpc (Kourkchi et al. 2020), CGCG 007-025 is a dwarf galaxy (M ★ = 1.2 × 10 8 M , Marasco et al. 2022) in interaction with another dwarf -UGC 5205-at a projected separation of 8.3 kpc (van Zee 2000). The interaction between them has triggered recent star formation episodes in both galaxies. This galaxy has been extensively studied as a prototypical nearby, metal-poor starburst (Izotov & Thuan 2004;van Zee & Haynes 2006;Senchyna et al. 2017Senchyna et al. , 2022Berg et al. 2022, among others), and as a local analogue to high-z low mass galaxies (Brammer et al. 2012). However, recent ground-based Subaru observations taken under the Hyper Suprime-Cam (HSC) Survey make more evident the presence of the underlying old component of CGCG 007-025 (Figure 1). In a forthcoming study (del Valle-Espinosa et al., in prep.) we present a detailed study of the host galaxy, which is characterised by old stellar populations and boxy isophotes. The photometric centre of the isophotes does not match the location of the brightest star forming region in the galaxy. In fact, as revealed by HST/WFC3 imaging and VLT/MUSE spectroscopy, the isophotes share a common centre at the location of a nuclear star cluster (NSC, Sánchez-Janssen et al. 2019 CGCG 007-025 is therefore best classified as an early-type, nucleated dwarf spheroidal galaxy exhibiting an off-centre starburst, very much akin to the well-known class of BCD galaxies with prominent underlying hosts (Amorín et al. 2007(Amorín et al. , 2009. Figure 1 shows the HSC giy-colour image, in which the underlying stellar component of the galaxy, as well as the NSC, can be observed. The over-plotted solidwhite contour represent the H emission measure in the available IFS data for this galaxy (see Sec 3.1 for the description). The dottedwhite contour corresponds to the isophote at = 26 mag arcsec −2 . In this paper, we aim to study the properties of the off-centre starburst episode of the galaxy, and we will present the characteristics of the underlying population in a companion paper.

Observational data
Our study is based on archival observations collected with the Multi Unit Spectroscopic Explorer (MUSE, Bacon et al. 2010) at the VLT (Very Large Telescope; ESO Paranal Observatory, Chile) 1 . Observations were carried out with the nominal Wide Field Mode (WFM-NOAO-N), corresponding to a field of view (FoV) of 1 × 1 with a spatial sampling of 0.2 . The spectral coverage ranges from 4800 to 9300 Å with a resolving power of ≈ 2600 around the H line. The observing conditions yield a FWHM effective spatial resolution of 0.7 . Data were already reduced using the MUSE pipeline (Weilbacher et al. 2012).
A few spectral pixels present NaN values around the peak of the H emission line. The location of these spaxels coincides with the brightest SF region in the galaxy. Upon inspection of the spectra associated to these NaN values we discovered that the 1 Program ID 0102.B0325 [O ] 5007Å/[O ] 4959Å flux ratio is below the theoretical emissivity ratio (≈ 2.97). We associate this behaviour to saturation and non-linearity of the CCD. In Section 2.4 we describe the procedure follow to account for these effects.

Methods
To extend the analysis of the data over the largest possible area voronoi tessellation was applied to the MUSE cube using the V B routine from Cappellari & Copin (2003). The tessellation wavelength definition was adapted according to the characteristic to study. For the detailed analysis of the H emission (Sec. 3) we only keep spaxels with / ≥ 3 in the 6560 -6566 Å rest-frame wavelength range and target a minimum / = 15 for the resulting tessellated bin. This allowed us to reach the faintest nebular emission in the galaxy. When studying the ionisation and chemistry of the object (Sec. 4) we performed the voronoi tessellation around the faintest set of lines we need for this purpose: the sulphur doublet [S ] 6717,6732Å. In this case, the previous / cuts were defined over the range 6710 -6740 Å.
The procedure for the fitting routine starts by determining the stellar continuum. We mask every emission line and model the stellar continuum using PXF (Cappellari & Emsellem 2004), which finds the best fit for each spectrum-bin considering a linear combination of spectral templates. We use the single stellar populations (SSPs) templates from the E-library (Vazdekis et al. 2016), which cover the full MUSE spectral range.
Once we model the stellar continuum we subtract from the original MUSE datacube, creating a pure-gas cube. We split the gas cube in three spectral regions: from 4750 to 5200 Å, gathering the H and  Å, which covers the emission of [S ] 6717,6732Å. We expect the nebular continuum to have little contribution to the overall line fluxes (Byler et al. 2017). Nevertheless, we add a one degree polynomial to model the continuum level in these three windows separately that should account for the smooth nebular continuum contribution. Due to the presence of broad wings in the base of the brightest emission lines, sometimes we use the sum of up to three gaussian profiles to model the overall emission line.
Some of the species only require a narrow component for the fit (such as the nitrogen doublet), while the others show more complex kinematic components. For these cases, we distinguish between lines that may require one additional kinematic component (hereafter intermediate component) and lines where even a third component is needed to fit the profile (the broad component). The latter case is only relevant for the H emission line, whereas an intermediate component is frequently necessary in H and in the oxygen and sulphur doublets. The fitting routine starts with narrow-only components and successively adds as many components as described above. We use the Bayesian Information Criterion (BIC) to decide how many components are required to model the data (Schwarz 1978). The model with the smallest BIC value is considered to provide the best description of the emission. Additionally, we only preserve kinematic components whose peak S/N is above 5. Lines in doublets share kinematic properties (i.e. red/blue-shifting and velocity dispersion), and their height ratios are fixed to the theoretical values when appropriate (3 for the oxygen doublet and 2.95 for the nitrogen doublet; Osterbrock & Ferland 2006). The fits to the H and H lines are fully independent, since they are fitted in different spectral windows and their ratios contain information on the obscuration by dust. Typical flux uncertainties for the narrow component are ∼ 5%, while for the intermediate and broad components are ∼ 20%. The velocity uncertainties are are 0.1%, 2.5% and 10% for the narrow, intermediate and broad components, respectively. Uncertainties in the velocity dispersion are 3.5% for the narrow component while for the intermediate and broad components are ∼ 20%. For the spaxels in which we require the presence of three kinematic components in the H line, the narrow component is typically always dominant (97% of the total line flux), with the intermediate and broad components representing a much smaller contribution (2% and 1%, respectively). Unless explicitly mentioned, the analyses in this study always refer to the narrow component. We correct the velocity dispersion from instrumental broadening using the derived values from Mentz et al. (2016), which correspond to ∼50 km/s for the H emission line.
Once all the lines are fitted and the number of components is defined, we correct all line fluxes of Milky Way extinction using the Cardelli attenuation law (Cardelli et al. 1989). At the location of CGCG 007-025, the adopted E(B-V) value for the Milky Way corresponds to 0.035 (Green et al. 2018).  Table A1 are fitted for a better constrain of the continuum and/or the additional kinematic components. Fernández et al. (2023) (hereafter F22) obtained the ( ) map of CGCG 007-025 using the Large Magellanic Cloud (LMC) reddening law by Gordon et al. (2003) with a = 3.1. When the signal-tonoise was high enough, they computed the ( ) coefficient using H and the available Paschen lines, whereas when the signal-to-noise decreases, they used H and H for the computation of the extinction. With this approach they can map ( ) in the brightest spaxels of the galaxy, where the H and [O ] 5007 lines were saturated. They found the central clump of CGCG 007-025 has a higher dust content than the rest of the galaxy, with a peak of extinction at ( − ) = 0.35 ± 0.04 and dropping to ( − ) = 0.15 ± 0.06. When moving to the outskirts, a mean value of ( − ) = 0.08 ± 0.08 was found. For a more detailed description of the method and inferred extinction we refer the reader to F22. For the spaxels in which the ( ) is known, we use the values derived in F22 to correct from intrinsic extinction. Otherwise, we corrected from extinction assuming the outskirts value ( − ) = 0.08. Finally (Osterbrock & Ferland 2006).

Star formation rate
One of the brightest emission lines widely used to estimate the star formation rate (SFR) in the literature is the H recombination line, which traces the ionised gas in the surroundings of massive stars. Since H is directly linked to the presence of massive stars younger than ∼10 Myr, the derived SFR represents the instantaneous value (Kennicutt 1998;Kennicutt & Evans 2012). Using a Kroupa initial mass function (Kroupa 2001) and the S 99 stellar models (Leitherer et al. 1999) at 1 , Murphy et al. (2011) updated the SFR(H ) calibration originally presented in Kennicutt (1998) . We convert the extinction-corrected H fluxes into luminosities, and derive the SFR surface density (Σ SFR ) by dividing the SFR by the area of each bin (in kpc 2 ) for all the spaxels. In the upper panel of Figure 3 we present the Σ SFR map, with the different contours representing the log(Σ SFR [M yr −1 kpc −2 ]) levels at [-2.5, -2, -1.5, -1, -0.5, 0]. Σ SFR derived from the H emission line can suffer from systematic uncertainties, mainly coming from the adopted dust correction and from variations in the predicted ionising flux of massive stars. However, Catalán-Torrecilla et al. (2015) studied the reliability of the Σ SFR obtained from the dust-corrected H surface brightness in 104 SF galaxies from the CALIFA survey. They compared this value with the hybrid Σ SFR derived using the FUV and 22 emission, finding a good agreement between them especially when We compute the global instantaneous SFR of the galaxy by integrating the Σ SFR distribution, which yields a value of log(SFR/(M yr −1 )) = −0.67 ± 0.09. Marasco  GALEX and the W4 (22 ) from WISE. They found the total SFR of the galaxy to be log(SFR/(M yr −1 )) = −0.64 ± 0.13. Berg et al. (2022) also compute the SFR of CGCG 007-025 by modelling the SED of the object using the FUV and NUV data from GALEX, and the ugriz photometry from SDSS, corrected from aperture effects. Their modelling recovers a SFR value of log(SFR/(M yr −1 )) = −0.78 +0.19 −0.16 . Therefore, our estimation computed using the total H luminosity is in good agreement with the values in the literature.

H gas mass
In order to calculate the total mass of ionised gas present in the galaxy we follow the procedure in Bik et al. (2018, see also Carniani et al. (2015), Olmo-García et al. (2017)). The mass of the ionised gas is directly proportional to the ratio of the H luminosity, H , and the electron density, , following: where is the atomic weight (assumed to be 1.0), H is the hydrogen mass, ℎ is Planck's constant, H is the H frequency and eff H is the coefficient for H in the Case B recombination. The latter coefficient depends on the electron temperature of H , so we need to derive first the electron density and temperature in order to apply the formula.
F22 inferred the electron density and two electron temperatures (for the high and low ionisation species) of individual spaxels where the [S ] 6312Å line was reliable enough to apply the direct method. Briefly, the [S ] 6717,6732Å ratio was used to obtain the electron densities, the ratio between [S ] 6312,9069Å for the low ionisation temperature, and the electron temperature of the high ionisation species was derived using the empirical relation from Hägele et al. (2006). The low ionisation species were found to have a temperature of (1.5 ± 0.1) × 10 4 K, which corresponds to a value of eff H = (8.0 ± 0.7) × 10 −14 cm 3 s −1 (Pequignot et al. 1991). The electron densities of the brightest star forming clumps in CGCG 007-025 were also obtained in F22. The density of the brightest star-forming knot is as high as = 380 ± 60 cm −3 , with lower densities in the other three knots reaching = 90 ± 50 cm −3 . We further assume that this value applies to the outer parts of the galaxy.
With the electron density values on hand as well as the H luminosity, we create the ionised gas mass surface density map. By integrating this map, we estimate the total amount of ionised gas is H = (1.3 ± 0.9) × 10 6 . A comparison of this value with the stellar mass estimated by Marasco et al. (2022) for the galaxy reveals the ionised gas mass is around 1% of the stellar mass.

EW map
The equivalent widths (EWs) of the strongest hydrogen recombination lines correlate closely with the age of the most recent burst (Dottori & Bica 1981;Stasińska & Leitherer 1996;Terlevich et al. 2004;Levesque & Leitherer 2013). On the one hand, the most massive and hottest stars produce a nebular component that will drop rapidly with time due to their short lifetimes. On the other hand, the underlying absorption will slowly decrease because of the significant contribution from low mass stars. By studying the ratio of these two features we can measure the ratio between the young and old population, giving an estimation of the age of the most recent burst.
Among the different hydrogen recombination lines, the EW of the H line is a proven age indicator up to ∼ 10 Myr independently of the star formation history adopted (Martín-Manjón et al. 2008). Although the EW of the hydrogen recombination lines is affected by the metallicity of the population, several authors (Schaerer & Vacca 1998;Levesque & Leitherer 2013;Eldridge et al. 2017, among others) have characterised its dependence using different stellar libraries and evolutionary synthesis codes. Thanks to these studies, the EW(H ) is a reliable proxy to determine the age of recent starburst.
We build the EW map using the narrow component of the H recombination line since the continuum around this line is better determined than in the case of H . As a sanity check we also compute the H EW map, which leads to the same distribution.  Figure 4 shows the maps for the flux (top row), velocity (middle row) and velocity dispersion (bottom row) of the three kinematic components of the H line: narrow (left column), intermediate (middle column) and broad (right column). The contours displayed in all the maps are the same Σ SFR levels as in Figure 3.

Kinematics
The flux distribution map of the narrow component expands 3.1 kpc in the NW-SE direction and 1.6 kpc in the NE-SW. When compared with the area encompassed by the isophote at = 26 mag arcsec −2 (Figure 1), the ionised gas covers around 45% of the galaxy. The ionised gas emission is concentrated in four main starforming regions with an additional diffuse extended component. As highlighted in Figure 1, the morphology of the ionised gas is significantly different in all four main star-forming regions. From the S-E region to the central clump the gas morphology changes from being noticeably dispersed (as are the stars) to being significantly more concentrated (with the stars still embedded). Moreover, this change in the gas morphology seems to be linked to the age of the regions: the S-E region has the lowest EW(H ) -the oldest SF region-while the central region has the largest one -the youngest SF region-. In this scenario, the arc has an intermediate value of EW(H ) when compared with these other two regions, with its stars more concentrated than in the S-E region and its gas forming a cavity around them. The intermediate component is present in the four main star forming knots, whereas the broad is just in the brightest one.
The velocity field of the narrow component -middle left panel of Figure 4-does not exhibit any coherent kinematic structure, with instrumental-limited velocity dispersion values almost constant all across the galaxy -lower left panel of Figure 4-. Our findings for the narrow component are in agreement with the maps present in Marasco et al. (2022). In the southwest of the galaxy, the velocity map seems to be chaotic, with small patches with opposite velocity values. Across these patches, the velocity dispersion is slightly higher than at the position of the main bursts. Another region with higher velocity dispersion is located at the east side of the gas distribution, almost coincident with the position of the NSC. We cannot detect a rotational profile for the gas component in GCGC 0007-025 due to its highly perturbed velocity field.
The velocity field of the intermediate component -middle central panel of Figure 4-presents, in the brightest SF clump, a velocity pattern which is compatible with a bipolar structure: a local velocity maximum and minimum are observed along the NE-SW direction, with relative velocities of 50 km s −1 . Moreover, the internal velocities of the gas in this region are higher than the others -note the higher velocity dispersion values in the lower central panel of Figure 4-.
On the central clump of CGCG 007-025 we uncover the presence of a broad component (right column of Figure 4). The gas associated with this component exhibits substantial internal motions, with velocity dispersions exceeding ≈ 1000 km s −1 (lower right panel of Figure 4). The velocity field (middle right panel of Figure 4) is compatible with a dusty expanding bubble, where its receding parts are not observed because of the high extinction of the region. We are not able to detect this broad component in the H emission line. This could be related to the addition of two effects: (1) the high extinction of the central clump, and (2)   A similar behaviour in low metallicity galaxies has being found in other galaxies (Izotov et al. 2010;Cann et al. 2020). This absence of broad emission in the forbidden lines could be explained if the broad lines are originated in a gas whose density is well above the critical density of the [O ], i.e. ∼ 10 6 cm −3 .

Chemistry
In F22 we not only computed the electron density and the electron temperatures (as briefly explained in Sec. 3.2) but also the ionic and total abundances of the gas for the spaxels in which [S ] 6312Å was reliably measured. However, and since the

Ionising mechanisms
The different excitation and ionisation mechanisms in the nebula are disentangled by studying different emission line ratios (Dopita & Sutherland 2003;Osterbrock & Ferland 2006). Among the most common ratios we have access to The [O ]/H distribution peaks in the positions where the SF knots are located, and decreases towards the edges of the gas distribution. This distribution is expected when a recent star formation event has occurred, since [O ] is generated by high ionising sources like O-type stars. On the other hand, the distribution of the low ionisation line ratios present the opposite behaviour, with the minimums located at the position of the SF knots. This behaviour is also expected when the ionising source produces large quantities of UV photons, which also is the case for O-type stars.
The structure of an ionisiation-bounded SF cloud can be described by two ionisation regimes: the innermost region is where highlyionised species reside, whereas the outer parts are dominated by low ionisation species. By studying the relative emission of two ions with different ionisation potential, we can disentangle the ionisation structure of our SF clumps. The right panel of Fig. 6

displays the [S ]/[O ] ratio, where [S ] is our low ionisation element and [O ] represents the highly ionised one. For an ionisation-bounded nebula the inner regions correspond to small values of [S ]/[O ]
, whereas the outer parts should reach values larger than one (Pellegrini et al. 2012).
The brightest SF knot in CGCG 007-025 has an average [S ]/[O ] value of 0.03, indicating the entire region is highly ionised. None of the 4 main star forming clumps present an [S ]/[O ] distribution expected for a classical ionised bounded nebula. This suggest the four SF clumps in the galaxy are optically thin, allowing the LyC photons to escape and ionise the surrounding gas (Pellegrini et al. 2012;Bik et al. 2015).
Our emission-line ratio maps allow us to explore spatially resolved diagnostic diagrams (the so called BPTs, Baldwin et al. 1981;Veilleux & Osterbrock 1987). By comparing two emission line ratios of species with different ionisation levels, BPTs manifest the origin of the emission lines since the position of the spaxels reflects the gas ionisation state and source of ionisation. Using a set of photoionisation models, Kewley et al. (2001) conclude that, if the emission is coming purely from starbursting episodes, the points should always fall below an empirical limit. On the other hand, if other ionisation mechanisms are involved in the generation of the emission lines -such as AGN or shocks-the points will move above this empirical limit. The position of this boundary in the BPT diagram depends mainly on the metallicity of the population and the ionisation parameter.
In Figure 7 Kewley et al. (2001) dividing lines are too conservative for a starburst galaxy with a metallicity as low as that for CGCG 007-025. In order to account for the metallicity deviation of CGCG 007-025 with respect to the solar value, we overplotted in both panels of Figure 7    grid used to derive the chemistry of the gas in Sec. 4.1. We limit the grid to the values obtained, so that log(N/O) = -1.5 and the oxygen abundance ranges from 7.6 to 8.4, while the ionisation parameter spans from -3.0 to -2.0. As shown in both diagnostic plots, the line ratios recovered for the narrow component can be entirely explained by star formation activity.

SPECTROSCOPY OF INDIVIDUAL STAR-FORMING CLUMPS
In F22, we identified the four main SF knots of CGCG 007-025. These regions were delimited at the 99.5 of percentile the flux distribution of all the spaxels in the datacube for the [S ] 6312Å emission line, the faintest line needed for the direct method. This percentile corresponds to a flux of 9.3 × 10 −20 erg s −1 cm −2 , implying all the spaxels have a flux of at least 1 × 10 −17 erg s −1 cm −2 for the H emission. We add the spaxels with flux values above this H flux limit to get the integrated spectrum of each SF clump. Figure 8 zooms into the H flux distribution of each individual clump. The names displayed on top of the zoom images are first adopted in F22. Note that the so-called central region is our brightest SF clump and the closest to the NSC, even though it does not overlap with the  NSC. With the data from Subaru/HSC it is evident that this region is not at the photometric centre of the galaxy, but we keep the name for consistency. The cyan contours mark the extent of the main SF clumps.
We analysed the spectra of the four clumps in the same way as explained in Sec. 2.3: we estimate the stellar continuum using PXF and fit the different emission lines using the same multi-component gaussian modelling. From the emission line modelling we obtain the EW(H ), the velocity, the total H luminosity and the SFR for each individual region. We also compute the integrated chemistry using HII-CHI-with the same photoionisation models as before. The derived characteristics of the four SF knots are collected in Table 1.

N-W region
We found the chemical properties of the SF knots are similar among them, with no O/H or N/O differences from clump to clump. However, the values of EW(H ) increase from the southeast region to the northwest one. Given that the EW(H ) traces the age of the burst, with higher values of EW(H ) meaning younger stellar populations, the star forming episodes of the galaxy seem to have a delay between them.
The brightest SF clump -central region in Table 1-has spectrophotometric data from SDSS (Ahumada et al. 2020). Berg et al. (2022) fit the ugriz SDSS photometry of the galaxy, complemented with FUV and NUV data from GALEX, and obtained a SFR value of log(SFR/(M yr −1 )) = −1.49 +0.21 −0.15 for this region, in agreement with our value from the extinction-corrected H emission.
The integrated spectra of the knots allow to better detect the intermediate components of A detailed visual inspection of the spectrum for the brightest complex of the galaxy (Figure 9) reveales the presence of numerous high excitation lines, such as [Fe ] or He , as well as the presence of the Wolf-Rayet (WR) red bump. A more detailed analysis of this feature is presented in the next section.

The ionising sources
All the emission lines across the galaxy suggest they are generated by strong star formation. As the EW(H ) values suggest, the age of the burst is no older than = 10 Myr. We now turn our attention to the nature of the ionising sources.
The integrated spectrum of the brightest SF complex of the galaxy presents a broad emission feature at 5808 Å which is typically associated with the presence of WR stars. WR stars are a late evolutionary stage of massive stars in which the strong stellar winds expel the outer layers of the star (e.g., Crowther 2007), exposing the innermost regions and reaching effective temperatures as high as 100,000 K (Smith et al. 2002). For single starburst episodes, the WR features appear between 3 to 5 Myr of the initial burst (Leitherer et al. 1999;Smith et al. 2002). In the optical part of the spectrum the two most prominent WR features are broad bumps which originate from the emission of He at 4686 Å-the so called blue bumpand the emission of C at 5808 Å-known as the red bump- (Schaerer et al. 1999;Smith et al. 2016). These broad emissions are generated by two different types of WR stars: nitrogen-rich WR (WN) are responsible for the blue bump whereas the red bump is generated by carbon-rich WR (WC) stars (Schaerer et al. 1999). In Figure 9, we zoom in the spectrum of the brightest clump around the WR red bump region (black solid line). The different emission lines are annotated with dotted lines. In order to determine the presence of WR in MaNGA galaxies, Liang et al. (2020) (2) where bump is the continuum-subtracted flux of the WR bump and rms is the noise of the spectrum around the baseline calculated in two windows at both sides of the red bump: 5660-5710Å and 5900-5950Å. Using the same formula in the collapsed brightest knot, the 5808 Å bump is detected with a significance of bump = 3.6. Based on the gas-phase metallicity derive, we assume the metallicity of the stars which ionised the gas should be 1/10 (Eldridge et al. 2017). At this metallicity, the presence of WR stars could be explained with models that incorporate binary populations. We use the models in order to reproduce this feature. We fit the region displayed in Figure 9 assuming the extinction follows the LMC (Gordon et al. 2003) and recover the best fit corresponds to an age of t = 5.1 ± 1.2 Myr. Senchyna et al. (2022) used deep UV spectroscopic data for the same region of the galaxy to derive an age of log(t/yr) = 6.7 +0.5 −0.2 , in agreement with our age estimation. The spectrum of this population is displayed in Figure 9 with a red solid line, with a vertical offset for better visualisation.   The best fit to the models is shown in orange, offset for visualisation purposes. High excitation lines, such as [Fe ] or He , are also annotated with dotted lines.
WR stars, as well as O-type stars, deposit a large amount of high energy photons into the surrounding H region. In an optically thick nebula, all the ionising photons emitted by these stars will be absorbed. However, in an optically thin nebula, a fraction of these photons is able to reach the surrounding gas, ionising the interstellar medium (ISM). The ratio between the expected and the observed H luminosity give us information about the fraction of photons escaping the ionising region. This escape fraction not only affects the available energy budget but also can probe the stars in the H region are the responsible of the ionisation of the ISM. The models provide the expected H luminosity in the case that all the ionising photons are absorbed by the nebula (Stanway & Eldridge 2018). For our population, with an age of 5 Myr and 1/10 , the expected luminosity is log( (H exp )) = 40.51 ± 0.01. The observed H luminosity of the brightest complex is log( (H obs )) = 40.23 ± 0.02. This implies that around 50% of the total ionising photons are absorbed, so that the total escape fraction of the central clump is esc = 0.5. Following a similar methodology, Della Bruna et al. (2021) inferred 0.2 < esc < 0.8 for typical H regions in NGC 7793, in agreement with our estimation for CGCG 007-025.

Estimation of the number of WR stars
The detection of the WR red bump allows us to estimate the total number of WC stars present in this complex. López-Sánchez & Esteban (2010) computed an expression of the luminosity of one WC star depending on the oxygen abundance measured. Considering our region has 12+log(O/H) = 7.9, the luminosity of the red bump generated by one WC star is WC = 1.56 × 10 36 erg s −1 . Since the measured bump has an extinction-corrected luminosity of red bump = (3.61 ± 0.56) × 10 37 erg s −1 , the number of WC stars in the region is (WC) = 23 ± 4. To obtain the total number of WR stars (WC+WN) in the region we need first to estimate the number Normalise flux blue bump red bump Figure 10. Comparison between the SDSS (orange) and MUSE (black) spectra using the same aperture. The red line represents the best fit model shown in Fig. 9. Both WR bumps are highlighted in grey. The presence of the red bump is clear in the MUSE spectrum, but its limited wavelength range does not cover the position of the blue bump. The SDSS spectrum covers both WR bumps, but its lower S/N does not allow a reliable detection (although both highlighted regions are systematically above the continuum).
of WN stars. As mentioned above, WN stars generate a broad bump emission at 4686Å, a region of the optical spectrum not accessible with the MUSE datacube but available in the SDSS spectrum. Despite covering a larger optical wavelength range, the WR bumps have not been reported in the SDSS spectrum. In Figure 10, we compare the SDSS spectrum of the galaxy (orange solid) with the MUSE spectrum (black solid) of the central knot, using the same aperture as in SDSS (i.e., 3 ). Both spectra have been smoothed with a median filter to an effective resolution of 3 Å, with the original spectra for each case in light grey. In red, the best fit model to the MUSE spectrum is displayed, as in Fig. 9. The location of the two WR bumps at 4680 Å and 5800 Åis highlighted in grey. The detection of the red bump in the MUSE spectrum is clear, but much less so in the SDSS spectrum. However, the flux in the two bump regions always is above the continuum level in the SDSS spectrum. We measure a low detection significance of blue bump = 1.3 and red bump = 1.2. Knowing from the MUSE data that the red bump is real, we model both bumps in the SDSS spectrum with Gaussian components and infer a flux ratio blue bump / red bump = 2.2 ± 0.7.
López-Sánchez & Esteban (2010) also derive a metallicitydependent relation for the luminosity of one WN, which corresponds to WN = 9.85 × 10 35 erg s −1 at the metallicity of this region. Using our measured value for the red bump in the MUSE spectrum, in addition with the ratio between the bumps derived from the SDSS spectrum, we estimate the number of WN stars as (WN) = 80 ± 25. The ratio of WC to WN stars then corresponds to WC/WN = 0.29 ± 0.09 and the estimated total number of WR in this region of the galaxy is N(WR) = 103 ± 26. With this value, we can perform a tentative estimation of the ratio of WR stars over O-type stars, (WR)/ (O). The number of O stars is related to the luminosity of the H line. For reference, a single O7V star produces a H luminosity of (H ) = 4.76 × 10 36 erg s −1 (Vacca & Conti 1992). O7V stars represent a fraction of the total population of O stars, 0 , which depends on the shape of the IMF, upper mass limit and age of the burst. Schaerer & Vacca (1998) calculate the time evolution of 0 ( ) for a Salpeter IMF with an upper mass limit of 120 , taking into account the evolution of the O star population as well as the emergence of WR stars. For a population with an age of 5 Myr, this value is 0 ≈ 0.25. The total H luminosity of the region is (H ) = (6.04 ± 0.06) × 10 39 erg s −1 , which yields to a total of (O) = 5080 ± 50. The ratio (WR)/ (O) derived for this region is (WR)/ (O) = (20 ± 6) × 10 −3 . As reference, Senchyna et al. (2017) used the SDSS spectrum of the central clump to analyse the possible WR population. Although they do not detect any broad feature in any of the two WR bump locations, they give a tentative value of the ratio (WR)/ (O) to be lower than 0.02, compatible with our findings. However, the BPASS models predict a total of 63 ± 1 WR stars as well as a ratio of N(WR)/N(O) = 0.012. Our estimations for both, the total number of WR stars and the WR stars over O-type stars ratio, are higher than the model predictions, mainly due to the large uncertainty in the estimation of WN stars.

The gas accretion scenario
From Figure 5 it is clear that the metallicity of the gas in CGCG 007-025 is not homogeneous across the galaxy. While the metallicity of the different star forming clumps is roughly the same (12 + log(O/H) ∼ 7.9), they stand out of the mean metallicity of the ISM (∼ 8.3). Remarkably, the SF clumps are at least 400 pc apart from the photometric centre of the galaxy.
In closed-box chemical evolution models (Tinsley 1980) the available gas reservoir in the galaxy is contaminated by previous SF episodes, so that each new burst of star formation has higher metallicity than the previous one. Moreover, the gas will fall to the centre of the potential well, generating SF regions close to the centre of the galaxy. Accordingly, the metallicity will be higher where the stellar mass density is also higher, creating metallicity gradients that typically decrease inside out (Vilchez et al. 1988;van der Kruit & Freeman 2011;Moran et al. 2012, among others). Thus, off-centred metallicity drops associated with star-forming regions cannot be explain within this scenario. Metallicity inhomogeneities can be generated via three different processes, that can operate simultaneosly: the mixing of metals within the ISM (Werk et al. 2011); strong outflows of metals (Chisholm et al. 2015); and inflows of metal-poor gas in tidally interacting systems (Pearson et al. 2019) or via clumpy gas accretion from the cosmic web (e.g., Sánchez Almeida et al. 2015;Ceverino et al. 2016). Below we give a general view of the three scenarios, and expose how our findings favour a metal-poor gas accretion event.

Mixing of metals
Metal production in galaxies traces star formation, and is highly concentrated toward the centres of galactic discs. In the absence of large-scale metal mixing, the density of metals should be directly proportional to the density of stars. However, several studies have found that metallicity gradients are flatter than density gradients (e.g., Queyrel et al. 2012). Petit et al. (2015) simulated an isolated disc galaxy and studied the effects of shear and turbulence in its metal distribution. They found the metal distribution is highly affected by both, with turbulence being the major agent in small scale structures and shear transferring momentum to turbulence from large scale to small scale dissipation. In the absence of perturbations, such as bars or spiral arms, large scale inhomogeneities are dissipated slower than small scale ones mainly by the turbulence of the gas. Petit et al. (2015) also shown that, in a axisymmetric system, the mixing timescales are typically a fraction of the orbital period, i.e. a few hundred of Myr. Galactic fountains also redistribute the gas mass on large scale and so smear out the metallicity inhomogeneities (e.g., Marasco et al. 2012;Sánchez Almeida et al. 2014a) Reconstruction of the star formation history for the underlying host galaxy (del Valle-Espinosa et al., in prep) indicates that star formation ceased over 1 Gyr ago. In the absence of any other mechanisms, metal mixing should produce a more homogeneous and flat oxygen abundance distribution than the one observed in the left panel of Figure 5.

Outflow of metals
The evolution of the most massive stars generates energetic winds and supernova explosions. These winds and blasts are the two main sources responsible of gas expulsion (Silich & Tenorio-Tagle 2017. In low-mass galaxies, the metal-enriched gas existing in the ISM can be easily swept out and returned to the IGM due to their shallow potential wells. The optical tracers of the stellar winds and supernova explosions are broad emission lines. As mentioned in Section 2.3, and displayed in Figure 4, the H emission in the brightest SF knot is described with a three-gaussian model. Our intermediate component around this region has an extension of 0.7×0.4 kpc 2 , with a median internal velocity dispersion of ∼200 km s −1 . Marasco et al. (2022) studied in more detail the broad components of a sample of star forming dwarf galaxies which includes CGCG 007-025. Their analysis presents two different approaches for the treatment of the broad components. Their Scenario 2 is more similar to our fitting procedure, and they found a wind component with an radial extension of 0.5 kpc and a velocity dispersion of 167 km s −1 . Their analysis, in alignment with other studies (Lelli et al. 2014;McQuinn et al. 2019) concludes that the baryonic feedback stimulates a gentle gas cycle rather than producing a large-scale blow out. Ejection of metals from the centre of the galaxy will produce -anomalous-enhanced O/H but also N/O ratios in the outer parts of the galaxy (Belfiore et al. 2015;Luo et al. 2021). With the outflow characterisation by Marasco et al. (2022), and since the N/O distribution remains flat across the galaxy, we can rule out the outflow of metals as the main explanation for the metallicity inhomogeneities mapped. . From top to bottom: oxygen abundance, nitrogen-to-oxygen and ionising parameter versus SFR surface density for every spectral bin. Markers are colour-coded by the distance to the centre of the brightest clump. While the oxygen abundance (ionising parameter) anti-correlates (correlates) with the SFR surface density, the nitrogen-to-oxygen shows no correlation with the Σ SFR .

The accretion of metal-poor gas
The simplest explanation for localised metallicity drops in a single system is the accretion of external metal-poor gas, induced by a recent merger, interaction or cosmological gas accretion. The same kind of metallicity deficit associated with bright star-forming regions has also been observed in several starbursting dwarf galaxies (Sánchez Almeida et al. , 2014b. This behaviour is fairly common at high redshift (Straughn et al. 2006;Windhorst et al. 2006;Elmegreen et al. 2007;Elmegreen & Elmegreen 2010), and in the local universe is usually associated to low mass ( ★ < 10 10.5 , e.g. Sánchez-Menguiano et al. 2019) and metal-poor galaxies (Papaderos et al. 2008;Morales-Luis et al. 2011;Filho et al. 2013).
Sánchez Almeida et al. (2018) studied the spaxel-to-spaxel anticorrelation between the gas-phase metallixity and the SFR in a sample of 14 local star-forming dwarf galaxies, expected if variable external metal-poor gas fuels the star-formation process. Similarly, in our Figure 11 we show the oxygen abundance against the SFR surface density for all the analysed bins (upper panel). The quantities are clearly anti-correlated, but it is necessary to discard variations in the physical conditions of the gas to confirm this anticorrelation is triggered by metal-poor gas accretion. We also display the log(N/O) (log(U)) versus the SFR surface density in the middle (lower) panel of Figure 11. Changes in the nitrogen-to-oxygen ratio are linked with changes in the chemical composition of the gas, whereas the ionisation parameter changes with the age and mass of the stellar population. The N/O remains constant across the galaxy, while the log(U) correlates with the Σ SFR . These three relations are naturally explained if a significant amount of low metallicity gas was accreted, reducing the gas-phase metallicity and enhancing star formation, but not altering the relative metal content of the ISM, including N/O. Moreover, the underlying population of CGCG 007-025 is much older than the current SF burst (1 Gyr versus 5 Myr, respectively), leaning the scale to the gas accretion scenario rather than gas consumption as the main responsible of the recent SF episode.

CONCLUSIONS
Local starbursting galaxies present a unique environment to study star formation processes akin to those at higher redshift with much greater sensitivity and spatial resolution. In this paper we present a study of the chemodynamical properties of the ionised gas in CGCG 007-025. CGCG 007-025 is a dwarf spheroidal galaxy undergoing an off-centre starburst event-most likely as a result of an interaction with a nearby dwarf companion. Our main results, based on the analysis of the ionised gas present in the galaxy, are summarise below: • We derive the velocity maps for the different kinematic components of the H line. The velocity field of the narrow component is highly perturbed, with no indication of rotation. The intermediate component presents a velocity profile compatible with a bipolar structure in expansion. The broad component, only detected in H and in the brightest SF knot, displays a velocity distribution which resembles to a dusty expanding bubble, and could be originated in denser gas ( ∼ 10 6 cm −3 ).
• Using HII-CHI-(Pérez-Montero 2014), we obtain the chemical properties of the gas. We find a homogeneous distribution of N/O, whereas the distribution of oxygen abundance decreases at the position of the SF clumps. The ionisation parameter distribution presents the opposite behaviour.
• We compute line ratios between species with different ionisation potentials. The • Thanks to the depth of the MUSE datacube, we are able to detect the WR red-bump in the brightest SF complex of the galaxy.
The detection of the bump allows to make a better estimation of the age of this burst, 5.1 ± 1.2 Myr, as well as the number of WR stars present in it, 103 ± 26.
• The overall relation between the distribution of metals and the Σ SFR suggest the accretion of metal poor gas to be responsible for the starbursts we are observing now in this galaxy. This paper is devoted to the physical properties and origin of the ionised gas in CGCG 007-025. However, this galaxy presents an underlying old population as well as a NSC. The stellar properties of this galaxy will be discussed in a companion paper (del Valle-Espinosa et al., in prep.).

APPENDIX A: LIST OF EMISSION LINES
This paper has been typeset from a T E X/L A T E X file prepared by the author.