© The Authors 2024

1. Introduction

The ability to resolve stellar populations in nearby galaxies allows us to study the physical processes that govern their evolution. The dwarf galaxies of the Local Group (LG), in particular, provide an excellent laboratory thanks to their diverse characteristics, including size, gas content, star formation history (SFH), chemical abundances, and kinematics (e.g. Tolstoy et al. 2009; McConnachie 2012; Kirby et al. 2013; Gallart et al. 2015; Simon 2019; Putman et al. 2021; Battaglia & Nipoti 2022). Since these characteristics can be influenced by the proximity of a larger galaxy, such as the Milky Way and Andromeda, isolated LG dwarf galaxies provide valuable insights into the intrinsic mechanisms that regulate their evolution. In the gas-rich ones, the interplay between gravitational instabilities and stellar feedback processes can be reconstructed by combining the kinematics and spatial distribution of stars of different ages with the SFH and gas content of the galaxy (e.g. Leaman et al. 2017; Collins & Read 2022). Gathering this kind of information is usually an observational challenge, but the advent of integral field spectroscopy (IFS) has been transformative in this respect, due to the ability to inspect different galactic tracers with high observational efficiency (Roth et al. 2018, but also Sánchez et al. 2012; Cortese et al. 2014; Bundy et al. 2015; McLeod et al. 2020; Zoutendijk et al. 2020; Júlio et al. 2023; Vaz et al. 2023).

Among LG dwarf galaxies, IC 1613 is an ideal IFS target. On the one hand, it is a typical low-mass star-forming dwarf irregular galaxy (Skillman et al. 2014). On the other hand, stellar feedback processes have generated voids and bubbles in the interstellar medium whose impact on stellar and gas kinematics needs to be explored (Read et al. 2016). First discovered by Wolf (1906), IC 1613 was recognised afterwards as an extragalactic object by Baade (1935, as reviewed by Sandage 1971), who first determined its distance using Cepheid variable stars. Since then, several other measurements of its distance have been made using different indicators, including RR-Lyrae variable stars and the tip of the red giant branch (RGB). A compilation of literature distance values (including their own) was reported by Bernard et al. (2010), who obtained a statistical average value of (m − M)₀ = 24.400 ± 0.014, or 759 ± 5 kpc.

IC 1613 therefore lies well beyond the virial radius of the Milky Way (MW) and that of M 31 (being at ∼520 kpc distance from it, McConnachie 2012). The systemic proper motion of IC 1613 obtained from McConnachie et al. (2021) with Gaia DR2 data does not rule out the possibility that the galaxy was ever within 300 kpc of M 31, for M 31 virial masses of ≳1.3 × 10¹² M_⊙. However, the most recent estimate of the systemic proper motion of IC 1613 by Bennet et al. (2023), combining data from Gaia eDR3 and the Hubble Space Telescope (HST), makes such an association unlikely when assuming an M 31 virial mass of 2 × 10¹² M_⊙. Interestingly, Buck et al. (2019) assign to IC 1613 a high probability of being a backsplash galaxy (i.e. a currently isolated system that may have once passed close to a large host, but without becoming bound) of M 31.

Its low luminosity (an absolute magnitude of M_V = −15.2 mag) corresponds to a total stellar mass of M_* ∼ 10⁸ M_⊙, similar to that of other LG dIrrs like NGC 6822 and WLM (McConnachie 2012). The SFH and stellar content of IC 1613 have been studied in detail over the years (see e.g. Cole et al. 1999; Skillman et al. 2003, 2014; Bernard et al. 2007). Based on deep HST Advanced Camera for Surveys (ACS) imaging, and after a review of archival HST data, Skillman et al. (2014) concluded that the SFH of IC 1613 has been close to constant on average throughout its life, with no evidence of an early dominant episode of star formation. We note that the HST/ACS field analysed by Skillman et al. (2014), roughly located at the half-light radius of the galaxy, was considered as representative of the global SFH of IC 1613. The comparison with archival HST observations located at different radii supported this assumption. The average metallicity was linearly increasing through time, with values ranging from [Fe/H] ∼ −2.0 dex to −0.8 dex. Stars were formed at an average rate of ψ(t) = 0.081 ± 0.001 M_⊙ yr⁻¹.

The structural properties of IC 1613 have been determined by studying the spatial distribution of stellar tracers in different evolutionary phases (see e.g. Albert et al. 2000; Borissova et al. 2000; Bernard et al. 2007; Garcia et al. 2009; Sibbons et al. 2015; McQuinn et al. 2017; Pucha et al. 2019; Higgs et al. 2021). In general, it was found that the young stars are more centrally concentrated than the intermediate-age and old stars, a feature that is common in dwarfs. However, Pucha et al. (2019) using deep and wide Subaru/Hyper-SuprimeCam observations of IC 1613, showed that its young main sequence (MS) stars, along with the RGB and ancient horizontal branch (HB) stars, all extend to the outskirts of the galaxy, up to ∼24 arcmin (i.e. ∼4 effective radii R_e). In particular, the young stars are found well beyond the currently active star-forming regions (within ∼1.5 × R_e), although with a much lower density compared to the intermediate and old age components. They also showed a steeper radial density profile in the inner regions than in the outer ones, in contrast to the RGB and HB stars. This seems to imply a different formation channel between the younger (e.g. from gas pushed outwards through stellar feedback) and older stars (probably via gas accretion) in the galaxy’s outskirts. In this work, we assume the structural parameters from Higgs et al. (2021) who conducted a homogeneous analysis of the isolated LG dwarf galaxies.

As other dIrrs, IC 1613 contains an extended component of neutral hydrogen (HI) gas, with a clumpy distribution on the inside, but showing regular contours at larger radii (see van den Bergh 2000, and references therein). In particular, the HI distribution is rich of shells and voids around the currently active star forming regions, where the ionised gas is also distributed (e.g. Lozinskaya et al. 2003; Silich et al. 2006; Moiseev & Lozinskaya 2012; Pokhrel et al. 2020; Yarovova et al. 2024). The HI kinematics shows a linearly increasing rotation curve with a v_max ∼ 20 km s⁻¹ (Lake & Skillman 1989; Oh et al. 2015; Read et al. 2016). The rotation value, however, may represent a lower limit due to the fact that this galaxy is probably seen close to face-on. This translates into a significant uncertainty in the determination of the dynamical mass being between M_dyn ∼ 5 × 10⁸ M_⊙ and ∼8 × 10⁹ M_⊙, depending on the assumed inclination angle (Read et al. 2016).

The first spectroscopic measurements of stars in IC 1613 were obtained from some of the brightest ones, that is, from nine early B-type young supergiants (Bresolin et al. 2007) and three evolved M-type supergiants (Tautvaišienė et al. 2007). Both studies obtained compatible metallicity values, 12 + log(O/H) = 7.90 ± 0.08 dex and [Fe/H] = −0.67 ± 0.09 dex, respectively¹, later confirmed by Berger et al. (2018) who studied 21 young BA-type supergiant stars. Interestingly, these authors found a bimodal metallicity distribution that appears to be correlated with their spatial location (i.e. lower metallicity stars are found in high density HI regions). A comparable bi-modal distribution was also observed by Chun et al. (2022) for a sample of 14 red supergiants. Furthermore, there is evidence from the youngest stellar population and evolved red supergiants of IC 1613 that the present-day [α/Fe] ratio could be sub-solar (i.e. ∼ − 0.1 dex, Tautvaišienė et al. 2007; Garcia et al. 2014). The chemical abundance of the interstellar medium has been obtained spectroscopically from the central HII regions and is in good agreement with results from the supergiant stars (Bresolin et al. 2007, but see also Lee et al. 2003; Tautvaišienė et al. 2007).

Kirby et al. (2013) were the first to present a statistically significant spectroscopic analysis of a sample of 125 RGB stars in IC 1613, which led to an average metallicity of [Fe/H] = −1.19 ± 0.01 dex, in agreement with the fact that older stars are in general more metal-poor than the younger ones, and with expectations from the age-metallicity relation found by Skillman et al. (2014). On the same sample, Kirby et al. (2014) conducted a kinematic analysis which led to a systemic velocity of −231.6 ± 1.2 km s⁻¹, in agreement with values from the HI component (Lake & Skillman 1989; Oh et al. 2015), and a velocity dispersion of $10.8_{-0.9}^{+1.0}$ km s⁻¹. They did not find signs of rotation for the stellar component; this was also later confirmed by Wheeler et al. (2017) who re-analysed the Kirby et al. (2013) dataset. This is not surprising considering that the spatial distribution of this dataset is roughly perpendicular to the HI major kinematic axis (Lake & Skillman 1989; Oh et al. 2015), assuming that the stellar rotation follows that of neutral gas. The dynamical mass at the half-light radius reported by Kirby et al. (2014) was M_1/2 = 1.1 ± 0.2 × 10⁸ M_⊙.

In this paper we present a study of the stellar kinematic and chemical properties of IC 1613 from spectroscopic data taken with the MUSE integral field spectrograph on the Very Large Telescope (VLT). We took advantage of the unique capabilities of this instrument, that combines high spatial resolution with a 1 arcmin² field of view and a wide wavelength range at medium spectral resolution (Bacon et al. 2014). We then performed a kinematic analysis as a function of stellar age and obtained metallicities for the largest spectroscopic stellar sample obtained to date for this galaxy. IC 1613 is located at high Galactic latitude (b ∼ 61 deg), so its low values for both foreground (Schlafly & Finkbeiner 2011) and internal reddening (Georgiev et al. 1999) make it an ideal laboratory to compare the observed properties of its diverse stellar content.

The article is structured as follows. In Sect. 2 we give details on the data acquisition and reduction processes we conducted. Section 3 is dedicated to spectral classification and velocity determination. In Sect. 4 we report details on the determination of likely member stars and their subsequent kinematic analysis, while in Sect. 5 we show results of the chemical analysis of the RGB stars in our sample. Section 6 is dedicated to the discussion of the age-kinematic trends of the different tracers inspected in this work. In this section, we also discuss the dynamical mass estimation and the role played by the galaxy’s inclination angle. Finally, Sect. 7 is dedicated to the summary and perspectives of future work, while in the appendices we report details on the sanity checks conducted during our analysis. The parameters adopted for IC 1613 throughout the text are summarised in Table 1.

2. MUSE data of IC 1613

2.1. Observations and data reduction

The data were acquired with VLT/MUSE in service mode². MUSE is an integral field spectrograph with a spatial sampling of 0.2″ and a spectral resolution varying between R = 1500 − 3000 along its wavelength coverage. We used the instrument in nominal wide field mode, with a field of view of ∼1′×1′ and a wavelength coverage between 4800 − 9300 Å. We observed three fields (hereafter F1, F2, F3, moving outwards from the centre; see Fig. 1 for their location) approximately along the projected semi-major axis of IC 1613, on the west side of the galaxy. We purposely avoided regions of ongoing star formation.

Fig. 1. Finding chart showing the location of the three MUSE fields, marked as white boxes, overlaid over a wide-field image of IC 1613 (Credit: ESO – VST/Omegacam Local Group Survey). The white ellipse marks the half-light radius. The smoothed HI column density map from the Little-THINGS survey (Hunter et al. 2012) is marked with logarithmically spaced isodensity contours (silver lines) starting at 3σ (white-smoke lines). North is up, east to the left.

Each field was observed with a total exposure time of 11 080 s on source, split in 4 × 2770 s exposures, obtained by rotating the position angle of the spectrograph by 90° with respect to the previous exposure. No separate sky field was acquired since the stellar component of IC 1613 is resolved in the magnitude range of interest and the sky contribution is taken care of during the de-blending and spectral extraction phase of the data analysis (see Sect. 2.2). The data were acquired in dark time, clear sky, with a request for a seeing in V-band at zenith ≤0.9″, fulfilled for all but one exposure of F1 (for which it deviated within 10% from the request). The observing log is reported in Table C.1.

All data cubes were reduced using version 1.6 of the official MUSE pipeline (Weilbacher et al. 2012). The basic steps of the reduction cascade (bias subtraction, slice tracing, wavelength calibration, and basic processing) were performed using the default settings³. This resulted in 24 pixel tables (one for each unit spectrograph) for each individual exposure. When combining the pixel tables for each exposure, flux calibration was also applied using standard stars observed on the same nights. Furthermore, the subtraction of sky emission lines was performed by determining their intensity directly from the scientific data, which contained sufficiently large patches of (almost) blank sky for this purpose. The sky emission lines were also used to quantify the quality of the wavelength calibration, which had an average accuracy of 0.03 Å, or 1.5 km s⁻¹. The final step of the data reduction was the combination of all individual exposures obtained for a pointing. The end products of the data reduction were three data cubes, one for each field. Each cube has a dimension of approximately 300 × 300 × 3680 pixels. In Fig. 2, we show the broad- and narrow-band images obtained from the reduced data cubes.

Fig. 2. Broad- and narrow-band images obtained from the MUSE data cubes. Top row: Combined images in the VRI-bands highlighting stellar sources and background galaxies. Bottom row: Combined images highlighting the emission-line ionised gas with Hα coloured in red, [SII] (6713 Å) in blue, and [OIII] (5007 Å) in green. In all panels, North is up and east is to the left.

2.2. Source extraction and photometric cleaning

The PampelMUSE software described in Kamann et al. (2013) was used to extract the stellar spectra from the reduced data cubes. To run PampelMUSE, an input catalogue of source positions and an estimate of their initial magnitudes is required. Since the MUSE data are only moderately crowded, the input catalogues were directly created from the data cubes. We produced synthetic broad-band V-, R- and I-band as well as emission line colour-composite images (see Fig. 2), using the procedures of Roth et al. (2018), and extracted photometric catalogues using the DAOPHOT code (Stetson 1987). The photometric catalogues were obtained as in Gallart et al. (1996). They proved to be a valuable resource for our analysis, as they contain the spatial and photometric information of all the extracted sources. The details of their astrometric and photometric calibration are reported in Appendix A.

The raw source catalogues (i.e. without astrometric solution or photometric calibration applied) were passed as input to PampelMUSE to first identify the positions of stars in the data cubes as a function of wavelength. The code then aimed to retrieve the point spread function (PSF) using an analytical Moffat profile. To obtain the final spectra, a simultaneous PSF fit was performed to all sources in each layer of the data cube; unresolved components were treated by including a background in the fit that was recalculated at fixed spatial offsets. As IC 1613 contains gaseous emission that varies on rather small spatial scales, well visible in Fig. 2, the background component was recomputed every 40 spaxels.

PampelMUSE extracted 2293 initial spectra, examples of which are shown in Figs. C.1, C.2 and C.3. For each extracted spectrum, PampelMUSE provided the associated flux uncertainty per pixel and an S/N value calculated around its central wavelength (hereafter S/N_C). However, as our sources cover a wide spectral range (see Sect. 3.1.1), we calculated two further S/N indicators: one around 5500 Å (S/N₅₅₀), and the other on the continuum around the CaT lines (S/N_CaT). This gave us a more accurate picture of the quality of our spectra depending on the spectral type of a given star.

The photometric catalogues were used to perform a first cleaning of the spectroscopic sample, excluding clearly non-stellar sources. We selected targets based on specific photometric parameters, including sharpness between −0.5 and 0.5, and a goodness-of-fit (CHI) of less than 1. We also excluded targets with magnitudes in only one band, which were generally low S/N targets identified only in the I band. As a result, the total number of extracted sources decreased from 2293 to 2053.

The quality of our data can be seen in Fig. 3. In particular, we obtained respectively in F1, F2, F3 about 65, 40, 10 sources with S/N₅₅₀ ≥ 20, and about 200, 120, 45 sources with S/N_CaT ≥ 10. We note that for this galaxy dust reddening and extinction along the line of sight (l.o.s.) are almost negligible, with average values of A_V = 0.067 and A_I = 0.038 (Schlafly & Finkbeiner 2011). Nevertheless, we have taken this into account.

Fig. 3. Extinction-reddening corrected colour-magnitude diagrams from MUSE images of the three fields, moving from the centre towards the outer regions from left to right; the filled circles represent the sources for which spectra have been extracted and they are colour-coded by their S/N at the central wavelength. We note that the galaxy’s distance modulus is (m − M)0 = 24.4 (Bernard et al. 2010).

3. Spectral types and velocity determination

We performed an analysis to determine the spectral type and l.o.s. velocity of our sources. For this purpose, we used, respectively, the spectral fitting codes ULySS (Koleva et al. 2009) and SPEXXY (Husser 2012). The reason for this dual approach is motivated as follows.

Crowded field IFS can be considered an accomplished technique for globular clusters (GCs; Giesers et al. 2018, 2019; Latour et al. 2019; Kamann et al. 2020; Saracino et al. 2022). However, for galactic systems at larger distances, involving young stellar populations and ongoing star formation, the method is still at a stage of exploration. Unlike GCs that are dominated by old stars, complications related to hot stars, emission line stars, carbon stars, but also to unresolved star clusters and the presence of strong emission lines from H II regions, make it difficult at present to blindly use automated fitting tools such as SPEXXY for nearby galaxies. The scheme of visually assisted spectral fitting described below was first developed in a pilot study of the disk galaxy NGC 300 at a distance of 1.88 Mpc (Roth et al. 2018, hereafter R18).

3.1. Spectral classification with ULySS

3.1.1. General method and classification results

Following R18, we used ULySS (Koleva et al. 2009)⁴ to fit the observed spectra to a linear combination of ten templates with different weights. Templates were extracted from the MIUSCAT empirical library⁵, rather than from a model grid. We note that we supplemented the MIUSCAT library with 19 spectra of carbon stars from the X-shooter library (Chen et al. 2014)⁶. We also added a set of mock spectra of Be stars composed ad hoc from library B-star spectra, combining two-component Gaussians at Hα and Hβ wavelengths with different equivalent widths. While these modifications did not aim to derive significant stellar parameters, they proved useful in identifying candidate stars that are not included in the MIUSCAT library for their later visual confirmation and velocity determination (see Sect. 3.2.2).

The linear combination fitting mode of ULySS has been employed in R18 with the experience that unresolved blends from different spectral types are flagged to prompt an inspection. However, we found that the mild crowding in IC 1613 did not present us with such cases, unlike the more distant galaxy NGC 300. Occasionally though, the fit returned implausible mixes of spectral types (e.g. mixing young early-types with evolved late-types) for low S/N stars. The apparent magnitude obtained by shifting the absolute magnitude in R-band of the stars in the MIUSCAT library to the distance of IC 1613 (distance modulus (m − M)₀ = 24.40, from Bernard et al. 2010), provided a way to weed out these clearly erroneous fits when they would lead to apparent magnitudes well below the detection limit of the MUSE data. We also found that good-quality fits to cool MS stars allowed us to confidently identify Milky Way foreground stars that would otherwise have been difficult to reject from photometry alone.

The ULySS spectral analysis was carried out on all those sources having a S/N_C >  2 and a clean photometry (see Sect. 2.2). Main outputs of the code were the best fitting templates, together with the associated spectral parameters (T_eff, log(g) and [Fe/H]) and l.o.s. velocity. We verified by eye the ULySS outcomes, storing the information of the most probable template chosen among those with the largest weights. We also assigned several quality flags to each analysed object, evaluating the quality of the input spectrum (QSP), the quality of the spectral fitting (QFT), the plausibility of the output l.o.s. velocity (PVR) and the plausibility of the spectral type classification (PCL). Flags were reported as integer numbers ranging from one to four, with the lower value meaning a useless measure, and the higher value implying an excellent fit. We refer to R18 for further details.

In this work, we used ULySS exclusively for spectral classification purposes, while we left the l.o.s. velocity determination part to the SPEXXY code. This was mainly because SPEXXY performed better than ULySS in the velocity determination task, especially in the low S/N regime (i.e. ≲10), as we verified in Appendices B.1 and B.2.

In Table 2, we report results from the spectral classification obtained for sources marked with a PCL flag between three and four (which implies an accuracy on the effective temperatures better than 500 K). The visual inspection during the spectral classification also allowed us to identify other contaminants that were not removed during the previous photometric selection (Sect. 2.2): they were mainly background galaxies (comprising high-z emitters) and low S/N targets highly polluted by diffuse ionised gas emission lines. We also visually identified carbon and B emission line stars that were missed by ULySS due to the lack of adequate stellar templates, in order to study them separately (see Sect. 3.2.2). We note that we define as bona-fide B emission-line stars (noted as Be hereafter) those sources whose spectral distribution is compatible with a star having T_eff >  10 000 K, but which show Hα and (occasionally) Hβ in strong emission (e.g. Porter & Rivinius 2003; Rivinius et al. 2013). The other seven stars showed only Hα in weak emission, but due to their low S/N their classifications were too uncertain and they were therefore discarded.

Overall, we obtained a reliable spectral classification for 808 stars⁷. Looking at Table 2, the large majority of sources are K-type giants, with an almost even representation of the other types, except for the M and O stars being a minority. In particular, among the M stars, only one was found to be a giant, while the other four were classified as MW-foreground dwarfs. We also report the presence of two other MW-contaminants classified as K V and G V stars, all confirmed in the l.o.s. velocity determination step. We have therefore excluded these MW-foreground stars from the counts in Table 2 and the following analysis. In Fig. 4 is shown the position of the classified stars on the colour-magnitude diagram, colour-coded according to their T_eff, visually confirming the general goodness of the spectral classification.

Fig. 4. Colour-magnitude diagram for data selected with PCL-flag between 3 and 4, as described in Sect. 3.1.1. Filled circles are colour-coded according to their effective temperature Teff as derived with ULySS; black triangles indicate the identified C-stars, while the grey squares mark the Be stars.

3.1.2. O and Be stars

The number of O-type stars identified is the most uncertain, since the poor coverage of the MIUSCAT library at the highest effective temperatures (i.e. for T_eff >  30 000 K) is a limiting factor for the spectral classification of these stars. We have identified 13 stars with a T_eff just above 30 000 K. This is in general agreement with the expected number of O stars in IC 1613. Indeed, considering the galaxy’s current star formation rate (Skillman et al. 2014) and assuming a Kroupa (2001) initial mass function, we would expect on the order of ten O stars in the surveyed area, given their average lifetime of 10 Myr and minimum mass of 16 M_⊙.

On the other hand, O stars are typically not uniformly distributed in space. They tend to form and evolve in OB associations. Our observations cover some of the OB associations identified by Garcia et al. (2009). These associations are low density groups of young stars with significant internal extinction, making the identification of O types more challenging. In Fig. 5 we show the spatial distribution of the OB associations reported by Garcia et al. (2009), together with that of our identified OB stars with T_eff >  20 000 K and also the Be stars. We see that they generally tend to be found where the OB associations are. Their spatial distribution also generally follows that of the ionised shells shown in Fig. 2, to which they are probably physically associated (Borissova et al. 2004; Garcia et al. 2009). At this stage, this evidence confirms the goodness of our initial classification of many of them as hot stars.

Fig. 5. Spatial distribution of hot OB stars (blue filled circles) and Be stars (darkslate grey filled squares), compared to that of the OB-associations (grey filled hexagons) reported by Garcia et al. (2009). The black contours indicate the Hα emission within the MUSE pointings, marked by black squares.

We further note that the identified Be stars account for up to 18% of the total OB sample. This is comparable to the observed Be fraction in nearby metal-poor dwarf galaxies (between 15% and 30%, Schootemeijer et al. 2022; Gull et al. 2022; Vaz et al. 2023). Since Be stars are likely to be highly rotating massive stars, our results add to the evidence that they are common in low-metallicity environments.

We are currently preparing a follow-up study using a grid of model atmospheres obtained from the NLTE code FASTWIND (Puls et al. 2005), whose spectra will allow us to fit the more massive stars of spectral type A to O that are sparsely covered by the empirical MIUSCAT library. Nevertheless, the wavelength range of MUSE, which does not cover the gravity and temperature sensitive lines below 4800 Å, is a limitation for the analysis of massive stars. An accurate determination of these parameters, as well as Fe and α-elements abundances, will be possible with the future BlueMUSE instrument (Richard et al. 2019).

3.1.3. C/M fraction

Looking at the coolest stars in our sample, we have identified a single M-type giant and 14 carbon (C) stars. The latter are a type of asymptotic giant branch (AGB) stars, characterised by having more carbon than oxygen in their atmospheres. This happens during the thermally pulsating (TP) phase, when C-rich material from the star’s interior is dredged up into the previously oxygen-rich atmosphere (a process known as the third dredge-up). A TP-AGB star is initially a O-rich M-type giant. Since it is easier to transition from an O-rich to a C-rich star when its metallicity is low, the C/M fraction is usually an indirect [Fe/H] indicator of intermediate-age stars in a galaxy. For IC 1613, literature works report an average value of C/M ∼ 0.6 (see e.g. Albert et al. 2000; Chun et al. 2015; Sibbons et al. 2015; Ren et al. 2022, based on wide-area JHK-photometric surveys).

We can qualitatively determine the C/M in our spectroscopic sample. The only M-type giant was clearly an O-rich AGB star. We also examined the K-type stars above the tip of the RGB and selected a sample with T_eff within 500 K of an early M-type star. Cross-correlating with the literature photometric catalogues (Albert et al. 2000; Chun et al. 2015; Sibbons et al. 2015), we found ten common targets classified as M stars (after excluding the possible red supergiants identified by Ren et al. 2022). As for our C stars, three have already been identified in the literature (see again Albert et al. 2000; Chun et al. 2015; Sibbons et al. 2015), while five have a more uncertain spectral classification due to their faint magnitudes (I >  21), well below the tip of the RGB. It should be noted that TP-AGBs are long-period variables (with periods of a few tens to hundreds of days), which can exhibit brightness excursions of several magnitudes (see e.g. Menzies et al. 2015). Alternatively, faint C stars can be binary systems in which the primary star has received C-enriched material from the secondary star, which was previously an AGB (also known as extrinsic C stars). For IC 1613, the expected fraction of extrinsic C stars is rather low (∼10%, e.g. Hamren et al. 2016). Considering between 9 and 14 C stars, the C/M ratio then ranges between 0.9 and 1.4, which is larger than the average value reported in the literature. However, taking into account the Poisson error associated with the small number statistics of our sample, we would remain in agreement with the literature.

3.2. Velocity determination

3.2.1. Velocity determination with SPEXXY

The determination of l.o.s. velocities was carried out using the SPEXXY code (v. 2.5, Husser 2012)⁸. The routine performs a full spectral fitting to the observed spectra using interpolated spectral templates generated from the PHOENIX library of high-resolution synthetic spectra (Husser et al. 2013)⁹. The library covers a wide wavelength range, going from 500 Å to 5.5 μm, and stellar parameters: 2300 < T_eff (K) < 15 000, 0.0 < log(g) (dex) <  + 6.0, −4.0 < [Fe/H] (dex) < +1.0 and −0.2 < [α/M] (dex) < +1.2. However, the PHOENIX models do not treat radiative transfer in the atmospheres of massive stars affected by strong stellar winds, so the library is of little use in the recovery of stellar parameters for T_eff >  10 000 K. We refer to Appendix B.1 for a detailed explanation of the spectral fitting steps and to B.2 for a comparison of its performance with the ULySS results presented above.

We ran SPEXXY on all targets having S/N > 2 (in at least one of the S/N indicators) and clean photometry (see Sect. 2.2). We also excluded all those targets that in the spectral classification step (see Sect. 3.1.1) were identified as contaminated by ionised gas, background galaxies, and Be stars (which we analysed separately). These steps reduced our sample to 1745 input sources. We were thus able to recover their l.o.s. velocity, along with estimates of their spectral parameters and associated errors. We note that only for 2% of stars SPEXXY failed to perform a successful fitting, mostly due to their low S/N.

The goodness of the recovered parameters was quantified in a series of tests, whose details are given in Appendix B.1. Here we report the essentials. We verified by performing several tests on mock spectra that the l.o.s. velocity errors are well estimated down to S/N_CaT ∼ 3.5 for red giant stars (S/N₅₅₀ >  5 for MS stars), and the velocity values are recovered without significant offsets. On the other hand, the recovery of the spectral parameters was more limited, with T_eff in particular being generally well constrained for red giants, but underestimated for the hottest MS stars due to the grid limitation of our spectral models. We emphasise again that our aim was to use SPEXXY mainly for the determination of l.o.s. velocities.

Since our main goal is to analyse not only the global kinematic properties of IC 1613, but also how they change as a function of stellar types, we decided to keep the sample obtained during the spectral classification step. From this sample we further removed those targets with a velocity error δ_v >  25 km s⁻¹, which meant excluding 10% of sources with the least reliable velocity values (mostly low S/N MS stars). We also excluded the six targets previously classified as MW contaminants, confirming that they are foreground sources based on their l.o.s. velocities, which were consistent with 0 km s⁻¹. Performing a cross-correlation with the third data release of the Gaia catalogue (Gaia Collaboration 2023), we further confirmed their foreground nature based on their proper motion values, significantly different from 0 km s⁻¹ in this case. With these conditions applied, our sample reduced to 727 sources to which we subsequently added the 24 Be stars whose velocity determination we describe in the next section. The distribution of l.o.s. velocity measurements for our sample is shown in Fig. 6, with histograms divided by pointing and stellar type. We note that the median velocity error was of δ_v ∼ 6.5 km s⁻¹, while the median S/N_C resulted around 10.

Fig. 6. Histogram of the l.o.s. velocity measurements divided by pointing (left panel) and stellar type (right panel).

3.2.2. Velocity determination for carbon and Be stars

The l.o.s. velocities of the C stars identified by spectral classification (see Sect. 3.1.1) were calculated using SPEXXY as for the main sample. However, the code was not always able to match the complex spectral features of these stars. Therefore, we double-checked the SPEXXY velocity determination with ULySS using stellar templates from the X-shooter library of carbon stars (Gonneau et al. 2016). We found a general agreement within the errors, confirming the goodness of the SPEXXY velocities. The median velocity error of the 14 C stars analysed (6 in both F1 and F2, and 2 in F3) was δ_v ∼ 4 km s⁻¹ for a median S/N_CaT ∼ 15.

On the other hand, the velocity determination for the Be stars was done separately, since SPEXXY failed in this task because such stars are not represented by the PHOENIX model atmospheres. We thus performed a cross-correlation with custom-made template spectra using the function correlation of the python package specutils, part of the Astropy project (Astropy Collaboration 2022). Templates were generated as ad-hoc proxies using the spectra of two B-type stars (classified as B9III and B3III) from the MIUSCAT library (Vazdekis et al. 2012), to which we added two Gaussian profiles reproducing the strength and width of the Hα and Hβ emission lines as detected in our sample. Since the spectra of our Be-star candidates in general did not show any split line profiles (except for two objects), the two-component Gaussian approximation seemed to be a good enough approach for the purpose of l.o.s. velocity measurements. We cross-correlated the observed spectra with the templates and assigned them the average of the derived velocities. The errors were instead assigned by a Monte Carlo process. For each template we generated mock spectra at different S/N (calculated around 5500 Å) values [2.5, 5, 7.5, 10, 15, 25, 40, 60]. For each S/N, we cross-correlated 250 mock spectra with the corresponding noise-free template and considered the median absolute deviation of the velocity distribution as the velocity error representative of that bin. The error distribution as a function of S/N was fitted by an exponential profile for each template and we used the mean profile to assign velocity errors to the observed spectra according to their S/N₅₅₀. For the 24 Be stars analysed (14, 6 and 4 in fields F1, F2 and F3, respectively), we recovered a median δ_v ∼ 7 km s⁻¹ for a median S/N₅₅₀ ∼ 20.

We note that if we perform the kinematic analysis described in the next section on both the Be and MS star samples, we find no significant differences between them in terms of systemic velocity and velocity dispersion. Therefore, we included them in the main sample when we performed the kinematic analysis.

4. Internal kinematics

We investigate the basic internal kinematic properties of the stellar component of IC 1613, such as its systemic l.o.s. velocity, its velocity dispersion and the possible presence of velocity gradients, indicating rotation of the stellar component¹⁰.

To proceed with the kinematic analysis, we first need to identify the possible contaminants left in our sample. We had already removed many of them during the spectral classification process, with the rest expected to be faint foreground stars. We used then the method outlined in Taibi et al. (2020), which allows a Bayesian kinematic analysis while assigning membership probabilities P_Mj, assumed to be Gaussian, to the individual targets. This approach is based on the expectation maximisation technique presented in Walker et al. (2009).

The P_Mj values depend on the l.o.s. velocities of the individual targets, but also on the prior information based on their radial distances from the galaxy’s centre and on the expected l.o.s. velocity distribution of the possible contaminants. The spatial prior takes into account that the probability of membership is higher, the closer to the galaxy’s centre. In this case, we simply required that our targets should follow a monotonically decreasing radial density profile applying an isotonic regression model, with no assumptions about its functional form¹¹.

We carry out a Bayesian analysis to explore and compare different kinematic models, considering one where the internal kinematics of IC 1613’s stellar component can be purely described by random motions, and another model which also contains a rotational term. For the rotational term we used the form: v_rot(R_j) cos(θ − θ_j), where R_j is the angular distance from the galaxy’s centre (i.e. $R_j=\sqrt{{\xi }_j^2+{\eta }_j^2}$, where (ξ_j, η_j) are the tangent plane’s coordinates), θ_j is the position angle (PA) of the j-target, and θ is the PA of the kinematic major axis (both PA measured from north to east). The rotational velocity term is assumed to linearly increase with radius, that is v_rot(R_j) = k R_j.

We note that, unless otherwise stated, for the stellar component we quote the observed rotational velocity v_rot and not its intrinsic value $v_{\mathrm{rot}}^{\mathrm{intrinsic}}$, as the latter requires knowledge of the inclination angle i of the angular momentum vector with respect to the observer, that is $v_{\mathrm{rot}}=v_{\mathrm{rot}}^{\mathrm{intrinsic}}\sin (i)$. The stellar component of LG dwarf galaxies is not in a thin disk; therefore its inclination cannot be determined without knowledge of the intrinsic 3D shape. In those systems that contain a clearly rotating HI component, as IC 1613, one can at least assume the inclination determined from the neutral gas component, under the assumptions that stars and gas share the same angular momentum vector. This is what we will do when quoting values for $v_{\mathrm{rot}}^{\mathrm{intrinsic}}$, or when deriving estimates for the circular velocity, using preferentially the value quoted in Table 1. In Sect. 6.2, we discuss the possibility that IC 1613 may be seen close to face-on (Read et al. 2016), and suggest a new value based on the combined kinematic properties of the different stellar tracers we analysed.

4.1. Kinematic properties of the main sample: detection of a rotation signal

For the kinematic analysis of IC 1613, the free parameters to be fitted were the systemic l.o.s. velocity ${\overline{v}}_{\mathrm{sys}}$, the l.o.s. velocity dispersion σ_v, and the linear gradient k for the rotation model. In the latter case, the particular spatial distribution of the MUSE pointings and the limited area they cover prevent us from finding θ. Therefore, we fixed the central coordinates and θ of the kinematic field to be either that of the stellar component or that of the HI gas (see values in Table 1).

We assumed the following prior ranges: $-50<({\overline{v}}_{\mathrm{sys}}-v_{\mathrm{glx}})[\mathrm{km}{\mathrm{s}}^{-1}]<50$, where v_glx is the mean value of the input velocity distribution, 0 < σ_v[ km s⁻¹]< 50, −10 < k[ km s⁻¹ arcmin⁻¹]< 10. The analysis also returned the model evidence Z which resulted useful for comparing the statistical significance of one model against another through the use of the Bayes factor ln(B_1, 2) = ln(Z₁/Z₂). We obtained the kinematic parameters and the membership probabilities using the MultiNest code (Feroz et al. 2009; Buchner et al. 2014), a multi-modal nested sampling algorithm.

We used the Besançon model of Galactic foreground stars (Robin et al. 2003) as a prior over the expected l.o.s. velocity distribution of the contaminants. We generated a catalogue along IC 1613’s direction, over an area up to its half-light radius, covering the range of colours and magnitudes of our targets. This resulted in a velocity distribution that was well sampled and representative of the spatial area around the MUSE pointings. We approximated it as the sum of two Gaussian profiles (with means ${\overline{v}}_{\mathrm{Bes},1}=-12$ km s⁻¹, ${\overline{v}}_{\mathrm{Bes},2}=-66$ km s⁻¹, and standard deviations σ_Bes, 1 = 34 km s⁻¹, σ_Bes, 2 = 114 km s⁻¹, with an amplitude ratio of k₁/k₂ ∼ 2). According to this model, considering the small area covered by the MUSE pointings and the spanned range of magnitudes and colours, we expect to find, in total, around five foreground contaminants in the velocity range between −400 km s⁻¹ and 200 km s⁻¹. If we restrict to the velocity range covered by the input targets (i.e. between −330 km s⁻¹ and −140 km s⁻¹, see Fig. 6), we would expect only one contaminant. We recall that already during the spectral classification step (see Sect. 3.1.1) we found and removed six likely contaminants, all being cold MS stars with velocities around ${\overline{v}}_{\mathrm{Bes},1}$.

Results of the probability-weighted analysis are reported in Table 3. Of the initial 751 targets, 746 have a P_Mi >  0.95. We found that rotation is strongly favoured over a dispersion-only model, with a stronger evidence for the case that assumes the central coordinates and orientation of the kinematic major axis from the HI. In this case, the l.o.s. systemic velocity of the model including a rotation term is −230.9 ± 0.7 km s⁻¹, the velocity dispersion 11.2 ± 0.4 km s⁻¹, and the velocity gradient was 1.1 ± 0.3 km s⁻¹ arcmin⁻¹ (see Fig. 7). As these values have the highest evidence, we adopt these as our reference values in the following analysis.

Fig. 7. MultiNest 2D and marginalised posterior probability distributions for the systemic velocity, velocity dispersion, and velocity gradient assuming a linear rotational model with the kinematic PA aligned with that of the HI kinematic major axis. Dashed lines in the histograms indicate the 16th, 50th, and 84th percentiles. Contours are shown at 1, 2, and 3σ level.

The systemic velocity and velocity dispersion values we found are within 1σ from those previously published by Kirby et al. (2014, i.e. −231.6 ± 1.2 km s⁻¹ and $10.8_{-0.9}^{+1.0}$ km s⁻¹, respectively). In contrast to this previous work, we detect a clear rotation signal. In particular, Wheeler et al. (2017), analysing data from Kirby et al. (2014), reported a rotation amplitude for IC 1613 that was largely unconstrained and consistent with a null value. This is not surprising as these data are spatially distributed along the kinematic minor axis (i.e. roughly perpendicular to the direction of our MUSE dataset). Therefore, this is the first time that stellar rotation is detected with high significance in IC 1613.

We note that the five low probability targets could still belong to IC 1613 based on their spectral type. Only one has P_Mi ≈ 0, but we checked that the velocity determination in this case was affected by the presence of residual emission lines from the surrounding ionised gas. By avoiding these lines, we recalculated its velocity, which was compatible with ${\overline{v}}_{\mathrm{sys}}$, although with too large an error (δ_v ∼ 30 km s⁻¹). The other stars have 0.3 < P_Mi < 0.95. Three of these are low S/N young stars with large velocity errors. Interestingly, the remaining star has P_Mi ∼ 0.8 and v_i = −184 ± 6 km s⁻¹, which puts it more than 3σ away from the expectation of the best fitting rotation model. However, it is found in F1 and has photometric properties and metallicity (see Sect. 5) in common with the RGB stars of the galaxy. Given the low probability of being a MW contaminant, this star could be a possible binary, but we cannot exclude that it could belong in projection to an extended hot stellar halo around IC 1613. Indeed, there is growing evidence that MW satellites have extended stellar halos (e.g. Chiti et al. 2021; Longeard et al. 2023; Sestito et al. 2023a,b; Jensen et al. 2024). Additionally, Pucha et al. (2019) have shown that the evolved stellar population of IC 1613 extends to 4 × R_e.

We continue our analysis by selecting sources with a P_Mi >  0.95, creating a high fidelity sample without the need to recalculate individual membership probabilities. This allows us to explore the kinematic properties of different subsamples and make a detailed comparison with the HI kinematic field.

4.2. Kinematic properties of each pointing

Here we determine the kinematic properties of IC 1613’s stellar component as a function of radius by analysing the l.o.s. velocities of probable members within each MUSE pointing independently using a dispersion-only model. Results are reported in Table 4 and shown in Fig. 8, where we also compare them to the kinematic properties of the HI component, derived in Read et al. (2016) from the Little-THINGS survey data (Hunter et al. 2012; Oh et al. 2015). To make a direct comparison between gas and stars, we use the same central coordinates, in this case those of the HI, as there is a slight offset with the optical values (see Table 1)¹². We verified that the bias introduced by the rotation in the recovery of the velocity dispersion values in each field is negligible in our case (see also Leaman et al. 2012).

Fig. 8. Average rotational velocity (left) and velocity dispersion (right) values for the full stellar sample (top) and dividing by spectral type (bottom) projected along the HI kinematic semi-major axis. Black circles represent the values obtained for the probable members per pointing, while coloured symbols indicate values obtained from the subsamples of red giant (red squares) and main sequence stars (blue triangles); also shown are the rotational velocity and dispersion values for the HI (green pentagons; Read et al. 2016). The circular velocity values for both stars and the HI are also shown (light black circles and light green pentagons, respectively), correcting for the asymmetric drift and assuming an inclination angle of iHI = 39.4°. The grey line represents the rotation curve from the linear model adopted in Sect. 4, with the grey bands indicating the 95% confidence interval. The lightgreen line is an exponential fit to the dispersion profiles of a set of dwarf galaxies with high-resolution HI data having similar mass as IC 1613 and radii normalised to the scale radius Rd (Iorio et al. 2017; Mancera Piña et al. 2021, 2024); the lightgreen bands indicates a typical uncertainty of 2.5 km s−1. For the stellar component, the error bars in the radial direction indicate that the rotational velocities are averages within the extent of each MUSE pointing.

It is evident that the typical velocities of stars in each field ${\overline{v}}_{\mathrm{Field}}$ change with respect to the systemic velocity ${\overline{v}}_{\mathrm{sys}}$ when moving outwards from the galaxy’s centre (see the left panels of Fig. 8, where we define the rotational velocity as v_rot as ${\overline{v}}_{\mathrm{sys}}-{\overline{v}}_{\mathrm{Field}}$ in order to show a positive trend). F1 and F2 contain almost 90% of the inspected stars and thus dominate the linear rotation signal recovered in the analysis of the entire sample (shown as a black solid line in figure). Indeed, the velocity difference between them is ∼4 km s⁻¹ which, considering that the average distance between F1 and F2 is ∼4 arcmin, is in agreement with the observed velocity gradient of k ∼ 1 km s⁻¹ arcmin⁻¹.

From Fig. 8, we see that F2 and F3 show approximately the same rotational velocity (since their ${\overline{v}}_{\mathrm{F}2}\approx {\overline{v}}_{\mathrm{F}3}\approx -235$ km s⁻¹). Due to the limited spatial sampling, it is difficult to determine whether the stellar velocity signal has reached a constant value or whether it will continue to vary radially. We recall that the only available spectroscopic dataset that we could use to improve the spatial sampling is that published by Kirby et al. (2014) which, unluckily, is mainly distributed along the optical minor axis of the galaxy providing, as already shown by Wheeler et al. (2017), a poor constraint to the rotation signal. A similar behaviour at these radii is also seen in the HI component, where the local depressions of the rotation curve may be due to the presence of large HI holes, particularly evident around 0.5 and 1.5 kpc from the centre, as shown in Fig. 1 and also discussed in Read et al. (2016) and Collins & Read (2022).

Inspection of the velocity dispersion of the stellar component reveals a radial decrease, with the σ_v of the central pointing being close to the value obtained from the full sample, while for F2 and F3 we find a significant decrease to a roughly similar lower value. The velocity dispersion profile of the HI starts at 3 − 4 km s⁻¹ near the centre of the galaxy, then rises up to ∼9 km s⁻¹ at R ∼ 0.5 R_e (i.e. around the position of F2), to finally decrease and stabilise around ∼5 km s⁻¹ for R ∼ R_e. The velocity dispersion of the HI is systematically lower than that obtained for the stars at the same radii, with this difference being much more marked in F1. Again, it is likely that the sudden drop in the central part of the HI velocity dispersion profile (marked with hollow pentagons in Fig. 8) is due to the much lower gas density in the inner region, and thus the presence of a shell-like structure within R ≲ 0.5 kpc; see again Fig. 1, but also Moiseev & Lozinskaya 2012; Stilp et al. 2013). The HI velocity dispersion is indeed correlated to the gas surface mass density, which regulates the amount of local turbulence (e.g. Bacchini et al. 2020).

To corroborate the idea that the sudden drop in velocity dispersion for the HI is related to the HI hole, we also overplotted in Fig. 8 (upper-right panel) an exponential fit to the dispersion profiles of a set of dwarf galaxies with high-resolution HI data having similar mass as IC 1613 (i.e. having v_circ < 35 km s⁻¹), as modelled in Iorio et al. (2017) and Mancera Piña et al. (2021, 2024). In fact, the dispersion profile is expected to increase exponentially towards the centre, and it seems that this was also the case for IC 1613 before the formation of the HI hole.

In order to estimate the galaxy’s circular velocity v_circ and verify that we can directly compare the kinematics of stars and gas, it is necessary to account for the inclination and remove the random motion component that suppresses the rotation curve. We followed the formalism of Read et al. (2016, see further details in Sect. 6.2) to apply the so-called asymmetric drift correction, assuming the same inclination angle (see Table 1) for both tracers. After applying this correction, the stellar and HI components resulted in excellent agreement, as shown in Fig. 8 (upper left panel). We discuss further in Sect. 6.2 the caveats related to the inclination value and their impact on the circular velocity and dynamical mass estimations.

4.3. Kinematic properties as a function of stellar type

We now divide the main sample of probable member stars into two subsamples, young MS stars and evolved red giant stars, to examine the kinematic properties as a function of stellar type. Relying on the analysis in Sect. 3.1.1, we perform a simple selection according to the measured T_eff: red giant stars (RGS, i.e. belonging to the RGB and AGB) as those with T_eff ≤ 7500 K, and MS stars (including the Be stars) as those with T_eff >  7500 K. We verified that the T_eff limit is in general effective in separating the two samples, comparing with the T_eff expected from the (V − I) colours of our targets using the empirical calibration for red giant stars from Alonso et al. (1999), valid in the colour range of 0.8 < (V − I) < 2.2 and mostly metallicity independent. We note that these two samples are also distinct in age. Based on stellar evolution theory, it is known that the RGS track stars of ages > 1 − 1.5 Gyr. While a simple comparison with a set of isochrones (with [Fe/H] = [ − 1.3; −0.7] dex, Girardi et al. 2000; Bressan et al. 2012) tells us that the stars we classify as MS are between 35 < t_age[Myr]< 560.

Results of the kinematic analysis are reported in Table 4 and Fig. 8. The outcomes for the RGS are very similar (within 1σ) to those of the full sample in terms of rotation velocities and velocity dispersion per field, as well as the global values of the linear velocity gradient. This is expected, as this subsample represents a large fraction (∼85%) of the main one. As is evident from Fig. 8, a model allowing for a flat rotation curve, or a more steeply rising rotation curve than the one we adopted, would have been a much better match to the kinematics of the MS stars. At all radii, the MS stars are showing a higher amplitude in their rotation with respect to the RGS (and the full sample), and in F2 and F3 also a significantly lower velocity dispersion. From the same figure, we can see that the MS stars follow the HI kinematic properties very closely (apart from the velocity dispersion in F1), while on the other hand, the kinematics of the older RGS, have decoupled from those of the HI. Therefore, the kinematics of the young MS stars appear to be still coupled to that of the HI component, although the large errors in the external pointings weaken this conclusion.

We have verified that the reported kinematic values are generally robust against the presence of possible biases, such as the presence of variable or binary stars, as reported in Appendix B.3. Only for the MS stars in F1, we found that the high fraction of OB stars present could inflate the value of the velocity dispersion towards a higher value. This would partly explain why this value is so close, in absolute terms, to that of the RGS. Nevertheless, it is insufficient to account for the large discrepancy with the HI central velocity dispersion. The simplest explanation for this discrepancy is again related to the presence of the central HI hole, where F1 is located. Extrapolating the HI velocity dispersion from the outer parts towards the centre on the basis of the trend observed in the literature for galaxies of similar mass as IC 1613 (see Fig. 8), we would expect the HI velocity dispersion at the location of the central pointing to be very similar to that measured for the RGS and MS stars. Therefore it is possible that the dispersion and spatial distribution of the MS stars in the central pointing reflects that of the HI gas from which they were born. Such a HI hole could have been caused through stellar feedback mechanisms like a supernova (SN) explosion. The shape and nature of the central ionised shells shown in Fig. 2 are typical of a SN remnant, which seem to support this explanation. In addition, the spatial distribution of the OB stars inside and just outside F1 (Fig. 5) seems to suggest not only a physical association with the ionised shells, but that they may have originated from the compression of cold gas due to SN shock waves (Borissova et al. 2004; Garcia et al. 2009).

5. Metallicity properties

The chemical analysis focuses only on the RGB stars in our sample, because, due to the low resolution of the MUSE data, it is difficult to make a quantitative estimate of spectral parameters other than T_eff using the spectral fitting techniques implemented in this work (see details in Appendix B.1). Therefore, we determined metallicities ([Fe/H]) for individual RGB stars from the equivalent widths of the near-IR Ca II triplet (CaT) lines. We restrict our analysis to this subsample by selecting stars with P_Mi ≥ 0.95 and 3500 < T_eff[K]< 6000, which helps to separate them from the hotter evolving giants. We verified the goodness of this selection using a set of isochrones (Bressan et al. 2012) with Z = 0.0001 ([Fe/H] ∼ −2.3 dex) and t_age >  1 Gyr, which traced the lower range of metallicities and stellar ages for the RGB stars as obtained from the SFH analysis of IC 1613 (Skillman et al. 2014). We also selected all stars with S/N_CaT >  10 to keep the uncertainties on the equivalent widths (EWs) less than 20%. This resulted in a subsample of 275 input sources.

We made use of the Starkenburg et al. (2010) calibration. As shown in Kacharov et al. (2017), this calibration can be safely applied to data of spectral resolution as low as R ∼ 2600, comparable to that of our MUSE spectra around the CaT lines. The Starkenburg et al. (2010) calibration combines the EWs of the two reddest Ca II lines with the (V − V_HB) term, where V is the visual magnitude of the selected star and V_HB that of the horizontal branch. The EWs were obtained by fitting a Voigt profile over a 15 Å window around the selected Ca II lines, weighting each pixel value for the flux uncertainty stored in the error spectrum. The EW uncertainties were then calculated from the covariance matrix of the fitted Voigt parameters. We adopted a V_HB value of 24.91 ± 0.12 (Bernard et al. 2010). Final uncertainties for the metallicity values were obtained by error propagation of the EW uncertainties. We verified that the input magnitudes’ errors do not have a significant impact on the final [Fe/H] uncertainties.

The uncertainties obtained in the EW measurements resulted to be about 0.5 Å, while the average [Fe/H] error was ∼0.25 dex. However, a few targets (< 5%) have metallicity values with errors δ_[Fe/H] >  0.5 dex. We note that a small fraction of the stars in our subsample are likely to be AGB stars, brighter on average, for which the calibration method can be applied anyway, since it does not introduce a significant systematic bias compared to RGB stars (see discussion in Pont et al. 2004).

From the analysis of the metallicity values we obtained a median [Fe/H] = −1.06 dex, a median absolute deviation σ_MAD = 0.29 dex, and an intrinsic scatter σ_intrinsic = 0.26 dex (following Eq. 8 of Kirby et al. 2013), as also reported in Table 5. For comparison, we also applied the Carrera et al. (2013) calibration using the fitting windows defined by Cenarro et al. (2001). The main difference with the Starkenburg et al. (2010) calibration is that this one is empirical, whereas the former also used synthetic spectra. Nevertheless, we obtained comparable results: median [Fe/H] = −1.14 dex, σ_MAD = 0.37 dex and σ_intrinsic = 0.34 dex (see again Table 5). The recovered median metallicity is ∼0.1 dex lower, but the offset is within the scatter of the data, in this case ∼0.1 dex higher.

Our results are in good agreement with those reported by Kirby et al. (2013): median [Fe/H] = −1.22 dex and σ_MAD = 0.23 dex, also obtained targeting the RGB stars of the galaxy. The small deviations we found, again of the order of 0.1 dex, can be attributed to the different spectral resolution between the data and implemented technique to obtain the metallicity values (namely by directly fitting the available Fe lines).

We further analysed our sample by considering each pointing separately. In Fig. 9, we show the metallicity distribution per field, for clarity only plotting values obtained with the Starkenburg et al. (2010) calibration. Visual inspection of the MDFs could give the impression of a decrease in the presence of metal-rich stars from F1 to F3 fields. Therefore, we ran a two-sample two-sided Kolmogorov–Smirnov test, comparing separately F2 and F3 against F1, which is the largest sample and the one with the widest range of metallicities. We found that it is likely that the samples are drawn from the same parent distribution (p-values higher than 0.05); therefore, the lack of metal-rich stars in F2 and F3 is due to the lower statistics. This is consistent with IC 1613 constantly forming stars in a spatially homogeneous manner (see Skillman et al. 2014), at least out to R_e, which is how far our data extend.

Fig. 9. Metallicity distributions. Left: Histograms of the metallicity values for the full sample and for each pointing coloured in purple, blue, green, and yellow, respectively. Right: [Fe/H] values as a function of the elliptical radius scaled with Re, represented as black dots. The black solid line represents the result of a Gaussian process regression analysis using a Gaussian kernel and taking into account an intrinsic scatter; the grey band indicates the corresponding 1σ confidence interval. The histogram on the right side represents the metallicity distribution of the full sample.

We found that the median [Fe/H] tends to decrease very slightly from the central to the outer field, as reported in Table 5 for both calibrations used (see also the right panel of Fig. 9). Thus, we investigated the presence of a radial metallicity gradient performing a Gaussian process regression (GPR) analysis using a Gaussian kernel together with a noise component to account for the intrinsic metallicity scatter of the data. Details of this type of analysis have been described extensively in Taibi et al. (2022). For both calibrations we obtained a value of the gradient ∇_[Fe/H] that is consistent with zero within the uncertainties (see Fig. 9 and Table 5), meaning that the bulk of the data, with [Fe/H] values of ∼ − 1.0 dex, do not show significant spatial variation. This result is consistent with the conclusions of Taibi et al. (2022), who found that LG dwarf galaxies, except for those likely to have undergone past mergers, have similar gradient values.

The CaT method only applies to RGB stars, while we have mentioned that SPEXXY provides spectral parameters, including [Fe/H], for all stars in our sample. However, we have shown in the Appendix B.1 that it has limitations in providing correct estimates. Nevertheless, we were able to show qualitatively that the RGB stars should be at least 0.2 dex more metal-poor than the MS stars, as expected from the SFH of IC 1613 (Skillman et al. 2014) and spectroscopic measurements of young supergiants (Bresolin et al. 2007; Tautvaišienė et al. 2007; Berger et al. 2018).

6. Discussion

In the previous sections, we have provided the results of the kinematic and chemical analysis of the MUSE data of IC 1613. The general picture we get is that both the gaseous and stellar components are rotating, with the evolved stars appearing to decouple from the motion of the HI mainly due to the asymmetric drift. The MS stars, on the other hand, appear to have maintained kinematics similar to that of the neutral gas. The velocity dispersion profile decreases with radius for all tracers, although the central MS value may be affected by the presence of binaries or the expansion motion of the central hole, which may bias it towards a higher value. The lack of neutral gas in the central parts is also the probable cause of the sudden decrease in the inner HI velocity dispersion. Nevertheless, there seems to be a correlation between the velocity dispersion of the different tracers and their age, with values increasing as we move from the gas towards the evolved stars at all radii. At the same time, the evolved stars seem to have a lower rotation support than the HI and MS stars.

In the following, we examine these age-kinematic relations in detail, discussing their implications for the formation mechanisms of IC 1613. At the same time, we address the fact that the inclination of this galaxy is rather uncertain, proposing an alternative way of calculating it using only stellar tracers of different ages. Although we do not obtain conclusive results, we find the multi-tracer analysis to be a promising way to resolve this issue.

6.1. Age-kinematics relations

In Sects. 4.2 and 4.3, we have shown how the different stellar tracers inspected exhibit different kinematic behaviour. Here we make more explicit the connection of these changes with stellar age, also taking into account information recovered from the chemical analysis. In particular, as shown in Fig. 10, we inspect the variation of the velocity dispersion and of the rotation-to-dispersion support ratio as a function of age.

Fig. 10. Age-kinematics relations of different tracers. From left to right, from top to bottom, are reported the σv (first three panels) and vrot/σv (last panel) for the different tracers analysed in this work, as a function of their average age (covered ranges are reported in the legend). Markers as in Fig. 8, to which we add as orange and maroon squares the values for the MR and MP subsamples of RGB stars, respectively.

In the first three panels of Fig. 10 (from left to right, from top to bottom), we show the age-velocity dispersion variation of the different kinematic tracers. We note that each point has an age value that is only the average of the probable age range covered by the individual tracers. To provide a finer sampling of the age axis, we added two additional subsamples selected from the RGB stars with metallicity measurements. We used their median [Fe/H] value (see Table 5) to divide them into a metal-poor (MP) and a metal-rich (MR) subsample. According to the age-metallicity relation obtained from the SFH of IC 1613 (Skillman et al. 2014), they cover an age range between 1 < t_age, MP[Gyr]< 5.5 and 3.5 < t_age, MR[Gyr]< 13.5, respectively. In all three panels, the σ_v values of each tracer are reported in Table 4, while for the HI we took the average of σ_v, HI over the entire radial range (top left panel), the average for R < 3 arcmin (top right panel) and the average for R >  3 arcmin (bottom left panel). The average of the σ_v, HI profile measured for other galaxies in the literature is also shown for comparison.

From the figure, it can be seen that the velocity dispersion tends to increase as we move towards higher ages, both when considering each tracer in its entirety as well as when considering the central and external pointings separately. In particular, we see that in the top two panels the stellar σ_v is higher than σ_v, HI, but increases only slightly towards older ages, being consistent with a constant trend when considering the errors. In the bottom-left panel, however, we see a clear linear increase, with the σ_v of the younger tracers being significantly lower than σ_v, MP. Considering the presence of the central HI hole, whose causes may have affected not only the kinematics of the HI itself but also that of the MS stars, the external pointings probably provide a more accurate picture of the age-velocity dispersion relation of IC 1613. In the bottom-right panel of Fig. 10, we show instead the ratio of rotational to dispersion support (without accounting for the inclination) for both stars and gas as a function of age. The HI value is calculated at the effective radius R_e, while the stellar values are the average of the external F2 and F3 pointings (see again Table 4). Again, a linear trend is evident, with v_rot/σ_v decreasing with age, with values between ∼0.5 − 1.2.

Such age-kinematics trends are in good agreement with those of two similarly bright, gas-rich, and isolated LG dwarf galaxies, WLM and NGC 6822. Both systems show a v_rot/σ_v that decreases with age and a σ_v that grows with age (Leaman et al. 2012; Belland et al. 2020)¹³. In particular, an age-velocity dispersion relation (AVR) has also been observed at higher luminosities for the LMC, M 33, M 31 and the MW, while the satellite dSphs of the MW seem to show a constant trend with age (e.g. Leaman et al. 2017, but also Dorman et al. 2015; Quirk et al. 2019, 2022; and e.g. Poci et al. 2019; Shetty et al. 2020; Pessa et al. 2023, for AVR studies outside the LG). The AVR has historically been interpreted by assuming that stars are born from a gas that is always in a dynamically cold state, with the increase in velocity dispersion with age due to various scattering mechanisms (e.g. giant molecular clouds or spiral arms) that secularly increase the orbital energy (Spitzer & Schwarzschild 1951, 1953; Carlberg et al. 1985; Kokubo & Ida 1992). Alternatively, the AVR could be interpreted by assuming that stars are born from an interstellar medium that becomes less turbulent over time due to a decreased gas fraction and less gravitational instability in the gas disk (e.g. Bournaud et al. 2009; Forbes et al. 2012). The analysis of Leaman et al. (2017) favours this latter scenario, in which a latent dynamical heating also seems to be at play moving towards higher galactic masses, while environmental effects may have affected the MW satellite systems. Although the coarse age resolution of our sample prevents us from performing a more quantitative analysis of the AVR of IC 1613, our results seem to support the conclusions of Leaman et al. (2017).

A similar argument could be made inspecting the variation of the rotation support with age. We have seen that the v_rot/σ_v(R_e) of IC 1613 decrease with age, similarly to WLM and NGC 6822. However, few other systems in the LG have kinematic studies covering the youngest stars, with the majority focusing on the more accessible red giant population. In this context, Wheeler et al. (2017) showed with an homogeneous kinematic analysis that in general the gas-rich systems of the LG are not strong rotators in their evolved stellar component (i.e. they have v_rot/σ_v ≲ 2). These results appear valid regardless of the environment, in agreement with recent results found for dwarf galaxies residing in voids (de los Reyes et al. 2023). The average rotation-to-dispersion support of the evolved stars in IC 1613 is v_rot/σ_v(R_e)∼0.5, which remains ≲2 even taking into account the uncertainties over the inclination value (see Sect. 6.2). This is in agreement with the results of Wheeler et al. (2017) who concluded, comparing with zoom-in cosmological simulations of isolated dwarf galaxies, that dwarf galaxies probably form as puffy dispersion-supported systems, rather than with a cold kinematics. Considering that IC 1613 evolved mostly in isolation, the age trends observed in v_rot/σ_v and σ_v appear to be intrinsic features of this system, the interpretation of which favours the Leaman et al. (2017) and Wheeler et al. (2017) scenarios.

The age-kinematic trends of IC 1613 are also in broad agreement with observations of galaxies at different redshift (z), where it can be seen that in general σ_v increases, while v_rot/σ_v decreases, as a function of z (e.g. Wisnioski et al. 2015, 2019; Übler et al. 2019). These studies, however, barely cover the mass regime of dwarf galaxies, while recent results seem to question whether high redshift galaxies are dispersion-supported, with a v_rot/σ_v similar to that of local spirals (Rizzo et al. 2020, 2021; Fraternali et al. 2021; Lelli et al. 2021). Such systems are more massive than local dwarf galaxies, and at present only a few facilities such as the James Webb space telescope, can help to understand if rotation-supported dwarfs were common in the early Universe (e.g. de Graaff et al. 2024). However, considering the kinematic state of local dwarf galaxies in isolation, including gas-poor ones (Taibi et al. 2018, 2020), this possibility seems less likely (Leaman et al. 2017; Wheeler et al. 2017).

Our results on IC 1613 show the importance of performing a multi-tracer chemo-kinematic analysis covering several age ranges. Ongoing and future spectroscopic surveys, such as DESI (Levi et al. 2013), WEAVE (Jin et al. 2024) and 4MOST (de Jong et al. 2019), will provide such an analysis for the nearest dwarf galaxies with large statistical samples. At the distance of IC 1613 and beyond, we have shown here the challenges inherent in obtaining large samples covering a wide-area, for which we should await the arrival of planned or proposed facilities, such as the Mauna Kea Spectroscopic Explorer (MSE, The MSE Science Team 2019) and the Wide-field Spectroscopic Telescope (WST, Mainieri et al. 2024).

6.2. Dynamical mass estimation and the role of inclination

We previously mentioned that the HI inclination value of IC 1613 is a matter of debate. Using low-resolution HI data, Lake & Skillman (1989) reported a value of i = 38° ±5°, in very good agreement with the best-fitting value from Read et al. (2016, see Table 1), which we have adopted. A similar value was also obtained from Oh et al. (2015), although unconstrained. However, Read et al. (2016) suggested that the true inclination could be as low as i = 15°, which allowed the galaxy to be brought into agreement with model expectations. In fact, the HI circular velocity curve of IC 1613 grows very slowly, even compared to expectations of a cored dark matter density profile and assuming the most favourable value of low inclination. This appears to be the consequence of depressions and holes in the HI field, with the stellar field probably having been affected by them too, as we saw in Sects. 4.2 and 4.3. For Collins & Read (2022) an inclination value of ∼20° would bring the galaxy in agreement with expectations from the baryonic Tully–Fisher relation (BTFR). We note that such a low inclination is at odds with Lake & Skillman (1989), who argue that a higher inclination value is to be preferred on the basis that the external isophotes of the HI (assumed as a thin disk) agree well in shape and orientation with those of the stellar component, as is also the case between the spatial and kinematic PA of the HI. Indeed from the stellar photometric axis ratio we would obtain an inclination of i = arccos(b/a) = arccos(1 − ϵ)≈37°. If the galaxy was seen almost face-on, it would be difficult to explain such features.

Here we attempt to derive the inclination angle entirely from the stellar component, using the RGS and MS stars as independent dynamical tracers, whose properties have been reported in Sect. 4.3. The idea is quite simple: since both tracers must share the same circular velocity, we can equalise the formula for the asymmetric drift correction for both of them and solve for the inclination. In practise, we equalise the right-hand side of the following equation: $v_{\text{circ}}^2=v_{\text{rot}}^2/{\text{sin}}^2(i)+{\sigma }_{\text{D}}^2$, where σ_D is the asymmetric drift correction term necessary to remove the random motion component that suppresses the rotation curve. Another implicit assumption is that both tracers also share the same inclination angle. The σ_D term could be rather complex, as it takes into account the shape of the velocity ellipsoid and the tracer density profile (see e.g. Weijmans et al. 2008).

As a first attempt, we follow the formalism of Read et al. (2016) which primarily assumes a kinematic tracer with an exponential surface brightness profile and a constant observed velocity dispersion profile. This reduces the above equation to: $v_{\text{circ}}^2\simeq v_{\text{rot}}^2/{\text{sin}}^2(i)+{\sigma }_{\text{v}}^{2}R/R_{\text{d}}$, where R_d is the exponential scale radius. We further assume that both tracers share the same scale radius (as R_d = R_e/1.68). To solve for the inclination, we draw 10⁵ random values from the distributions of the rotation and velocity dispersion values, assuming Gaussian errors. We find then an inclination of i = 36° ±11°. This value is dominated by the observational errors, in particular of the MS stars.

We also attempted a more sophisticated approach following, for instance, the prescription of Weijmans et al. (2008, using their Eqs. 6–12 and following considerations in their Appendix A.2). In particular, assuming that the galaxy is axisymmetric, that the rotation curve flattens at large radii and that the potential is spherical, the asymmetric drift correction term becomes:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{D}}^2=-{\sigma }_{\mathrm{R}}^2[\frac{\partial \ln \mathrm{\Sigma }}{\partial \ln R}+\frac{\partial \ln {\sigma }_{\mathrm{R}}^2}{\partial \ln R}+\frac{R^2}{2R_{\mathrm{mod}}^2}+\frac{\kappa R^2}{\kappa (2R_{\mathrm{mod}}^2-R^2)+R^2}]\end{array}$$

where the first term between parenthesis is the slope of the surface brightness profile, again assumed to be exponential, while the last two terms are connected to the shape of the velocity ellipsoid and the shape of the rotation curve. The former is modelled with a parameter κ ∈ [0; 1], while the latter assumes that the rotation curve follows a power-law model where the curve becomes flat at large radii. The radius at which the rotation curve starts to flatten is called core radius R_c and enters in the equation as $R_{\text{mod}}^2=R^2+R_{\text{c}}^2$. We calculate R_c from the HI rotation curve (7.2 ± 1.9 arcmin) and assume it is the same also for the stellar tracers. While calculating i, we also randomise over these two terms assuming a uniform distribution for κ and a Gaussian one for R_c. Nevertheless, their contribution to the final value resulted to be minor. The radial dispersion is instead parametrised as σ_R(R) = σ₀ + σ₁exp(−R_mod/R_d), which can be obtained from the observed σ_v following Eq. (11) of Weijmans et al. (2008). This equation depends on i, but also on κ and R_c. Therefore, in order to find a parametric solution of σ_R, we also randomise over i that will act as a prior distribution, which we assume uniform over u ∈ [0; 1], with i = arccos(u). Compared to the simpler approach described above, we find here a slightly lower inclination value, i = 32° ±11°, but a comparable uncertainty.

Both recovered inclinations lie in the middle of the literature values, with the error bars preventing one from favouring a high (Lake & Skillman 1989) or low value (Read et al. 2016; Collins & Read 2022). It seems though that a value as low as 15° is less favoured. If an inclination of ∼35° will be preferred, this will confirm IC 1613 as an outlier of the BTFR. Such behaviour seems unusual for a galaxy with a baryonic mass of M_b = 9.5 ± 0.5 × 10⁷ M_⊙ (Read et al. 2016), while it seems to resemble that of gas-rich ultra-diffuse galaxies found in the field (though an order of magnitude more massive; see Mancera Piña et al. 2019). We defer to a future dynamical study the task of better quantifying the value of inclination, where all available kinematic tracers, both gas and stars, are simultaneously taken into account (see e.g. Leung et al. 2021, for WLM).

We can now estimate IC 1613’s dynamical mass, based on the stellar rotational velocities, corrected for the galaxy’s inclination, and the asymmetric drift; therefore turned into circular velocities v_circ. For the inclination, we stick to the last calculated value of 32° ±11°, while for simplicity we assume again the asymmetric drift correction from Read et al. (2016) in order to estimate v_circ. We perform this estimate within the 3D de-projected half-light radius r_1/2 ≃ (4/3)R_e to make a direct comparison with values from mass estimators for dispersion supported systems. We take the average of the observed values of v_rot and σ_v between F2 and F3, assuming that both the rotation and dispersion profiles are slowly varying up to the half-light radius. The circular velocity is then v_circ(r_1/2) = 15 ± 2 km s⁻¹. The dynamical mass $M_{\text{dyn}}(r_{1/2})=v_{\text{circ}}^2{r}_{1/2}/G$ is then 1.1 ± 0.3 × 10⁸ M_⊙, while the corresponding mass-to-light ratio within the half-light radius would be M_dyn/L_V(r_1/2) = 2.2 ± 0.6 M_⊙/L_⊙.

For comparison, the Wolf et al. (2010) mass-estimator for dispersion-supported systems, that is $M_{1/2}=3G^{-1}{\sigma }_{\mathrm{v}}^2{r}_{1/2}\approx 4G^{-1}{\sigma }_{\mathrm{v}}^2{R}_{\mathrm{e}}$, where G is the gravitational constant, yields M_1/2 = 1.0 ± 0.2 × 10⁸ M_⊙ within the same radius. This is within the previous results accounting for inclination, but also within the value reported by Kirby et al. (2014). We note that this mass estimator is valid as long as the system is spherically symmetric, without rotation and with a flat l.o.s. velocity dispersion profile. To some degree, all these conditions are violated for IC 1613.

7. Summary and future work

In this paper we present results from the chemo-kinematic analysis of the stellar component of the isolated dwarf irregular galaxy IC 1613. The analysis is based on a new set of spectroscopic data obtained with three pointings of the VLT/MUSE instrument. We extracted ∼2000 sources, from which we separated the stellar objects for their subsequent spectral classification and analysis.

The quality of the dataset allowed for a generally accurate spectral classification (i.e. with T_eff determined to better than 500 K) of 808 stars. We found a majority of K-type stars, and representatives from all the other types (from O-type to M-type stars). The classification included the identification of a sample of 24 probable Be stars, together with a sample of 14 probable C stars.

For all types except the Be stars, the l.o.s. velocities were determined with a pixel-by-pixel analysis, resulting in measurements with a mean precision of δ_v = 7 km s⁻¹ and signal-to-noise ratio S/N_C ∼ 10. For the Be stars we have determined the velocities by cross-correlation with dedicated templates. From this sample we identified 746 probable members (i.e. having P_M >  0.95), which show a significant rotation pattern. When modelled as a linear velocity gradient this yields k = 1.1 ± 0.3 km s⁻¹ arcmin⁻¹, along the HI kinematic major axis (with PA = 72°). This is the first time that rotation for the stellar component has been detected with high significance for this galaxy.

We further analysed the probable members looking at the kinematic properties of each individual field and as a function of stellar types tracing different broad age ranges. The individual fields analysis confirmed the presence of a rotation signal, and a velocity dispersion profile that decreases with radius. On average, the stellar component shows higher velocity dispersion and a slower rotation than the HI gas; this is a manifestation of the asymmetric drift, as demonstrated by the fact that the kinematics of the two tracers comes in good agreement when applying an asymmetric drift correction. The kinematic analysis of the different selected stellar populations showed a coherent picture with the analysis of the main sample, with the MS stars following more closely the kinematics of the HI gas than the evolved red giant stars.

The chemical analysis was conducted on the subsample of RGB stars exploiting the wavelength region of the Ca II triplet lines. We obtained a median [Fe/H] = −1.06 dex, with scatter values of σ_MAD = 0.29 dex and σ_intrinsic = 0.26 dex, in good agreement with previous results reported by Kirby et al. (2013). The average [Fe/H] is well within the rms scatter of the stellar luminosity-metallicity relation found for the LG dwarf galaxies by the same authors. The spatial distribution of the [Fe/H] values showed no significant signs of a radial metallicity gradient, compatible with results from other LG dwarf galaxies (Leaman et al. 2013; Taibi et al. 2022).

Examining the kinematics of stars as a function of broad age ranges, we find that the velocity dispersion increases as a function of age, with the behaviour being very clear in the outermost pointings, while the rotation-to-velocity dispersion support decreases. On timescales of < 1 Gyr, the stellar kinematics still follow very closely that of the neutral gas, while the two components decouple on longer timescales. We interpret these results as a consequence of stars being born from a less turbulent gas over time. On the other hand, with the dynamical mass estimation we provide a new value i = 32° ±11° of the inclination angle of IC 1613. To this end, we use stellar tracers with different observed kinematics that must be consistent with the same underlying gravitational potential.

In future efforts, the data presented here could be used to perform a dedicated analysis of the massive stars in our sample using the FASTWIND stellar atmosphere code (Puls et al. 2005), which will allow us to obtain reliable stellar parameters from their spectra. We also aim to improve the identification process of the M-type stars in our sample, and therefore to constrain the spectroscopic C/M ratio. This will tell us more of the chemical properties of the intermediate-age stellar population of IC 1613. This MUSE dataset could also be used to analyse the ionised gas medium using available emission lines, like Hα or nebular oxygen lines in order to get an insight into the kinematics of the shell-like gas structure visible in field F1 and compare it with that of the youngest stars. Finally, we aim to perform a dedicated dynamical study where all available kinematic tracers, both gas and stars, are simultaneously taken into account in order to firmly establish the (dark) matter distribution of IC 1613.

Acknowledgments

We wish to thank the anonymous referee for the constructive comments that helped to improve the manuscript. We also thank M.-R. Cioni, L. R. Patrick, J. Read, F. Lelli, and I. Trujillo for useful discussions and comments at various stage of this project. ST acknowledges funding of a Leibniz-Junior Research Group (PI: M. Pawlowski; project number J94/2020) via the Leibniz Competition. GB and CG acknowledge support from the Agencia Estatal de Investigación del Ministerio de Ciencia en Innovación (AEI-MICIN) and the European Regional Development Fund (ERDF) under Grant Number PID2020-118778GB-I00/10.13039/50110001103 and the AEI under Grant Number CEX2019-000920-S. GI acknowledges financial support from the European Research Council (ERC) through the ERC Consolidator Grant DEMOBLACK, under contract No. 770017. SK acknowledges funding from UKRI in the form of a Future Leaders Fellowship (grant no. MR/T022868/1). PEMP acknowledges the support from the Dutch Research Council (NWO) through the Veni grant VI.Veni.222.364. This research has made use of NASA’s Astrophysics Data System, VizieR catalogue access tool (CDS, Strasbourg, France, DOI: 10.26093/cds/vizier), and extensive use of Python3.8 (Van Rossum & Drake 2009), including iPython (v8.12, Pérez & Granger 2007), Numpy (v1.24, Harris et al. 2020), Scipy (v1.10, Virtanen et al. 2020), Matplotlib (v3.7, Hunter 2007), Astropy (v5.2, Astropy Collaboration 2022) and Scikit-learn (v1.2, Pedregosa et al. 2011) packages.

References

Appendix A: Calibration of MUSE catalogues

The raw catalogues obtained from the MUSE VRI pseudo-images contain the spatial and photometric information of all extracted sources, which is needed for the kinematic and chemical analysis we conducted. However, they needed astrometric and photometric calibration. For this purpose, we used the publicly available HST/WFPC2 catalogues from Holtzman et al. (2006), which corresponded to the F1 and F3 fields, and a proprietary Subaru/SuprimeCam catalogue for the F2 field, kindly provided to us by M. Monelli (IAC).

We note that the field of view (FoV) of the Subaru catalogue (34′×27′) is wide enough to cover all the three MUSE fields at once. However, it lacked an astrometric solution at first, and it showed severe crowding towards the galaxy’s centre, where F1 is placed. The HST catalogues, on the other hand, did not show any of these problems and we used them as a baseline during the astrometric and photometric calibration processes.

We downloaded the PSF-extracted HST catalogues with transformed magnitudes in the Johnson’s system¹⁴. The covered photometric bands were V and I, with the addition of the B-band for the central HST pointing. The Subaru catalogue covered instead the Johnson’s B and V bands. These are deep photometric catalogues, reaching down to the horizontal branch (at V ∼ 25), which means that they cover a much wider photometric range than that of our MUSE data. Before proceeding, we performed a pre-cleaning, selecting all sources marked as stellar in the HST catalogues, while retaining those targets in the Subaru catalogue having DAOPHOT parameters sharpness (SHARP) between −0.5 and 0.5, and a goodness-of-fit (CHI) < 1.

The astrometric solution and photometric calibration were obtained using the suite of codes CataXcorr and CataComb, kindly provided to us by P. Montegriffo and M. Bellazzini (INAF-OAS). We started with the astrometric registration of the Subaru catalogue. Due to its large FoV compared to that of the HST catalogues, we first performed a pre-astrometric registration using the stars in common with a PanSTARRS catalogue (Chambers et al. 2016) generated within the FoV of the Subaru catalogue. We note that the PanSTARRS catalogue was shallower that the Subaru one, reaching only a magnitude lower the tip of the RGB (at I ∼ 21). We then used the external HST field catalogue to refine the astrometric solution, improving its accuracy to better than 0.05″.

The photometric bands of the MUSE catalogues we were interested in calibrating were the V and I bands. We recall that the raw catalogues extracted from the data cubes were made in the VRI bands. The B-band was not extracted because it is not completely covered by the wavelength range of MUSE. Since the HST catalogues covered the BVI and VI bands (central and outer pointing, respectively), while the Subaru only covered the BV bands, as a compromise, we decided to add the I-band to the Subaru catalogue and calibrate the MUSE catalogues in VI.

We inferred the Subaru I-band using the common targets with the central HST catalogue (which covered the B-band). We applied the following linear equation: I_HST − V_HST = c_I (B_HST − V_HST)+ZP_I. We obtained a colour term c_I = −0.99 and a zero-point ZP_I = −0.08, which allowed to calculate the I_Sub magnitudes, whose agreement with the I_HST of the common targets was good, showing a median absolute deviation (MAD) scatter of 0.08 mags. However, comparing the (V − I) colours of the common sources, we found a difference of ∼0.2 mags driven by the bluer objects. This implied the introduction of a colour bias for the probable MS stars of the F2 field, as can be appreciated in Fig. 3, whose calibration was based on the Subaru catalogue. However, their calibrated (V − I) resulted safely below 0.5 mags (i.e. the colour limit of the RGB).

Finally, we used the mentioned catalogues to astrometrically and photometrically calibrate the corresponding MUSE fields – the HST ones with the F1 and F3 fields, while the Subaru one with the F2 field. For the photometric calibration of the MUSE pointings using the reference catalogues, we applied the following equations: V_MUSE − V_ref = c_V (V_ref − I_ref)+ZP_V, and I_MUSE − I_ref = c_I (V_ref − I_ref)+ZP_I. The accuracy of the astrometric and photometric solutions resulted on average 0.08″ and 0.05 mags, respectively.

Appendix B: Consistency checks

B.1. Verify SPEXXY general performance

We verified the performance of SPEXXY in recovering l.o.s. velocities and associated errors. We recall that SPEXXY performs a full spectral fit to the observed spectra using interpolated templates generated from the PHOENIX library of high-resolution synthetic spectra (Husser et al. 2013). It requires a configuration file specifying initial values for velocity, effective temperature, surface gravity, metallicity and [α/M] from which to start the fit. We obtained the initial stellar parameters from a set of isochrones (Girardi et al. 2000; Bressan et al. 2012) with [Fe/H] = −1 dex (around the mean metallicity of IC 1613) by performing a linear regression between the observed colours and magnitudes and the theoretical T_eff and log(g). We also set an initial velocity guess of −230 km s⁻¹ (i.e. approximately IC 1613 systemic velocity). For the spectral fitting, we used the available PHOENIX templates, convolved with the MUSE line spread function, with solar scale chemical composition, thus fixing [α/M] = 0 dex. We let the code find the best fit of all other parameters (i.e. T_eff, log(g) and [Fe/H], together with the l.o.s. velocity). We also masked certain spectral regions to avoid possible emission lines due to the interstellar medium and residuals due to sky line subtraction (such as the persistent residual due to the [OI] line at 5577 Å).

To test for possible systematic velocity shifts we ran the code on several noiseless templates from the PHOENIX library of giant stars at different T_eff (from 2500 K to 15000 K) at fixed [Fe/H] = −1.5 dex. Templates were shifted at several l.o.s. velocities, from −500 km s⁻¹ to 50 km s⁻¹ at step of 50 km s⁻¹. As initial guess for the velocity fitting, we set again −230 km s⁻¹. We found that in general SPEXXY correctly recovers the assigned l.o.s. velocities (mean and scatter of the difference between input and output velocities being −0.06 ± 0.56 km s⁻¹), although we noticed some failure around the extremes (i.e. for −500 km s⁻¹ and 50 km s⁻¹) with velocities variations of up to 50%. This occurred in only 3% of the cases that are those differing the most from the initial guess for the velocity. We recall that based on the Besançon model we expect almost no contaminants with true velocities ∼ − 500 km s⁻¹, while the variation at 50 km s⁻¹ is small enough to keep such possible contaminants away from the bulk velocity of IC 1613.

Regarding the other parameters, T_eff values were also well recovered, except for some scatter of ∼2000 K at temperatures > 10000 K; log(g) was recovered with a significant scatter (> 2 dex) for the hottest stars and those with T_eff ∼ 6500 K (the latter a minority in our main sample); [Fe/H] was also well recovered, but with a large scatter (> 1 dex) for templates with extreme velocities. Again, this is mostly due to the large difference between the initial guess and the true velocity of the stellar template, which leads SPEXXY to an incorrect parameter estimate. Therefore, SPEXXY generally recovers accurate velocities and spectral parameters.

To test whether SPEXXY provides the correct uncertainties, we developed a series of Monte Carlo tests in which we generated mock spectra at different S/N, starting from noiseless templates of different spectral types and injecting noise directly from the error spectra associated with the observed ones. However, we found that the flux uncertainty associated with the observed spectra was underestimated by an average factor of 1.4 with respect to the measured noise on the observed flux. This was found systematically for all spectra, regardless of their spectral type or S/N. It is most likely a result of the resampling into a regularly sampled data cube during data reduction. This introduces covariances between pixels that are neglected by the pipeline. In principle, this factor is not a problem, as the uncertainties of SPEXXY are calculated from the residuals between the input spectrum and the best-fit template, but we had to take this into account when generating the mock spectra.

Since we have two main stellar populations in our observed sample, RGS and MS stars, we generated mock spectra representative of these two subsamples. For the RGS, we selected stars with T_eff ∼ 5000 K, with a difference between observed values found by SPEXXY and ULySS of no more than 500 K. From their best-fitting templates and scaled noise spectra, we generated 250 mock spectra in bins of S/N_CaT around [3.5, 5, 10, 15, 25]. The same was done for the MS stars, but for stars with T_eff ∼ 10000 K, generating mock spectra at S/N₅₅₀ around [5, 10, 20, 50]. Results are reported in Table B.1.

We found that for the RGS stars their velocities and associated errors are well recovered, with deviations between the mean of the errors and the measurement scatter within 0.5 km s⁻¹. The effective temperature was also well recovered, although the errors were as much as ∼50 K lower than the measurement scatter. On the other hand, log(g) and [Fe/H] were more uncertain, with errors underestimated between 0.1 and 0.2 dex (higher deviations at lower S/N), leading to an offset on the recovered values of up to 0.1 dex.

For the MS stars we recovered higher deviations from the input values than the RGS case. In particular, velocity errors at S/N₅₅₀ = 5 were underestimated by more than 15 km s⁻¹, while at the highest S/N the deviation was less than 1 km s⁻¹. In the main text we reported that we made a quality cut in our sample by excluding stars with observed velocity errors > 25 km s⁻¹. This cut mainly excluded MS stars with the lowest S/N, while ensuring that we kept those with well-estimated velocity errors within 10%. For the spectral parameters, the associated errors were generally more underestimated (> 10%) than in the RGS case.

Therefore, it can be concluded that SPEXXY is generally effective in accurately and precisely recovering stellar velocities. However, it struggles to correctly recover stellar parameters, with log(g) being the most uncertain parameter in our mock tests. This issue is primarily related to the low resolution of the MUSE data rather than a problem with the fitting routine. It is also important to note that our tests only explored a limited region of the parameter space. To further evaluate the performance of SPEXXY in recovering [Fe/H] values, we conducted an additional test by comparing the values obtained from SPEXXY with those obtained for the RGB stars using the CaT method (see details in Sect. 5).

We found a significant deviation between the distributions of [Fe/H] values from the two methods of ∼0.4 dex calculated between the medians, with those obtained with SPEXXY being metal richer. The recovered deviation is significantly larger than expected based on the mock tests (see again Table B.1) and did not show a dependence with S/N. This discrepancy can be partly attributed to the poor recovery of log(g) by SPEXXY, as shown by the mock tests. In fact, we re-run SPEXXY on our targets while fixing log(g) to the initial input values, thus reducing the number of free parameters to fit. As a result, the offset to the CaT-based values decreased to ∼0.2 dex. The observed difference remains significant, but the distributions now have comparable standard deviations.

Fixing the value of log(g) during the SPEXXY fitting did not significantly change the [Fe/H] distribution of the MS stars. This could be due to the lower sensitivity of the MS stars to log(g) variations when recovering metallicity values. The results of the mock tests show that the recovered values are quite accurate (for S/N> 10), but with slightly underestimated errors.

Assuming that fixing the log(g) yields metallicity values accurate to within a constant shift in the RGB case, we can make a relative comparison between the [Fe/H] distribution of RGB and MS stars using SPEXXY outputs. We further refined the selection of MS stars by keeping the probable kinematic members with T_eff < 15000 K (i.e. up to the grid limit of the PHOENIX library) and S/N₅₅₀ >  15 to ensure a similar sample quality as for the RGB. This reduced the MS sample to 30 stars. Their average metallicity was comparable to that of the RGB sample, but with a larger associated standard deviation. Taking into account that the average value for the RGB may be overestimated by 0.2 dex, we can conclude that the MS stars are expected to be more metal-rich than the RGB. This difference may increase if we take into account the [α/Fe] for the two populations, which we set to zero during the SPEXXY fitting. The fact that the MS stars are more metal-rich than the RGB would be expected (Skillman et al. 2014), since also spectroscopic measurements of young blue and red supergiants from the literature suggest an average [Fe/H] = −0.7 dex (Bresolin et al. 2007; Tautvaišienė et al. 2007; Garcia et al. 2014; Berger et al. 2018), while the ionised HII medium may be slightly metal-rich (by ∼0.1 dex, e.g. Bresolin et al. 2007).

B.2. Comparison between SPEXXY and ULySS

In our analysis, we used SPEXXY and ULySS in a complementary way: the former for l.o.s. velocity measurements, while the latter for spectral classification. However, since both codes provide similar outputs, we compared their performances in recovering the l.o.s. velocity and T_eff from the observed spectra.

We considered all targets in the sample used for the kinematic analysis (see Sect. 4). The l.o.s. velocity errors for the ULySS measurements were obtained from Monte Carlo simulations using several templates of giant stars at different intervals of temperature and S/N. We made this choice because the formal errors returned by ULySS were generally overestimated, without following any clear trend with the S/N or spectral type. For the Monte Carlo simulations we chose seven templates from the MIUSCAT library (see Sect. 3.1.1), with T_eff ranging from 3500 K to 20000 K. The simulations were performed by adding Poisson noise to each template to obtain simulated spectra at increasing S/N intervals, with values ranging from 4 to 30. We repeated each process 100 times. For each template, the recovered velocity scatter as a function of S/N was fitted by an exponential law. The ULySS measurements were then divided into effective temperature bins according to the templates considered, and the velocity errors assigned using the corresponding exponential fit.

First we compared the velocity differences for pair measurements obtained with ULySS and SPEXXY. We found that ∼75% of the measured velocities are in agreement within the combined uncertainties, obtained as $\sqrt{{\delta }_{\mathrm{rad},\mathrm{Spexxy}}^2+{\delta }_{\mathrm{rad},\mathrm{ULySS}}^2}$. The comparison is shown in the left panel of Fig. B.1 where for clarity we reported only those measurements with a combined error < 30 km s⁻¹. Nevertheless, the difference between the ULySS and SPEXXY measurements showed an average positive shift of ∼5 km s⁻¹. What is more, a higher shift of ∼10 km s⁻¹ is observed for those objects with a S/N_C < 10. This was somehow expected for the UlySS measurements, since the Monte Carlo tests showed average velocity shifts of up to 20 km s⁻¹ for templates with effective temperatures of ∼4500 K and greater than ∼10000 K having simulated S/N < 10. On the other hand, we have already showed in the previous section the capacity of SPEXXY to correctly recover l.o.s. velocities, even at the lowest S/N for red giant stars. This is a very important point because most of the stars in our dataset are RGS in a low S/N regime (≲10), where SPEXXY showed its strength in the recovery of l.o.s. velocities.

Fig. B.1. Comparison between SPEXXY and ULySS. Left: Velocity differences for pair measurements obtained with the two codes. The solid black line represents the zero velocity offset. Filled circles are colour-coded according to their S/NC, while the grey error-bars indicate the combined uncertainty, that is $\sqrt{{\delta }_{\mathrm{rad},\mathrm{Spexxy}}^2+{\delta }_{\mathrm{rad},\mathrm{ULySS}}^2}$. Right: Direct comparison between the effective temperature values obtained with the two codes; the solid black line is for reference only.

The comparison of effective temperatures, as shown in the right-hand panel of Fig. B.1, gave good results for T_eff < 7500 K, although the SPEXXY values tend to be systematically ∼500 K higher than the ULySS values. At higher T_eff, SPEXXY values started to deviate significantly. This is mainly due to the fact that SPEXXY measurements are based on the PHOENIX library, which is not optimised for very hot stars. In the ULySS case, instead, for each spectrum the stored temperature value is that of the best fitting template, which was chosen by eye among the most likely templates. We recall that this procedure had an accuracy of up to ∼500 K, for a quality PCL-flag between 3 and 4.

B.3. Biases in the recovery of kinematic properties

We checked for the presence of possible biases that could alter the recovery of the kinematic properties of the stellar samples, in particular when inspecting the properties of the RGS and MS stars. On one hand, the region approximately separating the MS from the RGS, 0.4 ≲ (V − I)≲0.8 (see Fig. 4), contains the instability strip of IC 1613, where it is expected to find mainly classical Cepheids with short-term variability (i.e. with log(P[days]) < 0.5 and t_age < 0.5 Gyr; see e.g. Udalski et al. 2001; Bernard et al. 2010). Since we observe them in random phases, we expect any differences to be averaged out. It is also true that, especially for Cepheids, hydrogen or metallic lines can lead to different velocities (see e.g. Vinko et al. 1998). Nonetheless, we find no significant differences when calculating the average velocities and velocity dispersion including or excluding stars in that colour range. Similar considerations can be made for the Be-star candidates in our sample, which should also be short-period variables (e.g. Porter & Rivinius 2003; Rivinius et al. 2013), although in this case it is mostly the intensity of the emission lines that is affected. We repeated the same test finding values compatible within the errors. As for the RGS sample, this contains a small number of C-star candidates, which are long period variables (e.g. Menzies et al. 2015). However, we have verified that the C-stars have similar properties to the rest of the sample within the errors.

We did not check for the presence of binary stars, which could particularly affect the velocity dispersion values. The low S/N of the individual exposures for each field in our sample and the short period over which they were observed (between a few days and a month), preclude a quantitative analysis. We can, however, make qualitative estimates. For a sample of mostly RGB stars with high measured velocity dispersion (> 10 km s⁻¹), as it is our case, the impact of unresolved binaries is expected to be minimal (i.e. ∼10 − 20% of inflation; see Minor et al. 2010; Minor 2013; Spencer et al. 2017, 2018; Arroyo-Polonio et al. 2023). However, the case is different when focusing on the MS stars sample. The multiplicity fraction (i.e. the fraction of stars having at least one companion) is generally expected to be ∼60% for A-types, which increase to 80%−100% for OB stars (e.g. Moe & Di Stefano 2017; Offner et al. 2023). In particular, these stars are often found in triples or higher order multiples. For Be stars such figures are less constrained, but mass accretion through binary interactions has been proposed as a possible channel for their formation (Pols et al. 1991; de Mink et al. 2013; Shao & Li 2014; Bodensteiner et al. 2020).

Considering that the multiple OBs cover a range of orbital periods, preferably a few days or a few decades (Moe & Di Stefano 2017), the stacking process to generate the spectra should already average the short-period ones. For the long-period ones (i.e. log(P[days]) >  3.5) we expect velocity variations of ∼3 − 10 km s⁻¹, which could translate into a velocity dispersion inflation of ≲20%−30% (following Minor et al. 2010, assuming that they behave as short-period low-mass binaries and have intrinsic σ_v, 0 = 7 − 10 km s⁻¹). This would especially affect the central pointing, where the fraction of OB stars is high (∼70%). It could thus partly explain why in this field the value of the velocity dispersion of the MS stars is comparable to that of the RGS, while in the other fields it tends to be closer to that of the HI. Another possible source of disturbance for the MS stars in F1 is that they may still be coupled to the motion of the ionised bubble, which is probably expanding (with a typical velocity of ∼10 − 20 km s⁻¹, Pokhrel et al. 2020). We also do not exclude the presence of possible runaway OB stars that could further impact the kinematics of the sample (e.g. Stoop et al. 2024).

Appendix C: Observing log and examples of extracted spectra

We report here a table with the observing log of our VLT/MUSE observations of IC 1613. We show as well some examples of the extracted spectra analysed in this work.

Fig. C.1. Example of sky-subtracted, wavelength calibrated, extracted spectra from our dataset. In red the best-fitting templates provided by SPEXXY, with the output spectral parameters (identification, Vrad, Teff, log(g), [Fe/H], S/N, binary quality flag) indicated in the labels. From top to bottom examples of an RGB star; a hot MS star; a foreground star.

Fig. C.2. Figure C.1. Continued. From top to bottom examples of an M giant star; a C star; a Be star.

Fig. C.3. Figure C.1. Continued. From top to bottom examples of a background galaxy; a blend with an high redshift emitting-galaxy; an example of ionised gas in IC 1613.

All Tables

All Figures

Fig. 1. Finding chart showing the location of the three MUSE fields, marked as white boxes, overlaid over a wide-field image of IC 1613 (Credit: ESO – VST/Omegacam Local Group Survey). The white ellipse marks the half-light radius. The smoothed HI column density map from the Little-THINGS survey (Hunter et al. 2012) is marked with logarithmically spaced isodensity contours (silver lines) starting at 3σ (white-smoke lines). North is up, east to the left.

In the text

	Fig. 2. Broad- and narrow-band images obtained from the MUSE data cubes. Top row: Combined images in the VRI-bands highlighting stellar sources and background galaxies. Bottom row: Combined images highlighting the emission-line ionised gas with Hα coloured in red, [SII] (6713 Å) in blue, and [OIII] (5007 Å) in green. In all panels, North is up and east is to the left.
In the text

Fig. 3. Extinction-reddening corrected colour-magnitude diagrams from MUSE images of the three fields, moving from the centre towards the outer regions from left to right; the filled circles represent the sources for which spectra have been extracted and they are colour-coded by their S/N at the central wavelength. We note that the galaxy’s distance modulus is (m − M)0 = 24.4 (Bernard et al. 2010).

In the text

	Fig. 4. Colour-magnitude diagram for data selected with PCL-flag between 3 and 4, as described in Sect. 3.1.1. Filled circles are colour-coded according to their effective temperature Teff as derived with ULySS; black triangles indicate the identified C-stars, while the grey squares mark the Be stars.
In the text

	Fig. 5. Spatial distribution of hot OB stars (blue filled circles) and Be stars (darkslate grey filled squares), compared to that of the OB-associations (grey filled hexagons) reported by Garcia et al. (2009). The black contours indicate the Hα emission within the MUSE pointings, marked by black squares.
In the text

	Fig. 6. Histogram of the l.o.s. velocity measurements divided by pointing (left panel) and stellar type (right panel).
In the text

	Fig. 7. MultiNest 2D and marginalised posterior probability distributions for the systemic velocity, velocity dispersion, and velocity gradient assuming a linear rotational model with the kinematic PA aligned with that of the HI kinematic major axis. Dashed lines in the histograms indicate the 16th, 50th, and 84th percentiles. Contours are shown at 1, 2, and 3σ level.
In the text

Fig. 8. Average rotational velocity (left) and velocity dispersion (right) values for the full stellar sample (top) and dividing by spectral type (bottom) projected along the HI kinematic semi-major axis. Black circles represent the values obtained for the probable members per pointing, while coloured symbols indicate values obtained from the subsamples of red giant (red squares) and main sequence stars (blue triangles); also shown are the rotational velocity and dispersion values for the HI (green pentagons; Read et al. 2016). The circular velocity values for both stars and the HI are also shown (light black circles and light green pentagons, respectively), correcting for the asymmetric drift and assuming an inclination angle of iHI = 39.4°. The grey line represents the rotation curve from the linear model adopted in Sect. 4, with the grey bands indicating the 95% confidence interval. The lightgreen line is an exponential fit to the dispersion profiles of a set of dwarf galaxies with high-resolution HI data having similar mass as IC 1613 and radii normalised to the scale radius Rd (Iorio et al. 2017; Mancera Piña et al. 2021, 2024); the lightgreen bands indicates a typical uncertainty of 2.5 km s−1. For the stellar component, the error bars in the radial direction indicate that the rotational velocities are averages within the extent of each MUSE pointing.

In the text

Fig. 9. Metallicity distributions. Left: Histograms of the metallicity values for the full sample and for each pointing coloured in purple, blue, green, and yellow, respectively. Right: [Fe/H] values as a function of the elliptical radius scaled with Re, represented as black dots. The black solid line represents the result of a Gaussian process regression analysis using a Gaussian kernel and taking into account an intrinsic scatter; the grey band indicates the corresponding 1σ confidence interval. The histogram on the right side represents the metallicity distribution of the full sample.

In the text

Fig. 10. Age-kinematics relations of different tracers. From left to right, from top to bottom, are reported the σv (first three panels) and vrot/σv (last panel) for the different tracers analysed in this work, as a function of their average age (covered ranges are reported in the legend). Markers as in Fig. 8, to which we add as orange and maroon squares the values for the MR and MP subsamples of RGB stars, respectively.

In the text

Fig. B.1. Comparison between SPEXXY and ULySS. Left: Velocity differences for pair measurements obtained with the two codes. The solid black line represents the zero velocity offset. Filled circles are colour-coded according to their S/NC, while the grey error-bars indicate the combined uncertainty, that is $\sqrt{{\delta }_{\mathrm{rad},\mathrm{Spexxy}}^2+{\delta }_{\mathrm{rad},\mathrm{ULySS}}^2}$. Right: Direct comparison between the effective temperature values obtained with the two codes; the solid black line is for reference only.

In the text

	Fig. C.1. Example of sky-subtracted, wavelength calibrated, extracted spectra from our dataset. In red the best-fitting templates provided by SPEXXY, with the output spectral parameters (identification, Vrad, Teff, log(g), [Fe/H], S/N, binary quality flag) indicated in the labels. From top to bottom examples of an RGB star; a hot MS star; a foreground star.
In the text

	Fig. C.2. Figure C.1. Continued. From top to bottom examples of an M giant star; a C star; a Be star.
In the text

	Fig. C.3. Figure C.1. Continued. From top to bottom examples of a background galaxy; a blend with an high redshift emitting-galaxy; an example of ionised gas in IC 1613.
In the text