11institutetext:  Leibniz-Institut für Astrophysik Potsdam (AIP), An der
Sternwarte 16, D-14482 Potsdam, Germany
11email<EMAIL_ADDRESS>22institutetext: Instituto de Astrofísica de Canarias,
Calle Vía Láctea s/n, E-38206 La Laguna, Tenerife, Spain 33institutetext:
Universidad de La Laguna, Avda. Astrofísico Fco. Sánchez, E-38205 La Laguna,
Tenerife, Spain 44institutetext: Universität Potsdam, Institut für Physik und
Astronomie, Karl-Liebknecht-Str. 24/25, D-14476 Potsdam, Germany
55institutetext: Dipartimento di Fisica e Astronomia “Galileo Galilei”,
Università di Padova, vicolo dell’Osservatorio 3, IT-35122, Padova, Italy
66institutetext: INAF - Osservatorio Astronomico di Padova, vicolo
dell’Osservatorio 5, IT-35122 Padova, Italy 77institutetext: INFN - Padova,
Via Marzolo 8, I–35131 Padova, Italy 88institutetext: Astrophysics Research
Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3
5RF, UK 99institutetext: Minnesota Institute for Astrophysics, University of
Minnesota, 116 Church St. SE, Minneapolis, MN 55455, USA 1010institutetext:
Department of Astrophysics, University of Vienna, Türkenschanzstraße 17, 1180
Vienna, Austria 1111institutetext: Centre for Astrophysics and Supercomputing,
Swinburne University, John Street, Hawthorn VIC 3122, Australia

# From gas to stars: MUSEings on the internal evolution of IC 1613

S. Taibi 11 G. Battaglia 2233 M. M. Roth 1144 G. Iorio 556677 S. Kamann 88 C.
Gallart 2233 E. D. Skillman 99 R. Leaman 1010 N. Kacharov 11 M. A. Beasley
11112233

(Received; accepted)

###### Abstract

Context. The kinematics and chemical composition of stellar populations of
different ages provide crucial information on the evolution of the various
components of a galaxy.

Aims. Our aim is to determine the kinematics of individual stars as a function
of age in IC 1613, a star-forming, gas-rich and isolated dwarf galaxy of the
Local Group.

Methods. We present results of a new spectroscopic survey of IC 1613 conducted
with MUSE, an integral field spectrograph mounted on the Very Large Telescope.
We extracted $\sim 2000$ sources, from which we separated stellar objects for
their subsequent spectral analysis. The quality of the dataset allowed us to
obtain accurate classifications ($T_{\rm eff}$ to better than $\pm 500$ K) and
line-of-sight velocities (with average $\delta_{v}\sim 7$ km s-1) for about
800 stars. Our sample includes not only red giant branch (RGB) and main
sequence (MS) stars, but also a number of probable Be and C stars. We also
obtained reliable metallicities ($\delta_{\rm[Fe/H]}\sim 0.25$ dex) for about
300 RGB stars.

Results. The kinematic analysis revealed for the first time the presence of
stellar rotation with high significance. We found general agreement with the
velocity field of the neutral gas component, although on average the stars
showed greater velocity dispersion and slower rotation, which can be
attributed to asymmetric drift. Examining the kinematics of the different
stellar components, the MS stars appear to follow closely that of the gas. The
chemical analysis of the RGB stars revealed average properties comparable to
other Local Group dwarf galaxies.

Conclusions. Our work represents a step forward in understanding the internal
processes that govern the dynamical evolution of a low-mass galaxy.

##  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 due 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 instability 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, 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, but
also).

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 soon 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)_{0}=24.400\pm 0.014$, or
$759\pm 5$ kpc.

IC 1613 therefore lies well beyond the Milky Way’s (MW) virial radius, as well
as that of M31 (being at $\sim 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 M31, for M31 masses of $\gtrsim 1.3\times
10^{12}\leavevmode\nobreak\ M_{\odot}$. 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 a mass of M31 of $2\times 10^{12}\leavevmode\nobreak\
M_{\odot}$. 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 M31.

Its low luminosity (an absolute magnitude of $M_{V}=-15.2$) corresponds to a
stellar mass of $M_{*}\sim 10^{8}\leavevmode\nobreak\ M_{\odot}$, 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; Bernard et al. 2007; Skillman et
al. 2014). 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 ${\rm[Fe/H]}\sim-2.0$ dex to
$-0.8$ dex. Stars were formed at an average rate of $\psi(t)=0.081\pm
0.001\leavevmode\nobreak\ M_{\odot}{\rm yr}^{-1}$.

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 $\sim 24$ arcmin
(i.e., $\sim 4$ effective radii $R_{e}$). In particular, the young stars are
found well beyond the currently active star forming regions (within $\sim
1.5\times 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 outward 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 an
homogeneous analysis for 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 also the ionized gas is
distributed (e.g., Lozinskaya et al. 2003; Silich et al. 2006; Moiseev &
Lozinskaya 2012; Pokhrel et al. 2020). The HI kinematics shows a linearly
increasing rotation curve with a $v_{\rm max}\sim 20$ km s-1 (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
face-on. This translates into a significant uncertainty in the determination
of the dynamical mass being between $M_{\rm dyn}\sim 5\times 10^{8}M_{\odot}$
and $\sim 8\times 10^{9}M_{\odot}$, 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 9 early B-type young supergiants
(Bresolin et al. 2007) and 3 evolved M-type supergiants (Tautvaišienė et al.
2007). Both studies obtained compatible metallicity values, $12+{\rm
log(O/H)}=7.90\pm 0.08$ dex and ${\rm[Fe/H]}=-0.67\pm 0.09$ dex,
respectively111Assuming ${\rm[O/Fe]}=-1$ dex and a solar oxygen abundance of
$12+{\rm log(O/H)}=8.69$ dex (Asplund et al. 2009)., 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. they have lower metallicity in
regions of high HI density). 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 ${\rm[\alpha/Fe]}$
ratio could be sub-solar (i.e., $\sim-0.1$ dex, Tautvaišienė et al. 2007 and
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 and 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 ${\rm[Fe/H]}=-1.19\pm 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\pm 1.2$ km s-1,
in agreement with values from the HI component (Lake & Skillman 1989; Oh et
al. 2015), and a velocity dispersion of $10.8^{+1.0}_{-0.9}$ km s-1. They did
not find signs of rotation for the stellar component; this was later confirmed
also by Wheeler et al. (2017) who re-analysed the Kirby et al. dataset. This
is not surprising considering that the spatial distribution of this dataset is
roughly along the optical minor axis of the galaxy, while in general it is
expected that the stellar rotation follows that of the gas, which in this case
rotates roughly along the optical major axis (Oh et al. 2015). The dynamical
mass at the half-light radius reported by Kirby et al. (2014) was of
$M_{1/2}=1.1\pm 0.2\times 10^{8}M_{\odot}$.

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). Taking advantage of the
unique capabilities of this instrument, that combines high spatial resolution
with a 1 arcmin2 field-of-view and a wide wavelength range at medium spectral
resolution (Bacon et al. 2014), we performed a kinematic analysis as a
function of stellar age and obtained metallicites for the largest
spectroscopic stellar sample obtained to date for this galaxy. IC 1613 is
located at high Galactic latitude ($b\sim 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. Section 5 shows results of the chemical analysis of the
RGB stars. Finally, Sect. 6 is dedicated to the summary and perspectives of
future work, while in the appendices we report details on the sanity checks
conducted on the spectral-fitting codes used during our analysis. The
parameters adopted for IC 1613 throughout the text are summarized in Table 1.

Table 1: Parameters adopted for IC 1613.

Parameter | Units | Value | Ref.
---|---|---|---
$\alpha_{\rm J2000}$ |  | $01^{h}04^{m}47.8^{s}$ | (1)
$\delta_{\rm J2000}$ |  | $+02^{\circ}07^{\prime}04^{\prime\prime}$ | (1)
$L_{\rm V}$ | $10^{6}L_{\odot}$ | $100\pm 15$ | (1)
$D_{\odot}$ | kpc | $759\pm 5$ | (2)
$R_{\rm e}$ | arcmin | $6.8\pm 0.05$ | (3)
| kpc | $1.50\pm 0.02$ | (3)
P. A. | deg | $90\pm 1$ | (3)
$\epsilon^{(a)}$ |  | 0.20$\pm$0.05 | (3)
$\alpha_{\rm J2000,HI}$ |  | $01^{h}04^{m}50.2^{s}$ | (4)
$\delta_{\rm J2000,HI}$ |  | $+02^{\circ}08^{\prime}26^{\prime\prime}$ | (4)
P. A.HI | deg | $72\pm 2$ | (4)
$i_{\rm HI}$ | deg | $39\pm 2$ | (4)
$\bar{v}_{\rm sys}$ | km s-1 | $-231\pm 1$ | (5)
$\sigma_{\rm v}$ | km s-1 | $11.2\pm 0.4$ | (5)
$\left\langle\left[{\rm Fe/H}\right]\right\rangle$ | dex | $-1.06$ | (5)
$\sigma_{\rm[Fe/H]}$ | dex | 0.26 | (5)
Conversion | pc arcmin-1 | 221 |
222The table lists: the coordinates of the galaxy’s optical centre; the
stellar luminosity in V-band; the heliocentric distance; the effective radius;
the position angle measured from north to east; the ellipticity
(${}^{(a)}\epsilon=1-b/a$); central coordinates, position angle and
inclination of the HI velocity field; the chemo-kinematic parameters obtained
in this work, i.e., the systemic velocity, the velocity dispersion, the median
metallicity, and the intrinsic metallicity scatter.
References. (1) McConnachie (2012); (2) Bernard et al. (2010); (3) Higgs et
al. (2021); (4) Read et al. (2016); (5) this work.

##  2 Observations and Data Reduction

Figure 1: Finding chart showing the location of the 3 MUSE fields, marked as
white boxes, overlaid over a wide-field image of IC 1613 (Credit: ESO -
VST/Omegacam Local Group Survey). The ellipse marks the half-light radius.
North is up, East to the left.

The data were acquired with VLT/MUSE in service mode333Under ESO programme
097.B-0373; PI: Battaglia.. MUSE is an integral field spectrograph with a
spatial sampling of 0.2 arcsec 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 $\sim 1\arcmin\times 1\arcmin$ and a
wavelength coverage between $4800-9300$ Å. We observed 3 fields (hereafter F1,
F2, F3, moving outward 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 on-going star formation.

Each field was observed with a total exposure time of 11080 s on source, split
in $4\times 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.1). The data were acquired in dark time, clear sky, with a request for a
seeing in V-band at zenith $\leq 0.9$ arcsec, fulfilled for all but one
exposure of F1 (for which it deviated within 10% from the request). The
observing log is reported in Table 7.

###  2.1 Data Reduction and Source Extraction

Figure 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.

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 settings444This included correcting for the
observer’s line of sight velocity with respect to the barycentre of the solar
system.. 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-1. 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\times 300\times 3680$ pixels.

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 color-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. 11, 12 and 13. For each extracted spectrum, PampelMUSE provided the
associated flux uncertainty per pixel and an S/N value calculated around its
central wavelength (hereafter S/NC). However, as our sources cover a wide
spectral range (see Sect. 3.1), we calculated two further S/N indicators: one
around 5500 Å (S/N550), and the other on the continuum around the CaT lines
(S/NCaT). 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${}_{550}\geq 20$,
and about 200, 120, 45 sources with S/N${}_{\rm CaT}\geq 10$. We note that for
this galaxy dust reddening and extinction along the line of sight 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.

Figure 3: Extinction-reddening corrected colour-magnitude diagrams from MUSE
images of the 3 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 line-of-sight
(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).

###  3.1 Spectral classification with ULySS

The spectral type classification was carried out with the same approach as in
Roth et al. (2018); we refer the reader to this article for full details. A
linear combination of up to 10 templates from the MIUSCAT stellar spectral
library555The MIUSCAT stellar spectral library fills in the gap between the
MILES (Sánchez-Blázquez et al. 2006) and the near-IR CaT library of Cenarro et
al. (2001), using the Indo-US spectral library of Valdes et al. (2004); it
also extends the wavelength coverage slightly blue-ward of MILES and red-ward
of the CaT library. MIUSCAT has been recently replaced by E-MILES (Vazdekis et
al. 2016). (Vazdekis et al. 2012) was fitted to the observed spectra using the
ULySS666http://ulyss.univ-lyon1.fr code (Koleva et al. 2009). We note that
ULySS is capable of identifying unresolved blended stars because the spectrum
is fitted with a linear combination of non-linear components, multiplied by a
polynomial continuum. In our case, no blends are found, not even in the
central field. This is mostly expected, given the density of sources visible
in Fig. 2, which was constructed from the MUSE data cubes.

We stored the information of the three spectral types having the largest
weights in the best fit. Occasionally, the fit returns implausible mixes of
spectral types (e.g., mixing young early-types with evolved late-types). 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)_{0}=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. However, 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${}_{\rm C}>2$ and a clean photometry (see Sect. 2.1). Main outputs of the
code were the best fitting templates, together with the associated spectral
parameters ($T_{\rm eff}$, log(g) and [Fe/H]) and l.o.s. velocity. We verified
by eye the ULySS outcomes, assigning several quality flags to each analysed
object. We evaluated 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. 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., $\lesssim
10$), as we verified in Appendix B.1 and B.2.

Table 2: Distribution of spectral type classification.

Sp. Type |  $T_{\rm eff}$ (K) | N
---|---|---
M | $2400-3900$ | 1
K | $3900-4900$ | 467
G | $4900-6000$ | 90
F | $6000-7500$ | 46
A | $7500-10000$ | 57
B | $10000-30000$ | 96
O | $>30000$ | 13
Be | $>10000$ | 24
C | $<3500$ | 14
Total |  | 808
777From inspection of ULySS outcomes, with PCL flag between 3 and 4, as
described in Sect. 3.1. In the table, C and Be indicate respectively the
identified number of carbon and B emission-line stars.

Figure 4: Colour-magnitude diagram for data selected with PCL-flag between 3
and 4, as described in Sect. 3.1. Filled circles are colour-coded according to
their effective temperature $T_{\rm eff}$; black triangles indicate the
identified C-stars, while the grey squares mark the Be stars.

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.1): 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.3). 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_{\rm eff}>10000$ 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 stars888The
spectroscopic $T_{\rm eff}$ for the red giant stars were found to be in fair
agreement with the empirical $T_{\rm eff}:(V-I)$ relation from Alonso et al.
(1999), which is valid in the colour range $0.8<(V-I)<2.2$ and mostly
metallicity independent.. 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. Fig. 4 shows the
position of the classified stars on the colour-magnitude diagram, colour-coded
according to their $T_{\rm eff}$, visually confirming the general goodness of
the spectral classification.

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_{\rm eff}>30000$ K) is a limiting factor for the spectral
classification of these stars. We have identified 13 stars with a $T_{\rm
eff}$ just above 30000 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_{\odot}$.

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_{\rm eff}>20000$ K and
also the Be stars. We see that they generally tend to be found where the OB
associations are. Their spatial distribution also coincides with 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, these evidences at
least confirm the goodness of our initial classification of many of them as
hot stars.

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$\dots$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
$\alpha$-elements abundances, will be possible with the future BlueMUSE
instrument (Richard et al. 2019).

Figure 5: Spatial distribution of hot OB stars (blue filled circles) and Be
stars (dark grey filled circles), compared to that of the OB-associations
(gray filled circles) reported by Garcia et al. (2009). The black squares mark
the MUSE pointings.

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 $\sim 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_{\rm 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 10 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 ($\sim 10\%$, Hamren et al. 2016, e.g.,).
Considering between nine 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 with spexxy

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

The determination of l.o.s. velocities was carried out using the
spexxy999https://github.com/thusser/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 spectra101010http://phoenix.astro.physik.uni-
goettingen.de/ (Husser et al. 2013). The library covers a wide wavelength
range, going from 500 Å to 5.5 $\mu$m, and stellar parameters: $2300<T_{\rm
eff}\,(\rm K)<15000$, $0.0<{\rm log}(g)\,(\rm dex)<+6.0$,
$-4.0<{\rm[Fe/H]\,(dex)}<+1.0$ and $-0.2<{\rm[\alpha/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_{\rm eff}>10000$ 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.1). We also excluded all those
targets that in the spectral classification step (see Sect. 3.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/NCaT $\sim 3.5$ for red giant stars
(S/N550 $>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_{\rm 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 $\delta_{\rm
v}>25$ km s-1, 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-1.
Performing a cross-correlation with the third data release of the Gaia
catalogue (Gaia Collaboration et al. 2023), we further confirmed their
foreground nature based on their proper motion values, significantly different
from 0 km s-1 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 $\delta_{\rm v}\sim 6.5$ km s-1, while the median S/NC resulted around
10.

###  3.3 Velocity determination for carbon and Be stars

The l.o.s. velocities of the C stars identified by spectral classification
(see Sect. 3.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 $\delta_{v}\sim
4$ km s-1 for a median S/N${}_{\rm CaT}\sim 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 et al. 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/N550. For the 24 Be stars analysed (14, 6 and 4 in fields F1, F2 and F3,
respectively), we recovered a median $\delta_{v}\sim 7$ km s-1 for a median
S/N${}_{550}\sim 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

Table 3: Parameters and evidences resulting from the Bayesian kinematic
analysis.

Sample | Nin | NP | Model | $\bar{v}_{\rm sys}$ | $\sigma_{\rm v}$ | k | Bayes factor
---|---|---|---|---|---|---|---
|  |  |  | (km s-1) | (km s-1) | (km s-1 arcmin-1) | ln($B_{\rm rot,disp}$)
Main | 751 | 746 | Linear rotation (HI) | $-230.9\pm 0.7$ | $11.2\pm 0.4$ | $1.1\pm 0.3$ | 4.5
|  |  | Linear rotation (optical) | $-232.2\pm 0.6$ | $11.2\pm 0.4$ | $1.2\pm 0.3$ | 3.5
|  |  | Dispersion only | $-233.3\pm 0.5$ | $11.3\pm 0.4$ |  |

111111The reported values of the kinematic parameters represent the median of
the corresponding marginalized posterior distributions, with 1-$\sigma$ errors
set as the confidence intervals around the central value enclosing $68\%$ of
each distributions. The values in the row labelled (HI) or (optical) where
determined assuming the central coordinates and orientation of the kinematic
major axis of the HI component or the stellar component, respectively.

Figure 7: MultiNest 2D and marginalized posterior probability distributions
for the for the systemic velocity, velocity dispersion, and velocity gradient
assuming a linear rotational model with the kinematic P. A. 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-$\sigma$
level.

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 component121212On the angular scales covered by the
MUSE sample, we find the perspective rotation due to the relative motion
between the galaxy and the Sun to play a negligible role with $\Delta v$ up to
$-0.2\pm 0.2$ km s-1 (following Eq. 6 of van de Ven et al. 2006 and assuming
the systemic proper motion from Battaglia et al. 2022). Velocity gradients
induced by tidal disturbances from the MW or M31 are also highly unlikely..

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_{\rm M_{i}}$ to the individual targets. This approach is
based on the expectation maximization technique presented in Walker et al.
(2009).

The $P_{\rm M_{i}}$ 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 form131313https://scikit-
learn.org/stable/modules/isotonic.html.

We used the Besançon model of Galactic foreground stars (Robin et al. 2003) to
obtain 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 $\bar{v}_{\rm Bes,1}=-12$
km s-1, $\bar{v}_{\rm Bes,2}=-66$ km s-1, and standard deviations $\sigma_{\rm
Bes,1}=34$ km s-1, $\sigma_{\rm Bes,2}=114$ km s-1, with an amplitude ratio of
$k_{1}/k_{2}\sim 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-1 and 200 km s-1. If we restrict to the
velocity range covered by the input targets (i.e. between $-330$ km s-1 and
$-140$ km s-1, see Fig. 6), we would expect only one contaminant. We recall
that already during the spectral classification step (see Sect. 3.1) we found
and removed six likely contaminants, all being cold MS stars with velocities
around $\bar{v}_{\rm Bes,1}$.

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:
$k\,R_{i}\,{\rm cos}(\theta-\theta_{i})$, where k is the velocity gradient,
$(R_{i},\theta_{i})$ are the angular distance from the galaxy’s centre (i.e.
$R_{i}=\sqrt{\xi_{i}^{2}+\eta_{i}^{2}}$, where ($\xi_{i},\eta_{i}$) are the
tangent plane’s coordinates) and position angle (P. A.) of the i-target, and
$\theta$ is the P. A. of the kinematic major axis (both P. A. measured from
North to East).

We note that, unless otherwise stated, for the stellar component we quote the
observed rotational velocity $v_{\rm rot}$ and not its intrinsic value $v_{\rm
rot}^{\rm intrinsic}$, as the latter requires knowledge of the inclination
angle i of the angular momentum vector with respect to the observer (i.e.,
$v_{\rm rot}=v_{\rm rot}^{\rm intrinsic}{\rm sin}(i)$). The stellar component
of Local Group dwarf galaxies is not in a thin disc, 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 inclination. This is what we
will do when quoting values for $v_{\rm rot}^{\rm intrinsic}$ or when deriving
estimates for the circular velocity. However, we note the possibility that IC
1613 may be seen close to face-on (Read et al. 2016), therefore even the HI
value has its limitations.

The particular spatial distribution of the MUSE pointings and the limited area
they cover prevent us from finding $\theta$. Therefore, we fixed the central
coordinates and $\theta$ of the kinematic field to be either that of the
stellar component or that of the HI gas (see values in Table 1). The free
parameters were then the systemic l.o.s. velocity $\bar{v}_{\rm sys}$, the
l.o.s. velocity dispersion $\sigma_{\rm v}$ common to both models, and the
linear gradient k. We assumed the following prior ranges: $-50<(\bar{v}_{\rm
sys}-v_{\rm glx})[{\rm km\,s^{-1}}]<50$, where $v_{\rm glx}$ is the mean value
of the input velocity distribution, $0<\sigma_{\rm v}[{\rm km\,s^{-1}}]<50$,
$-10<k[{\rm km\,s^{-1}\,arcmin^{-1}}]<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
${\rm ln}(B_{1,2})={\rm ln}(Z_{1}/Z_{2})$. 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.

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

The systemic velocity and velocity dispersion values we found are within
1-$\sigma$ from those previously published by Kirby et al. (2014, i.e.
$-231.6\pm 1.2$ km s-1 and $10.8^{+1.0}_{-0.9}$ km s-1, respectively). On the
other hand, Wheeler et al. (2017), analyzing the data by Kirby et al. (2014),
reported a rotation amplitude for IC 1613 which was largely unconstrained and
consistent with null. This is not surprising as their data are spatially
distributed roughly perpendicularly to the direction of our MUSE data-set, for
which we do detect a clear rotation signal. Therefore, this is the first time
that rotation of the stellar component 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_{\rm M_{i}}\approx 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
$\bar{v}_{\rm sys}$, although with too large an error ($\delta_{v}\sim 30$ km
s-1). The other stars have $0.3<P_{\rm M_{i}}<0.95$. Three of these are low
S/N young stars with large velocity errors. Interestingly, the remaining star
has $P_{\rm M_{i}}\sim 0.8$ and $v_{i}=-184\pm 6$ km s-1, which puts it more
than 3$\sigma$ 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, Pucha et al. (2019) have shown that its evolved
stellar population extends to $4\times R_{e}$, while 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).

We continue our analysis by selecting sources with a $P_{\rm M_{i}}>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 sub-samples and make a detailed comparison with the HI kinematic
field.

###  4.1 Kinematic properties along each pointing

Table 4: Results from the Bayesian kinematic analysis performed for each
pointing applying a dispersion-only model.

Sample | NP | $\bar{v}$ | $\sigma_{\rm v}$
---|---|---|---
|  |  | (km s-1) | (km s-1)
| F1 | 410 | $-231.3\pm 0.7$ | $12.1^{+0.6}_{-0.5}$
Main | F2 | 274 | $-236.2\pm 0.7$ | $9.0\pm 0.6$
| F3 | 62 | $-234.5\pm 1.3$ | $8.2^{+1.2}_{-1.0}$
| F1 | 331 | $-230.6\pm 0.7$ | $12.2\pm 0.6$
RGS | F2 | 230 | $-236.1\pm 0.8$ | $9.3\pm 0.6$
| F3 | 54 | $-234.4\pm 1.4$ | $8.6^{+1.2}_{-1.1}$
| F1 | 79 | $-234.9^{+1.6}_{-1.6}$ | $11.8^{+1.5}_{-1.3}$
MS | F2 | 44 | $-237.3\pm 1.6$ | $6.2^{+2.0}_{-1.6}$
| F3 | 8 | $-236.6^{+4.2}_{-2.8}$ | $5.2^{+5.2}_{-3.2}$
141414The values of the kinematic parameters and their associated
uncertainties are reported as in Table 3.

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)151515We note that throughout the text, velocities with
the sub-script rot are the observed rotational velocities, neither corrected
for the inclination ${\rm sin}(i)$ nor for the asymmetric drift; those with
the sub-script circ are circular velocities, thereby corrected for both
effects.. 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).

Figure 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 sub-samples of red giant
(red squares) and main sequence stars (blue triangles); also shown 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
$i_{\rm HI}=39.4^{\circ}$. The black line represents the rotation curve from
the model adopted in Sect. 4, with the grey bands indicating the 95%
confidence interval. 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 $\bar{v}_{\rm
Field}$ change with respect to the systemic velocity $\bar{v}_{\rm sys}$ when
moving outward from the galaxy’s centre (see the left panels of Fig. 8, where
we define the rotational velocity as $v_{\rm rot}$ as $\bar{v}_{\rm
sys}-\bar{v}_{\rm 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
$\sim 4$ km s-1 which, considering that the average distance between F1 and F2
is $\sim 4$ arcmin, is in agreement with the observed velocity gradient of
$k\sim 1$ km s-1 arcmin-1.

From Fig. 8, we see that F2 and F3 show approximately the same rotational
velocity (since their $\bar{v}_{\rm F2}\approx\bar{v}_{\rm F3}\approx-235$ km
s-1). Due to the limited spatial sampling, it is difficult to determine
whether the 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, and it appears to be only a local flattening
of the rotation curve, as it becomes clear when looking at the HI kinematics
over a larger radial range (see also Oh et al. 2015; Read et al. 2016).

Inspection of the velocity dispersion of the stellar component reveals a
radial decrease, with the $\sigma_{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-1 near the centre of the galaxy, then
rises up to $\sim 9$ km s-1 at $R\sim 0.5\leavevmode\nobreak\ R_{\rm e}$
(i.e., around the position of F2), to finally decrease and stabilise around
$\sim 5$ km s-1 for $R\sim R_{\rm e}$. The velocity dispersion of the HI is
systematically lower than that obtained for F2 and F3 at the same radii, while
this difference increases significantly when compared to the F1 value. It is
probable that the lower central HI density due to the presence of a shell-like
gas structure could have had an impact on the gas kinematics, and hence on the
velocity dispersion profile, inside $R\lesssim 0.5\leavevmode\nobreak\ R_{\rm
e}$ (see 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). We will
investigate further in the future whether this may be related to the greater
velocity dispersion measured from the stars in F1.

In order to estimate the galaxy’s circular velocity and directly compare the
kinematics of stars and gas, it was necessary to remove the random motion
component that suppresses the rotation curve. We followed the formalism of
Read et al. (2016, using their Eqs. 9 and 12) to apply the so-called
asymmetric drift correction. Although their analysis focused on the HI
component, for our purposes the same reasoning could be applied to the stellar
field. We also assumed the same inclination angle (see Table 1) for both
tracers. After applying this correction, the stellar and HI components showed
excellent agreement. This indicates that the differences observed in the
rotational velocities and velocity dispersion between the HI and stellar
components are indeed real.

We now examine the ratio of rotational to dispersion support for both stars
and gas in IC 1613. For the stellar sample, we obtain a $v_{\rm
rot}/\sigma_{\rm v}\sim 0.5-0.7$ at the projected half-light radius $R_{e}$,
depending on whether we adopt as $v_{\rm rot}$ the average of the values
between F2 and F3 or whether we extrapolate the values obtained for the linear
rotation model reported in Table 3. We can consider the values from these two
approaches as compatible with each other, taking into account the errors, as
shown in Fig. 8. This is the observed $v_{\rm rot}/\sigma_{\rm v}$ without
taking into account the inclination angle of IC 1613.

Using low-resolution HI data, Lake & Skillman (1989) reported a value of
$i=38^{\circ}\pm 5^{\circ}$, 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^{\circ}$, allowing them for a better dynamical modelling. This is
in contrast to Lake & Skillman (1989), who argue for an inclination angle of
$\sim 40^{\circ}$. This conclusion is based on the fact that the external
isophotes of the HI agree well in shape and orientation with those of the
stellar component, as is also the case between the spatial and kinematic P.A.
of the HI. If the galaxy was seen almost face-on, it would be difficult to
explain such features.

Depending on the inclination adopted, the stellar $v_{\rm rot}^{\rm
intrinsic}/\sigma_{\rm v}$ is between 0.8 and 1.1 for $i=39^{\circ}$ or
between 1.9 and 2.7 for $i=15^{\circ}$. The HI, for comparison, has an
observed $v_{\rm rot,HI}/\sigma_{\rm v}\sim 1.2$ at the same radius, which
rises to $2-4.5$ when the inclination angles are taken into account.

We can compare our results with those of two similarly bright LG dwarf
irregular galaxies, WLM and NGC 6822. Unlike IC 1613, these are seen almost
edge-on. WLM shows a stellar $v_{\rm rot}/\sigma_{\rm v}\sim 1$ at $R_{e}$,
with the rotation axis aligned with that of the HI (Leaman et al. 2012). In
contrast, the kinematics of NGC 6822 is more complex, as is the spatial
distribution of its tracers. In fact, the P. A. of the HI is roughly
perpendicular to the optical, and while both the stellar and HI components
show signs of rotation (Weldrake et al. 2003; Demers et al. 2006; Swan et al.
2016; Wheeler et al. 2017; Belland et al. 2020), it is not clear whether they
share the same kinematic axis. In any case, the stellar rotation signal seems
to be between $v_{\rm rot}/\sigma_{\rm v}\sim$ 0.4 and 1 (Swan et al. 2016;
Wheeler et al. 2017; Belland et al. 2020). For both systems, the HI rotation
resulted in much higher values with $v_{\rm rot,HI}/\sigma_{\rm v}\sim 5$
within $R_{e}$ for NGC 6822 (Weldrake et al. 2003), and $v_{\rm
rot,HI}/\sigma_{\rm v}\sim 4$ for WLM (Oh et al. 2015; Iorio et al. 2017),
similar to the IC 1613 case (assuming though the lowest inclination angle).
Our results seem to be in good agreement with the conclusions drawn by Wheeler
et al. (2017) that the gas-rich systems of the LG are not strong rotators in
their stellar components (i.e. they have $v_{\rm rot}/\sigma_{\rm v}\lesssim
2$).

###  4.2 Kinematic properties as a function of stellar type

We now divide our main sample of probable member stars into two sub-samples,
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, we perform a simple selection according to the measured $T_{\rm eff}$:
red giant stars (RGS, i.e., belonging to the RGB and AGB with stellar ages
$t_{\rm age}>1$ Gyr) as those with $T_{\rm eff}\leq 7500$ K, and MS stars
(having $t_{\rm age}<1$ Gyr) as those with $T_{\rm eff}>7500$ K. We verified
that the $T_{\rm eff}$ limit is in general effective in separating the two
samples, comparing with the $T_{\rm 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.

The results are reported in Table 4 and Fig. 8. The results for the RGS are
very similar (within 1-$\sigma$) 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 sub-sample
represents a large fraction 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.

We checked for the presence of possible biases that could alter the recovery
of the kinematic properties of the MS stars. On one hand, the region
approximately separating the MS from the RGS, $0.4\lesssim(V-I)\lesssim 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[{\rm days}])<0.5$ and $t_{\rm age}<0.5$ Gyr; Udalski et al. 2001;
Bernard et al. 2010, see e.g.). 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 color 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-1), as it is our case, the
impact of unresolved binaries is expected to be minimal (i.e. $\sim 10-20\%$
of inflation; Minor et al. 2010; Minor 2013; Spencer et al. 2017, 2018;
Arroyo-Polonio et al. 2023, see). 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 $\sim 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[{\rm days}])>3.5$) we expect velocity variations
of $\sim 3-10$ km s-1, which could translate into a velocity dispersion
inflation of $\lesssim 20\%-30\%$ (following Minor et al. 2010, assuming that
they behave as short-period low-mass binaries and have intrinsic
$\sigma_{v,0}=7-10$ km s-1). This would especially affect the central
pointing, where the fraction of OB stars is high ($\sim 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 $\sim 10-20$
km s-1, 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).

In summary, the different examined tracers point to a similar picture: both IC
1613’s gaseous and stellar components are rotating, with the older stars
decoupled from the HI motion mainly by the effect of asymmetric drift, while
on timescales of $<1$ Gyr, the stars have maintained a similar kinematics as
the neutral gas. The velocity dispersion profile decreases with radius for all
tracers, although in F1 the value of the MS stars may be affected by the
presence of binaries or by the expanding motion of the central hole, which may
have also affected the central 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 towards the older
stars. This is in qualitative agreement with the results of Leaman et al.
(2017), suggesting that stars are born from gas that becomes less turbulent
over time due to decreasing gas fractions.

###  4.3 Dynamical mass estimation

We can use the results from the previous sections to 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. We perform this estimate within the 3D de-projected half-light
radius $r_{1/2}\simeq(4/3)R_{\rm e}$ to make a direct comparison with values
from mass estimators for dispersion supported systems.

We took the average of the observed values of $v_{\rm rot}$ and $\sigma_{\rm
v}$ between F2 and F3, assuming that both the rotation and dispersion profiles
are rather flat up to the half-light radius. Assuming both values quoted in
Read et al. (2016) for the the galaxy’s inclination i, we would have $V_{\rm
c}(r_{1/2})=13.2\pm 1.0$ km s-1 and $20.5\pm 3.0$ km s-1, respectively. The
dynamical mass $M_{\rm c}(r_{1/2})$ is then between $0.8\pm 0.1\times
10^{8}M_{\odot}$ and $1.9\pm 0.6\times 10^{8}M_{\odot}$, respectively. The
corresponding mass-to-light ratio within the half-light radius $M_{\rm
c}/L_{V}(r_{1/2})$ would be between $1.6\pm 0.3\,M_{\odot}/L_{\odot}$ and
$3.8\pm 1.2\,M_{\odot}/L_{\odot}$.

For comparison, the Wolf et al. (2010) mass-estimator for dispersion-supported
systems, that is $M_{1/2}=3G^{-1}\sigma_{\rm v}^{2}r_{1/2}\approx
4G^{-1}\sigma_{\rm v}^{2}R_{\rm e}$, where G is the gravitational constant,
yields $M_{1/2}=1.0\pm 0.2\times 10^{8}M_{\odot}$ within the same radius. This
is within the previous results accounting for inclination, but also with 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.

It is important to note that the estimates of the dynamical mass, and
dynamical mass-to-light ratio, based on rotational velocities given in this
section should be considered as lower limits. The HI circular velocity curve
of IC 1613 is known to rise very slowly, even with respect to the expectations
of a cored dark matter density profile and assuming the more favourable
inclination of $i=15^{\circ}$. This appears to be the consequence of
depressions, particularly evident around 0.5 and 1.5 kpc, due to the presence
of large HI bubbles, as discussed in Read et al. (2016). If the stellar
velocities reflect a similar behaviour, this would imply that also the stellar
rotation curve suffers from these dips, thereby leading to an underestimate of
the circular velocity and of the dynamical mass.

##  5 Metallicity properties

Table 5: Results from the chemical analysis.

|  | Starkenburg et al. (2010) | Carrera et al. (2013)
---|---|---|---
Sample | N | ${\rm[Fe/H]}$ | $\sigma_{\rm MAD}$ | $\sigma_{\rm intrinsic}$ | $\nabla_{\rm[Fe/H]}$ | ${\rm[Fe/H]}$ | $\sigma_{\rm MAD}$ | $\sigma_{\rm intrinsic}$ | $\nabla_{\rm[Fe/H]}$
|  | (dex) | (dex) | (dex) | (dex $R_{\rm e}^{-1}$) | (dex) | (dex) | (dex) | (dex $R_{\rm e}^{-1}$)
Full | 275 | $-1.06$ | 0.29 | 0.26 | $-0.06\pm 0.08$ | $-1.14$ | 0.37 | 0.34 | $-0.05\pm 0.11$
F1 | 150 | $-1.03$ | 0.29 | 0.26 |  | $-1.09$ | 0.36 | 0.32 |
F2 | 93 | $-1.10$ | 0.31 | 0.26 |  | $-1.15$ | 0.35 | 0.36 |
F3 | 27 | $-1.15$ | 0.22 | 0.21 |  | $-1.22$ | 0.30 | 0.32 |

161616From left to right, columns are the considered sample, its size N, its
median metallicity together with its scaled median absolute deviation,
intrinsic scatter, and radial metallicity gradient. Values are given for both
calibration methods applied to obtain them.

Figure 9: Metallicity distributions. Left: histograms of the metallicity
values for each pointing colored in blue, green and red respectively for F1,
F2 and F3. Right: [Fe/H] values as a function of the elliptical radius scaled
with $R_{\rm e}$, 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-$\sigma$ confidence interval. The histogram on the right side
represents the metallicity distribution of the full sample.

We determined metallicities ([Fe/H]) for the individual stars from the
equivalent widths of the near-IR Ca II triplet (CaT) lines. Since the
available calibrations are developed for RGB stars, we will be restricting our
analysis to this sub-sample (i.e., $P_{\rm M_{i}}\geq 0.95$, with $3500<T_{\rm
eff}[{\rm K}]<6000$, which helps separating 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]$\sim-2.3$ dex) and $t_{\rm
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${}_{\rm CaT}>10$ to keep the
uncertainties on the equivalent widths (EWs) less than 20%. This resulted in a
sub-sample 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\sim 2600$, comparable to that of our MUSE spectra
around the CaT lines. The Starkenburg et al. calibration combines the EWs of
the two reddest Ca II lines with the $(V-V_{\rm HB})$ term, where V is the
visual magnitude of the selected star and $V_{\rm 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_{\rm HB}$ value of $24.91\pm 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 $\sim 0.25$ dex. However, a few targets
($<5\%$) have metallicity values with errors $\delta_{\rm[Fe/H]}>0.5$ dex. We
note that a small fraction of the stars in our sub-sample 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
${\rm[Fe/H]}=-1.06$ dex, a median absolute deviation $\sigma_{\rm MAD}=0.29$
dex, and an intrinsic scatter $\sigma_{\rm 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. calibration is that this one is empirical, whereas the former also used
synthetic spectra. Nevertheless, we obtained comparable results: median
${\rm[Fe/H]}=-1.14$ dex, $\sigma_{\rm MAD}=0.37$ dex and $\sigma_{\rm
intrinsic}=0.34$ dex (see again Table 5). The recovered median metallicity is
$\sim 0.1$ dex lower, but the offset is within the scatter of the data, in
this case $\sim 0.1$ dex higher.

Our results are in good agreement with those reported by Kirby et al. (2013):
median ${\rm[Fe/H]}=-1.22$ dex and $\sigma_{\rm 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. 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_{\rm e}$, which is how
far our data extend.

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 $\nabla_{\rm[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 $\sim-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 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 3 pointings of the
VLT/MUSE instrument. We extracted $\sim 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_{\rm 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- to M-). 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 have been determined
with a pixel-by-pixel analysis, resulting in measurements with a mean
precision of $\delta_{\rm v}\pm 7$ km-1 and S/N${}_{\rm C}\sim 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 $1.1\pm 0.3$ km s-1 arcmin-1). This is
the first time that rotation for the stellar component is detected with high
significance for this galaxy.

We further analysed the probable members looking at the kinematic properties
of each individual field and for different stellar populations. The individual
fields analysis confirmed the presence of a rotation signal, and the decrease
with radius of the velocity dispersion values that appears to flatten moving
outward. We found a general agreement with the velocity field of the neutral
HI component, as analysed in Read et al. (2016), though, on average, the
stellar values showed higher velocity dispersion and a slower rotation, which
can be ascribed to asymmetric drift. 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. It appears that the kinematics of
this younger stellar component is still coupled to that of the gas.

The chemical analysis was conducted on the sub-sample of RGB stars exploiting
the wavelength region of the Ca II triplet lines. We obtained a median
${\rm[Fe/H]}=-1.06$ dex, with scatter values of $\sigma_{\rm MAD}=0.29$ dex
and $\sigma_{\rm intrinsic}=0.26$ dex, in fairly 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).

In future dedicated efforts, we plan 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. Finally, this MUSE data-set could also be used
to analyse the ionized gas medium using available emission lines, like
$H_{\alpha}$ or nebular oxygen lines in order to get an insight on the
kinematics of the shell-like gas structure visible in field F1 and compare it
with that of the youngest stars.

###### Acknowledgements.

We thank M.-R. Cioni, L. R. Patrick and J. Read for useful discussions and
comments. 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) 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.3, Pérez & Granger
2007), Numpy (v1.21, Harris et al. 2020), Scipy (v1.7, Virtanen et al. 2020),
Matplotlib (v3.5, Hunter 2007), Astropy (v5.0, Astropy Collaboration et al.
2022) and Scikit-learn (v1.0, Pedregosa et al. 2011) packages.

## References

  * Albert et al. (2000) Albert, L., Demers, S., & Kunkel, W. E. 2000, AJ, 119, 2780
  * Alonso et al. (1999) Alonso, A., Arribas, S., & Martínez-Roger, C. 1999, A&AS, 140, 261
  * Arroyo-Polonio et al. (2023) Arroyo-Polonio, J. M., Battaglia, G., Thomas, G. F., et al. 2023, A&A, 677, A95
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167
  * Bacchini et al. (2020) Bacchini, C., Fraternali, F., Iorio, G., et al. 2020, A&A, 641, A70
  * Bacon et al. (2014) Bacon, R., Vernet, J., Borisova, E., et al. 2014, The Messenger, 157, 13
  * Battaglia & Nipoti (2022) Battaglia, G. & Nipoti, C. 2022, Nature Astronomy, 6, 659
  * Battaglia et al. (2022) Battaglia, G., Taibi, S., Thomas, G. F., & Fritz, T. K. 2022, A&A, 657, A54
  * Belland et al. (2020) Belland, B., Kirby, E., Boylan-Kolchin, M., & Wheeler, C. 2020, ApJ, 903, 10
  * Bennet et al. (2023) Bennet, P., Patel, E., Sohn, S. T., et al. 2023, arXiv e-prints, arXiv:2312.09276
  * Berger et al. (2018) Berger, T. A., Kudritzki, R.-P., Urbaneja, M. A., et al. 2018, ApJ, 860, 130
  * Bernard et al. (2007) Bernard, E. J., Aparicio, A., Gallart, C., Padilla-Torres, C. P., & Panniello, M. 2007, AJ, 134, 1124
  * Bernard et al. (2010) Bernard, E. J., Monelli, M., Gallart, C., et al. 2010, ApJ, 712, 1259
  * Bodensteiner et al. (2020) Bodensteiner, J., Shenar, T., & Sana, H. 2020, A&A, 641, A42
  * Borissova et al. (2000) Borissova, J., Georgiev, L., Kurtev, R., et al. 2000, Rev. Mexicana Astron. Astrofis., 36, 151
  * Borissova et al. (2004) Borissova, J., Kurtev, R., Georgiev, L., & Rosado, M. 2004, A&A, 413, 889
  * Bresolin et al. (2007) Bresolin, F., Urbaneja, M. A., Gieren, W., Pietrzyński, G., & Kudritzki, R.-P. 2007, ApJ, 671, 2028
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
  * Buck et al. (2019) Buck, T., Macciò, A. V., Dutton, A. A., Obreja, A., & Frings, J. 2019, MNRAS, 483, 1314
  * Bundy et al. (2015) Bundy, K., Bershady, M. A., Law, D. R., et al. 2015, ApJ, 798, 7
  * Carrera et al. (2013) Carrera, R., Pancino, E., Gallart, C., & del Pino, A. 2013, MNRAS, 434, 1681
  * Cenarro et al. (2001) Cenarro, A. J., Cardiel, N., Gorgas, J., et al. 2001, MNRAS, 326, 959
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints, arXiv:1612.05560
  * Chiti et al. (2021) Chiti, A., Frebel, A., Simon, J. D., et al. 2021, Nature Astronomy, 5, 392
  * Chun et al. (2015) Chun, S.-H., Jung, M., Kang, M., Kim, J.-W., & Sohn, Y.-J. 2015, A&A, 578, A51
  * Chun et al. (2022) Chun, S.-H., Yoon, S.-C., Oh, H., Park, B.-G., & Hwang, N. 2022, ApJ, 939, 28
  * Cole et al. (1999) Cole, A. A., Tolstoy, E., Gallagher, John S., I., et al. 1999, AJ, 118, 1657
  * Collins & Read (2022) Collins, M. L. M. & Read, J. I. 2022, Nature Astronomy, 6, 647
  * Cortese et al. (2014) Cortese, L., Fogarty, L. M. R., Ho, I. T., et al. 2014, ApJ, 795, L37
  * de Mink et al. (2013) de Mink, S. E., Langer, N., Izzard, R. G., Sana, H., & de Koter, A. 2013, ApJ, 764, 166
  * Demers et al. (2006) Demers, S., Battinelli, P., & Kunkel, W. E. 2006, ApJ, 636, L85
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1
  * Gallart et al. (1996) Gallart, C., Aparicio, A., & Vilchez, J. M. 1996, AJ, 112, 1928
  * Gallart et al. (2015) Gallart, C., Monelli, M., Mayer, L., et al. 2015, ApJ, 811, L18
  * Garcia et al. (2014) Garcia, M., Herrero, A., Najarro, F., Lennon, D. J., & Alejandro Urbaneja, M. 2014, ApJ, 788, 64
  * Garcia et al. (2009) Garcia, M., Herrero, A., Vicente, B., et al. 2009, A&A, 502, 1015
  * Georgiev et al. (1999) Georgiev, L., Borissova, J., Rosado, M., et al. 1999, A&AS, 134, 21
  * Girardi et al. (2000) Girardi, L., Bressan, A., Bertelli, G., & Chiosi, C. 2000, A&AS, 141, 371
  * Gonneau et al. (2016) Gonneau, A., Lançon, A., Trager, S. C., et al. 2016, A&A, 589, A36
  * Gull et al. (2022) Gull, M., Weisz, D. R., Senchyna, P., et al. 2022, ApJ, 941, 206
  * Hamren et al. (2016) Hamren, K., Beaton, R. L., Guhathakurta, P., et al. 2016, ApJ, 828, 15
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357–362
  * Higgs et al. (2021) Higgs, C. R., McConnachie, A. W., Annau, N., et al. 2021, MNRAS, 503, 176
  * Holtzman et al. (2006) Holtzman, J. A., Afonso, C., & Dolphin, A. 2006, ApJS, 166, 534
  * Hunter et al. (2012) Hunter, D. A., Ficut-Vicas, D., Ashley, T., et al. 2012, AJ, 144, 134
  * Hunter (2007) Hunter, J. D. 2007, Computing in science & engineering, 9, 90
  * Husser (2012) Husser, T.-O. 2012, PhD thesis, Georg August University of Gottingen, Germany
  * Husser et al. (2013) Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6
  * Iorio et al. (2017) Iorio, G., Fraternali, F., Nipoti, C., et al. 2017, MNRAS, 466, 4159
  * Jensen et al. (2024) Jensen, J., Hayes, C. R., Sestito, F., et al. 2024, MNRAS, 527, 4209
  * Júlio et al. (2023) Júlio, M. P., Brinchmann, J., Zoutendijk, S. L., et al. 2023, A&A, 678, A38
  * Kacharov et al. (2017) Kacharov, N., Battaglia, G., Rejkuba, M., et al. 2017, MNRAS, 466, 2006
  * Kamann et al. (2013) Kamann, S., Wisotzki, L., & Roth, M. M. 2013, A&A, 549, A71
  * Kirby et al. (2014) Kirby, E. N., Bullock, J. S., Boylan-Kolchin, M., Kaplinghat, M., & Cohen, J. G. 2014, MNRAS, 439, 1015
  * Kirby et al. (2013) Kirby, E. N., Cohen, J. G., Guhathakurta, P., et al. 2013, ApJ, 779, 102
  * Koleva et al. (2009) Koleva, M., Prugniel, P., Bouchard, A., & Wu, Y. 2009, A&A, 501, 1269
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
  * Lake & Skillman (1989) Lake, G. & Skillman, E. D. 1989, AJ, 98, 1274
  * Leaman et al. (2017) Leaman, R., Mendel, J. T., Wisnioski, E., et al. 2017, MNRAS, 472, 1879
  * Leaman et al. (2013) Leaman, R., Venn, K. A., Brooks, A. M., et al. 2013, ApJ, 767, 131
  * Leaman et al. (2012) Leaman, R., Venn, K. A., Brooks, A. M., et al. 2012, ApJ, 750, 33
  * Lee et al. (2003) Lee, H., Grebel, E. K., & Hodge, P. W. 2003, A&A, 401, 141
  * Longeard et al. (2023) Longeard, N., Jablonka, P., Battaglia, G., et al. 2023, MNRAS, 525, 3086
  * Lozinskaya et al. (2003) Lozinskaya, T. A., Moiseev, A. V., & Podorvanyuk, N. Y. 2003, Astronomy Letters, 29, 77
  * McConnachie (2012) McConnachie, A. W. 2012, AJ, 144, 4
  * McConnachie et al. (2021) McConnachie, A. W., Higgs, C. R., Thomas, G. F., et al. 2021, MNRAS, 501, 2363
  * McLeod et al. (2020) McLeod, A. F., Kruijssen, J. M. D., Weisz, D. R., et al. 2020, ApJ, 891, 25
  * McQuinn et al. (2017) McQuinn, K. B. W., Boyer, M. L., Mitchell, M. B., et al. 2017, ApJ, 834, 78
  * Menzies et al. (2015) Menzies, J. W., Whitelock, P. A., & Feast, M. W. 2015, MNRAS, 452, 910
  * Minor (2013) Minor, Q. E. 2013, ApJ, 779, 116
  * Minor et al. (2010) Minor, Q. E., Martinez, G., Bullock, J., Kaplinghat, M., & Trainor, R. 2010, ApJ, 721, 1142
  * Moe & Di Stefano (2017) Moe, M. & Di Stefano, R. 2017, ApJS, 230, 15
  * Moiseev & Lozinskaya (2012) Moiseev, A. V. & Lozinskaya, T. A. 2012, MNRAS, 423, 1831
  * Offner et al. (2023) Offner, S. S. R., Moe, M., Kratter, K. M., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 275
  * Oh et al. (2015) Oh, S.-H., Hunter, D. A., Brinks, E., et al. 2015, AJ, 149, 180
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Pérez & Granger (2007) Pérez, F. & Granger, B. E. 2007, Computing in Science & Engineering, 9
  * Pokhrel et al. (2020) Pokhrel, N. R., Simpson, C. E., & Bagetakos, I. 2020, AJ, 160, 66
  * Pols et al. (1991) Pols, O. R., Cote, J., Waters, L. B. F. M., & Heise, J. 1991, A&A, 241, 419
  * Pont et al. (2004) Pont, F., Zinn, R., Gallart, C., Hardy, E., & Winnick, R. 2004, AJ, 127, 840
  * Porter & Rivinius (2003) Porter, J. M. & Rivinius, T. 2003, PASP, 115, 1153
  * Pucha et al. (2019) Pucha, R., Carlin, J. L., Willman, B., et al. 2019, ApJ, 880, 104
  * Puls et al. (2005) Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669
  * Putman et al. (2021) Putman, M. E., Zheng, Y., Price-Whelan, A. M., et al. 2021, ApJ, 913, 53
  * Read et al. (2016) Read, J. I., Iorio, G., Agertz, O., & Fraternali, F. 2016, MNRAS, 462, 3628
  * Ren et al. (2022) Ren, T., Jiang, B., Ren, Y., & Yang, M. 2022, Universe, 8, 465
  * Richard et al. (2019) Richard, J., Bacon, R., Blaizot, J., et al. 2019, arXiv e-prints, arXiv:1906.01657
  * Rivinius et al. (2013) Rivinius, T., Carciofi, A. C., & Martayan, C. 2013, A&A Rev., 21, 69
  * Robin et al. (2003) Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523
  * Roth et al. (2018) Roth, M. M., Sandin, C., Kamann, S., et al. 2018, A&A, 618, A3
  * Sánchez et al. (2012) Sánchez, S. F., Rosales-Ortega, F. F., Marino, R. A., et al. 2012, A&A, 546, A2
  * Sánchez-Blázquez et al. (2006) Sánchez-Blázquez, P., Peletier, R. F., Jiménez-Vicente, J., et al. 2006, MNRAS, 371, 703
  * Sandage (1971) Sandage, A. 1971, ApJ, 166, 13
  * Schlafly & Finkbeiner (2011) Schlafly, E. F. & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Schootemeijer et al. (2022) Schootemeijer, A., Lennon, D. J., Garcia, M., et al. 2022, A&A, 667, A100
  * Sestito et al. (2023a) Sestito, F., Roediger, J., Navarro, J. F., et al. 2023a, MNRAS, 523, 123
  * Sestito et al. (2023b) Sestito, F., Zaremba, D., Venn, K. A., et al. 2023b, MNRAS, 525, 2875
  * Shao & Li (2014) Shao, Y. & Li, X.-D. 2014, ApJ, 796, 37
  * Sibbons et al. (2015) Sibbons, L. F., Ryan, S. G., Irwin, M., & Napiwotzki, R. 2015, A&A, 573, A84
  * Silich et al. (2006) Silich, S., Lozinskaya, T., Moiseev, A., et al. 2006, A&A, 448, 123
  * Simon (2019) Simon, J. D. 2019, ARA&A, 57, 375
  * Skillman et al. (2014) Skillman, E. D., Hidalgo, S. L., Weisz, D. R., et al. 2014, ApJ, 786, 44
  * Skillman et al. (2003) Skillman, E. D., Tolstoy, E., Cole, A. A., et al. 2003, ApJ, 596, 253
  * Spencer et al. (2018) Spencer, M. E., Mateo, M., Olszewski, E. W., et al. 2018, AJ, 156, 257
  * Spencer et al. (2017) Spencer, M. E., Mateo, M., Walker, M. G., et al. 2017, AJ, 153, 254
  * Starkenburg et al. (2010) Starkenburg, E., Hill, V., Tolstoy, E., et al. 2010, A&A, 513, A34
  * Stetson (1987) Stetson, P. B. 1987, PASP, 99, 191
  * Stilp et al. (2013) Stilp, A. M., Dalcanton, J. J., Skillman, E., et al. 2013, ApJ, 773, 88
  * Stoop et al. (2024) Stoop, M., Derkink, A., Kaper, L., et al. 2024, A&A, 681, A21
  * Swan et al. (2016) Swan, J., Cole, A. A., Tolstoy, E., & Irwin, M. J. 2016, MNRAS, 456, 4315
  * Taibi et al. (2022) Taibi, S., Battaglia, G., Leaman, R., et al. 2022, A&A, 665, A92
  * Taibi et al. (2020) Taibi, S., Battaglia, G., Rejkuba, M., et al. 2020, A&A, 635, A152
  * Tautvaišienė et al. (2007) Tautvaišienė, G., Geisler, D., Wallerstein, G., et al. 2007, AJ, 134, 2318
  * Tolstoy et al. (2009) Tolstoy, E., Hill, V., & Tosi, M. 2009, ARA&A, 47, 371
  * Udalski et al. (2001) Udalski, A., Wyrzykowski, L., Pietrzynski, G., et al. 2001, Acta Astron., 51, 221
  * Valdes et al. (2004) Valdes, F., Gupta, R., Rose, J. A., Singh, H. P., & Bell, D. J. 2004, ApJS, 152, 251
  * van de Ven et al. (2006) van de Ven, G., van den Bosch, R. C. E., Verolme, E. K., & de Zeeuw, P. T. 2006, A&A, 445, 513
  * van den Bergh (2000) van den Bergh, S. 2000, The Galaxies of the Local Group
  * Van Rossum & Drake (2009) Van Rossum, G. & Drake, F. L. 2009, Python 3 Reference Manual (Scotts Valley, CA: CreateSpace)
  * Vaz et al. (2023) Vaz, D., Brinchmann, J., Zoutendijk, S. L., et al. 2023, A&A, 678, A59
  * Vazdekis et al. (2016) Vazdekis, A., Koleva, M., Ricciardelli, E., Röck, B., & Falcón-Barroso, J. 2016, MNRAS, 463, 3409
  * Vazdekis et al. (2012) Vazdekis, A., Ricciardelli, E., Cenarro, A. J., et al. 2012, MNRAS, 424, 157
  * Vinko et al. (1998) Vinko, J., Remage Evans, N., Kiss, L. L., & Szabados, L. 1998, MNRAS, 296, 824
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Walker et al. (2009) Walker, M. G., Mateo, M., Olszewski, E. W., Sen, B., & Woodroofe, M. 2009, AJ, 137, 3109
  * Weilbacher et al. (2012) Weilbacher, P. M., Streicher, O., Urrutia, T., et al. 2012, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8451, Design and capabilities of the MUSE data reduction software and pipeline, 84510B
  * Weldrake et al. (2003) Weldrake, D. T. F., de Blok, W. J. G., & Walter, F. 2003, MNRAS, 340, 12
  * Wheeler et al. (2017) Wheeler, C., Pace, A. B., Bullock, J. S., et al. 2017, MNRAS, 465, 2420
  * Wolf et al. (2010) Wolf, J., Martinez, G. D., Bullock, J. S., et al. 2010, MNRAS, 406, 1220
  * Wolf (1906) Wolf, M. 1906, MNRAS, 67, 91
  * Zoutendijk et al. (2020) Zoutendijk, S. L., Brinchmann, J., Boogaard, L. A., et al. 2020, A&A, 635, A107

##  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\arcmin\times
27\arcmin$) 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 toward 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 system171717Catalogues can be found at
http://astronomy.nmsu.edu/holtz/archival/ic1613/html/ic1613.html. 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\sim 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\sim 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_{\rm HST}-V_{\rm HST}=c_{I}\,(B_{\rm HST}-V_{\rm 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_{\rm Sub}$ magnitudes, whose agreement with the
$I_{\rm 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 $\sim 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_{\rm MUSE}-V_{\rm ref}=c_{V}\,(V_{\rm
ref}-I_{\rm ref})+ZP_{V}$, and $I_{\rm MUSE}-I_{\rm ref}=c_{I}\,(V_{\rm
ref}-I_{\rm 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 [$\alpha$/M] from which to start
the fit. We obtained the initial values from a set of isochrones (Girardi et
al. 2000; Bressan et al. 2012) with ${\rm[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_{\rm eff}$ and ${\rm log}(g)$.
For the spectral fitting, we used the available PHOENIX templates, convolved
with the MUSE line spread function, with solar scale chemical composition,
thus fixing ${\rm[\alpha/M]}=0$ dex. We let the code find the best fit of all
other parameters (i.e., $T_{\rm eff}$, ${\rm log}(g)$ and ${\rm[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_{\rm eff}$ (from 2500 K to 15000 K) at fixed ${\rm[Fe/H]}=-1.5$ dex.
Templates were shifted at several l.o.s. velocities, from $-500$ km s-1 to
$50$ km s-1 at step of 50 km s-1. As initial guess for the velocity fitting,
we set $-230$ km s-1 (i.e. approximately IC 1613 systemic velocity). 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\pm 0.56$ km s-1), although we noticed some failure around the extremes
(i.e. for $-500$ km s-1 and $50$ km s-1) 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 $\sim-500$ km s-1,
while the variation at $50$ km s-1 is small enough to keep such possible
contaminants away from the bulk velocity of IC 1613.

Regarding the other parameters, $T_{\rm eff}$ values were also well recovered,
except for some scatter of $\sim 2000$ K at temperatures $>10000$ K; ${\rm
log}(g)$ was recovered with a significant scatter ($>2$ dex) for the hottest
stars and those with $T_{\rm eff}\sim 6500$ K (the latter a minority in our
main sample); ${\rm[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 to 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 sub-samples. For
the RGS, we selected stars with $T_{\rm eff}\sim 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/NCaT around $[3.5,5,10,15,25]$. The same was done for the
MS stars, but for stars with $T_{\rm eff}\sim 10000$ K, generating mock
spectra at S/N550 around $[5,10,20,50]$. Results are reported in Table 6.

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-1. The effective temperature was also
well recovered, although the errors were as much as $\sim 50$ K lower than the
measurement scatter. On the other hand, ${\rm log}(g)$ and ${\rm[Fe/H]}$ were
relatively 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${}_{550}=5$ were
underestimated by more than $15$ km s-1, while at the highest S/N the
deviation was less than 1 km s-1. 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-1. 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 ${\rm 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 ${\rm[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 ${\rm[Fe/H]}$
values from the two methods of $\sim 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 6)
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 $\sim 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 ${\rm[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
${\rm[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_{\rm eff}<15000$ K (i.e. up to the grid limit of the PHOENIX
library) and S/N${}_{550}>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 ${\rm[\alpha/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 ${\rm[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 $\sim 0.1$ dex, Bresolin et al. 2007, e.g.).

Table 6: Results from the Monte Carlo tests.

S/N | $V_{\rm rad,0}$ | $T_{\rm eff,0}$ | ${\rm log}(g)_{0}$ | ${\rm[Fe/H]}_{0}$ | $\frac{(V_{\rm rad}-V_{\rm rad,0})}{\delta_{\rm V_{\rm rad}}}$ | $\frac{(T_{\rm eff}-T_{\rm eff,0})}{\delta_{\rm T_{\rm eff}}}$ | $\frac{({\rm log}(g)-{\rm log}(g)_{0})}{\delta_{\rm{\rm log}(g)}}$ | $\frac{({\rm[Fe/H]}-{\rm[Fe/H]}_{0})}{\delta_{\rm{\rm[Fe/H]}}}$
---|---|---|---|---|---|---|---|---
(pxl-1) | (km s-1) | (K) | (dex) | (dex) |  |  |  |
3.5 | $-211\pm 14$ | $5523\pm 457$ | $3.3\pm 1.2$ | $-0.2\pm 0.4$ | $0.0,1.0$ | $0.0,1.2$ | $0.1,1.2$ | $0.0,1.1$
5 | $-247\pm 9$ | $5065\pm 289$ | $2.7\pm 0.9$ | $-0.8\pm 0.3$ | $0.1,0.9$ | $0.0,1.1$ | $-0.1,1.2$ | $0.2,1.2$
10 | $-235\pm 5$ | $5499\pm 145$ | $2.3\pm 0.5$ | $-0.3\pm 0.2$ | $0.0,1.0$ | $0.3,1.1$ | $0.2,1.2$ | $0.2,1.1$
15 | $-223\pm 4$ | $4851\pm 109$ | $1.3\pm 0.4$ | $-0.7\pm 0.2$ | $0.0,0.9$ | $0.1,1.3$ | $-0.1,1.2$ | $0.2,1.3$
25 | $-234\pm 2$ | $4965\pm 64$ | $2.3\pm 0.3$ | $-0.6\pm 0.1$ | $0.1,1.0$ | $0.0,1.1$ | $0.0,1.3$ | $0.0,1.0$
5 | $-249\pm 29$ | $9090\pm 1077$ | $3.9\pm 0.9$ | $1.0\pm 2.2$ | $-0.2,1.5$ | $-0.7,1.2$ | $-0.3,1.1$ | $-0.7,0.5$
10 | $-223\pm 13$ | $7206\pm 187$ | $3.5\pm 0.5$ | $-0.8\pm 0.3$ | $0.0,1.1$ | $-0.3,1.2$ | $-0.3,1.8$ | $0.1,1.4$
20 | $-228\pm 7$ | $7581\pm 94$ | $3.0\pm 0.2$ | $-1.2\pm 0.2$ | $-0.1,1.1$ | $0.0,1.3$ | $-0.1,1.5$ | $0.1,1.5$
50 | $-224\pm 3$ | $11213\pm 117$ | $2.9\pm 0.1$ | $-0.6\pm 0.1$ | $0.0,1.2$ | $0.0,1.2$ | $0.1,1.3$ | $0.1,1.2$

181818From left to right, the first column indicates the S/N obtained from the
continuum around the CaT lines, for the first five rows, and around 5500 Å,
otherwise; the next fours columns are the l.o.s. velocity, effective
temperature, surface gravity and metallicity recovered by spexxy for the
selected stars; the last four columns are the mean and standard deviation of
the difference between the mock values and the measured ones, scaled for the
mock errors of the considered parameters.

###  B.2 Comparison between SPEXXY and ULySS

Figure 10: 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^{2}_{\rm rad,Spexxy}+\delta^{2}_{\rm rad,ULySS}}$. Right: direct
comparison between the effective temperature values obtained with the two
codes; the solid black line is for reference only.

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_{\rm 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),
with $T_{\rm 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 $\sim 75$% of the measured velocities are in
agreement within the combined uncertainties, obtained as
$\sqrt{\delta^{2}_{\rm rad,Spexxy}+\delta^{2}_{\rm rad,ULySS}}$. The
comparison is shown in the left panel of Fig. 10 where for clarity we reported
only those measurements with a combined error $<30$ km s-1. Nevertheless, the
difference between the ULySS and spexxy measurements showed an average
positive shift of $\sim 5$ km s-1. What is more, a higher shift of $\sim 10$
km s-1 is observed for those objects with a S/N${}_{\rm 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-1 for templates with effective
temperatures of $\sim 4500$ K and greater than $\sim 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 ($\lesssim 10$), where
spexxy showed its strength in the recovery of l.o.s. velocities.

The comparison of effective temperatures, as shown in the right-hand panel of
Fig. 10, gave good results for $T_{\rm eff}<7500$ K, although the spexxy
values tend to be systematically $\sim 500$ K higher than the ULySS values. At
higher $T_{\rm 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 optimized 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 $\sim 500$ K, for a
quality PCL-flag between 3 and 4.

##  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.

Table 7: Observing log of the VLT/MUSE observations of IC 1613 under ESO
programme 097.B-0373 (PI: G. Battaglia).

Field | Position (RA, Dec) | Date / Hour | Exp. | Airmass | DIMM Seeing
---|---|---|---|---|---
| (J2000) | (UT) | (s) |  | (arcsec)
F1 | 01:04:49.00, +02:07:16.3 | 2016-08-12 / 09:01 | 1385 | 1.14 | 0.79
| 01:04:49.03, +02:07:16.0 | 2016-08-12 / 09:26 | 1385 | 1.17 | 0.68
| 01:04:49.00, +02:07:15.6 | 2016-08-13 / 07:57 | 1385 | 1.12 | 0.65
| 01:04:48.97, +02:07:15.8 | 2016-08-13 / 08:23 | 1385 | 1.12 | 0.64
| 01:04:49.03, +02:07:15.6 | 2016-08-14 / 08:29 | 1385 | 1.12 | 0.67
| 01:04:49.00, +02:07:15.9 | 2016-08-14 / 08:54 | 1385 | 1.14 | 0.68
| 01:04:48.97, +02:07:16.3 | 2016-09-02 / 03:40 | 1385 | 1.76 | 0.60
| 01:04:48.98, +02:07:15.3 | 2016-09-02 / 04:05 | 1385 | 1.56 | 0.62
F2 | 01:04:40.69, +02:04:45.2 | 2016-09-28 / 04:00 | 1385 | 1.19 | 0.73
| 01:04:40.70, +02:04:46.1 | 2016-09-28 / 04:25 | 1385 | 1.15 | 0.54
| 01:04:40.67, +02:04:45.6 | 2016-09-29 / 05:40 | 1385 | 1.13 | 0.52
| 01:04:40.73, +02:04:45.9 | 2016-09-29 / 06:05 | 1385 | 1.15 | 0.56
| 01:04:40.73, +02:04:45.4 | 2016-09-30 / 03:56 | 1385 | 1.18 | 0.69
| 01:04:40.70, +02:04:45.7 | 2016-09-30 / 04:21 | 1385 | 1.14 | 0.71
| 01:04:40.72, +02:04:45.2 | 2016-09-30 / 04:52 | 1385 | 1.12 | 0.76
| 01:04:40.68, +02:04:45.9 | 2016-09-30 / 05:17 | 1385 | 1.12 | 0.73
F3 | 01:04:29.44, +02:03:35.3 | 2016-10-04 / 04:41 | 1385 | 1.12 | 0.67
| 01:04:29.39, +02:03:35.4 | 2016-10-04 / 05:06 | 1385 | 1.12 | 0.55
| 01:04:29.42, +02:03:35.5 | 2016-10-04 / 05:36 | 1385 | 1.14 | 0.70
| 01:04:29.42, +02:03:36.0 | 2016-10-04 / 06:01 | 1385 | 1.17 | 0.83
| 01:04:29.41, +02:03:35.1 | 2016-12-21 / 01:16 | 1385 | 1.21 | 0.79
| 01:04:29.39, +02:03:35.9 | 2016-12-21 / 01:41 | 1385 | 1.28 | 0.64
| 01:04:29.43, +02:03:35.7 | 2016-12-22 / 01:05 | 1385 | 1.20 | 0.92
| 01:04:29.44, +02:03:36.0 | 2016-12-22 / 01:30 | 1385 | 1.26 | 0.92
191919From left to right, column names indicate: the pointing field name; the
field center coordinates; observing date and starting time of the scientific
exposure; the exposure time in seconds; the starting airmass; the average DIMM
seeing during the exposure in arcsec.

Figure 11:  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, $V_{\rm rad}$,
$T_{\rm eff}$, 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.

Figure 12: Figure 11. Continued. From top to bottom examples of: an M giant
star; a C star; a Be star.

Figure 13: Figure 11. Continued. From top to bottom examples of: a background
galaxy; a blend with an high redshift emitting-galaxy; an example of ionized
gas in IC 1613.