© ESO 2020

1. Introduction

Dwarf galaxies are the least massive, yet the most dark-matter-dominated galactic systems observed (e.g., Mateo 1998; Battaglia et al. 2013; Walker 2013). In the Local Group (LG), the nearest ones to the largest spirals, meaning the Milky Way (MW) and M 31, are gas-poor dwarf spheroidal systems (dSphs) with no on-going star formation (Tolstoy et al. 2009). Although they share the same morphology, their full star formation histories show complex evolutionary pathways (Gallart et al. 2015).

Due to their small masses, the formation and evolution of these galaxies could be strongly influenced by environmental effects. The orbital properties of the dwarf galaxy satellites of the MW, obtained by integrating, back in time, their present-day systemic motion (see e.g., Gaia Collaboration 2018; Fritz et al. 2018; Simon 2018, for Gaia-DR2 based determinations), are consistent with repeated pericentric passages for several of these objects. In practice, (unknown) factors such as the triaxiality of the MW’s potential, or interactions between satellites, among others, introduce significant uncertainties in the reconstruction of their full orbital history (e.g., Lux et al. 2010), but a comparison with the properties of dark matter sub-halos around MW-sized hosts does suggest that most MW satellites fell in the MW halo at intermediate to early times (see Rocha et al. 2012; Wetzel et al. 2015). A large body of work has therefore focused on exploring the impact of tidal and/or ram-pressure stripping caused by a MW-sized host onto various properties of the dSph satellites of the MW (see e.g., Piatek & Pryor 1995; Mayer et al. 2001a,b, 2006; Read et al. 2006; Muñoz et al. 2008; Klimentowski et al. 2009; Kazantzidis et al. 2011; Pasetto et al. 2011; Battaglia et al. 2015; Iorio et al. 2019). However, as was recently shown by Hausammann et al. (2019), the ram-pressure stripping induced by a host halo has its limits in the actual quench and gas depletion of a dSph. Internal effects, like stellar feedback resulting from episodic star formation and supernova-driven winds, play an important role too (see e.g., Sawala et al. 2010; Bermejo-Climent et al. 2018; Revaz & Jablonka 2018). Recent hydro-dynamic cosmological simulations have indeed shown that both internal and environmental mechanisms are necessary to the reproduction of the observed properties of the LG dwarf galaxies (see e.g., Brooks & Zolotov 2014; Garrison-Kimmel et al. 2014; Wetzel et al. 2016; Sawala et al. 2016). Stellar feedback, for example, was also found to be an important ingredient in enhancing tidal-stirring effects (Kazantzidis et al. 2017).

Given the multitude of physical mechanisms affecting low-mass galaxy evolution, data on the ages, chemical abundances, spatial distribution, and kinematics of the stellar component of LG dwarf galaxies are needed to understand the observed diversity of these systems. Detailed observations of MW satellites have been accumulated over the past years including large spectroscopic datasets (e.g., Tolstoy et al. 2004; Battaglia et al. 2006, 2008, 2011; Walker et al. 2009a; Kirby et al. 2011; Lemasle et al. 2012, 2014; Hendricks et al. 2014; Spencer et al. 2017). The observations showed little to no sign of internal rotation (e.g., Battaglia et al. 2008; Walker et al. 2008; Leaman et al. 2013; Wheeler et al. 2017), pronounced radial metallicity gradients (e.g., Battaglia et al. 2006, 2011; Walker et al. 2009a; Kirby et al. 2011; Leaman et al. 2013) and multiple populations with different chemo-dynamical properties (e.g., Tolstoy et al. 2004; Battaglia et al. 2006, 2008). Nonetheless, very few simulations have explored the link between the formation scenarios and the internal kinematic and chemical status of dwarf galaxies (Revaz et al. 2009; Schroyen et al. 2013; Revaz & Jablonka 2018).

It is likely that the MW influenced at least in part the evolution of its satellites. Therefore, isolated dSphs represent a tool that is valuable in the understanding of the intrinsic properties of the systems that have spent all or most of their life in a more benign environment. They are crucial to our understanding of the possible formation scenarios. In the LG, just a handful of dSphs are found in isolation: namely And XVIII, Cetus, and Tucana. Although it is probable that in the past they may have interacted just once with the MW or M 31 based on their line-of-sight (l.o.s.) systemic velocity (Lewis et al. 2007; Fraternali et al. 2009; but see also Sales et al. 2010; Teyssier et al. 2012), the fact that they spent most of their lives in isolation makes them ideal targets to contrast with the satellite dwarfs.

This work is part of a larger body of studies aiming to improve our knowledge of the observed properties of a selected sample of isolated LG dwarf galaxies – Phoenix, (Kacharov et al. 2017); Cetus, (Taibi et al. 2018); Aquarius, (Hermosa Muñoz et al. 2020) – mainly exploiting VLT/FORS2 multi-object spectroscopic data. Here we focus on the Tucana dSph.

Tucana is an early-type dwarf galaxy found in extreme isolation with a heliocentric distance of D_⊙  =  887  ±  49 kpc (Bernard et al. 2009), which places it at more than 1 Mpc away from M 31 and the LG center, with only the Phoenix dwarf found within ∼500 kpc of it (McConnachie 2012). Photometric observations have shown that the galaxy is mainly old and metal-poor, with an extended horizontal branch (HB) and a population of variable stars (see e.g., Saviane et al. 1996; Bernard et al. 2009). The structural analysis by Saviane et al. (1996) showed a highly flattened system (e ∼ 0.5) with a surface density profile well-described by an exponential fit. The recovery of the full star formation history (SFH) from deep HST/ACS observations reaching the oldest main sequence turn-off by Monelli et al. (2010) showed that Tucana formed the majority of its stars more than 9 Gyr ago. It experienced a strong initial period of star formation (SF) starting very early on (∼13 Gyr ago). Tucana harbors at least two stellar sub-populations based on observed splitting of the HB, double red giant branch (RGB) bump, and the luminosity-period properties of the RR-Lyrae, which imply that this system experienced at least two early phases of SF in a short period of time. Using the same HST/ACS dataset, Savino et al. (2019) refined the HB analysis, showing that Tucana experienced two initial episodes of sustained SF followed by a third less intense, but more prolonged one, ending between 6 and 8 Gyr ago. The spatial analysis of the same dataset indicates the presence of a population age gradient inside ∼4 R_e (Monelli et al. 2010; Hidalgo et al. 2013; Savino et al. 2019).

The first spectroscopic study of individual stars in the Tucana dSph was conducted by Fraternali et al. (2009), with the VLT/FORS2 instrument obtaining a relatively small sample of ∼20 RGB probable member stars. They reported the systemic velocity and velocity dispersion values (${\overline{\mathit{v}}}_{\mathrm{sys}}=194.0\pm 4.3$ km s⁻¹ and ${\sigma }_{\mathrm{v}}=15.8_{-3.1}^{+4.1}$ km s⁻¹) for the galaxy, together with the presence of a maximum rotation signal of ∼16 km s⁻¹. They also determined a mean metallicity value of [Fe/H] = − 1.95 ± 0.15 dex with a dispersion of 0.32 ± 0.06 dex. The determination of the systemic velocity ruled out an association with a nearby HI cloud, confirming that Tucana is devoid of neutral gas (down to a HI mass of 1.5 × 10⁴ M_⊙), and in addition it is moving away from the LG barycenter (Fraternali et al. 2009). If bound, Tucana has not reached its apocenter yet. It is indeed possible that a past interaction between Tucana and the MW happened around 10 Gyr ago (roughly coinciding with the major drop in its SF; Sales et al. 2007; Fraternali et al. 2009; Teyssier et al. 2012).

A new spectroscopic study of Tucana’s RGB stars has recently been conducted by Gregory et al. (2019) using VLT/FLAMES data. Their sample of probable members is slightly larger than that of Fraternali et al. (2009), but covers a much more extended spatial area (up to ∼10 R_e  ∼  2 R_tidal). Their velocity dispersion value (${\sigma }_{\mathrm{v}}=14.4_{-2.3}^{+2.8}$ km s⁻¹) obtained for 36 probable members is similar to that of Fraternali et al. (2009), although their systemic velocity shows a significant offset (${\overline{\mathit{v}}}_{\mathrm{sys}}=216.7_{-2.8}^{+2.9}$ km s⁻¹); they also detect a velocity gradient of $k=7.6_{-4.3}^{+4.2}$ km s⁻¹ arcmin⁻¹ along the optical major axis. Performing a dynamical modeling of Tucana’s kinematic properties based on the FLAMES/GIRAFFE l.o.s. velocities of the probable members, the authors found a massive dark matter halo with a high central density. The implied dark matter halo mass profile is much denser than the other dwarfs of the LG, making Tucana the first exception of the too-big-to-fail problem (see e.g., Boylan-Kolchin et al. 2012). In fact, the pure N-body simulations in the framework of Λ cold dark matter (Λ-CDM) predict that the dwarf galaxies of the LG should live in denser halos than those inferred from observations. The fact that Tucana is found to reside in such a massive halo in agreement with Λ-CDM predictions seems to indicate that, during its evolution, it has been able to maintain its initial content of dark matter, independently of the internal and environmental mechanisms that have driven its evolution.

Motivated by the unique internal kinematics of Tucana and the importance of isolated dwarfs in disentangling evolutionary processes, in this study we present results from a new investigation of the kinematic and chemical properties of the stellar component of the Tucana dSph. We analyzed a new dataset of multi-object spectroscopic observations of 50 individual RGB stars taken with the VLT/FORS2 instrument targeting the near-IR wavelength region of the Ca II triplet (CaT) lines. To understand how systematics may influence the key results of our study, we further re-reduced the original datasets presented in Fraternali et al. (2009) and Gregory et al. (2019) and performed a combined analysis together with our own data. In this work, we present an in-depth and homogeneous analysis of all currently available spectroscopic data for the Tucana dSph.

The article is structured as follows: in Sect. 2, we present the data acquisition and reduction processes for all the datasets analyzed in this work. Section 3 is dedicated to the determination of the l.o.s. velocity and metallicity measurements. In Sect. 4, we describe the criteria applied to select likely member stars in the different datasets we present here. Section 5 shows the results from the kinematic analysis, presenting the determination of the galaxy systemic velocity and velocity dispersion, along with the search for a possible rotation signal and the implication for the dark matter halo properties of Tucana. In Sect. 6, we describe the determination of metallicities ([Fe/H]) and the subsequent chemical analysis. Finally, Sect. 7 is dedicated to our summary and conclusions, while in the appendices we report the detailed comparison between the measurements obtained for our dataset with those reported in the literature, along with supplementary material from the kinematic analysis. The parameters adopted for the Tucana dSph are summarized in Table 1.

2. Data acquisition and reduction processes

The analysis we present in this work required combining data obtained with different instruments and observational set-ups. Here, we provide details about the data acquisition and reduction processes we followed.

2.1. The P91 FORS2 dataset

The primary dataset analyzed in this work was obtained with the FORS2 instrument mounted at the UT1 (Antu) of the Very Large Telescope (VLT) at the ESO Paranal observatory. Observations were taken in service mode over several nights between July 2013 and July 2014 as part of the ESO program 091.B-0251, PI: M. Zoccali (see Table 2). The instrument was used in multi-object spectroscopic mode (MXU), which allows the observer to employ exchangeable masks with custom-cut slits. We used pre-imaging FORS2 photometry taken in Johnson V- and I-band, to allocate slits to stars with colors and magnitudes compatible with Tucana’s RGB. Slits that would otherwise have remained empty (five of them) were allocated to random targets in the same magnitude range. We selected 50 objects distributed over two overlapping masks of 27 slits each; the observation of four objects was repeated on purpose for internal accuracy measurements. Therefore, the selected targets covered an area up to the nominal King tidal radius of Tucana (R_tidal = 3.7′), as can be seen in Fig. 1.

Fig. 1. Spatial distribution (left) and color-magnitude diagram (right) of stars along the line of sight to the Tucana dSph. Black points represent the objects classified as stars in the FORS2 photometric data (see main text); red and blue dots indicate the P91-FORS2 MXU targets classified as probable members (i.e., with P >  0.05) and non-members, respectively. Yellow triangles and green squares represent the probable member stars from the P69-FORS2 and FLAMES datasets, respectively, which were added to the P91-FORS2 and analyzed through the text. The two observed FORS2 pointings are represented as large squares, while the ellipses denote the galaxy half-light radius and the spatial extension of the dataset (i.e., up to 5 × Re ∼ Rtidal). We note that the photometric data are not corrected for reddening.

The adopted instrumental set-up and observing strategy were the same as in our previous studies (Kacharov et al. 2017; Taibi et al. 2018; Hermosa Muñoz et al. 2020), so we only report the essentials here. We used the 1028z+29 holographic grism together with the OG590+32 order separation filter in order to cover the wavelength range between 7700 and 9500 Å. Slits had spatial sizes of 1″ × 10″ (8″ in some cases to avoid overlaps) in the first mask (Tuc0) and of 1″ × 8″ (7″ for overlaps) in the second one (Tuc1). This led to a binned spectral dispersion of 0.84 Å pxl⁻¹ and a resolving power of R = λ_cen/Δλ ∼ 2600 at λ_cen = 8600 Å (equivalent to a velocity resolution of 28 km s⁻¹ pxl⁻¹). Ten identical observing blocks (OBs¹) were taken for each pointing in order to reach the necessary signal-to-noise ratio (S/N) for velocity and metallicity measurements. In Table 2, we report a complete observing log.

The data were provided by ESO as individual OB-datasets within the FORS2 standard delivery plan. We adopted the same data-reduction process as in Taibi et al. (2018, hereafter, T18), based on IRAF² routines and custom-made python scripts. Briefly, our pipeline was developed to organize and reduce each OB-dataset independently. After making bias and flat-field corrections on the two-dimensional (2D) multi-object scientific and lamp-calibration frames, these were also cleaned from cosmic rays and bad rows. The 2D images are corrected for distortions in the spatial direction by rectifying their slit traces in order to cut them into individual 2D spectra. The arc-lamp spectra are then used to find the wavelength solution to calibrate the scientific exposures, with a typical RMS accuracy of 0.05 Å. The wavelength calibration also has the important effect of rectifying the sky lines, which had initially been curved by the instrument disperser, helping to reduce the residuals during the sky-subtraction part. The rectified wavelength-calibrated 2D individual scientific exposures were finally background subtracted, optimally extracted into 1D spectra and normalized by fitting the stellar continuum. The median S/N around the CaT for the individual exposures was ∼8 Å⁻¹.

We then followed the approach already presented in T18 and Hermosa Muñoz et al. (2020) to stack the repeated individual exposures together for each target. To do so, we needed to account for possible zero-point displacements in the wavelength calibration, small slit-centering shifts, and the different dates of observation. We used the numerous OH emission lines in the extracted sky background to refine the wavelength calibration of the individual spectra, using the IRAF fxcor task to cross-correlate with a reference sky spectrum over the region 8250 − 9000 Å. The calculated off-sets roughly varied between 5 and 25 km s⁻¹ with an average error of 2 km s⁻¹. The correction for slit-centering shifts was done using the through-slit frames typically taken before each scientific exposure. The offsets were calculated as the difference in pixels between the slit center and the star centroid for every target per mask. The correction for each target was taken as the median value of all its slit shifts; the associated error was the scaled median absolute deviation (MAD) of those values. We found slit shifts in the range ±0.1 − 9.5 km s⁻¹ with errors of ±1 − 5 km s⁻¹. Finally, we obtained the heliocentric correction for all the individual exposures using the IRAF rvcorrect task. All the above shifts were applied to the individual spectra using the IRAF dopcor task. We then regrouped and averaged the repeated exposures of each target, weighting them with their associated σ-spectra, given by the extraction procedure. Final error spectra were obtained accordingly. The median S/N around the CaT for the stacked spectra is ∼27 Å⁻¹ (see also Table C.1 for the properties of each individual observed star).

2.2. A new reduction for the Fraternali et al. (2009) FORS2 dataset

Fraternali et al. (2009, hereinafter, F09) presented l.o.s. velocities and metallicities for 23 individual stars with magnitudes and colors compatible with Tucana’s RGB (of which 17 were classified as members) from an earlier FORS2 MXU dataset (ESO program 69.B-0305(B), PI: E. Tolstoy). We wished to use this catalog to increase our sample size; however, a comparison with the l.o.s. velocities and metallicities between the targets in common (3) with the F09 catalog showed significant systematic shifts in both quantities (see Appendix A.3 for details). For the sake of homogeneity, we therefore reduced the F09 dataset again following the same procedure as adopted for the P91 sample. Hereinafter, we refer to this additional sample as the P69 FORS2 dataset.

The instrumental set-up of the P69 program was similar to the one we adopted for our observations, with the difference being that the size of each slit was of 1.2″ × 8″, which translates into a slightly lower spectral resolution (≈30 km s⁻¹ pxl⁻¹). A total of 45 initial targets were assigned to an equivalent number of slits packed into a single pointing centered on Tucana. The total exposure time of the observations was 5.2 h.

The reduction of some slits, in particular the background subtraction step, turned out to be particularly problematic, since these targets were not well-centered along the slit spatial direction, but placed at their edges. This led to 15 extracted spectra with high sky residuals, which made them unreliable for l.o.s. velocity and metallicity estimates. In addition, these targets had magnitudes and colors outside the RGB of Tucana. Therefore, we excluded them from the sample, along with two further objects whose extracted spectra did not show any CaT lines.

The final P69 sample was reduced to 28 objects, with five more than the initially published F09 catalog. Once stacked together, reliable spectra had an S/N  ∼  22 Å⁻¹ (see also Table C.2 for the properties of each individual observed star).

Another issue was the lack of through-slit images in the provided data, although in F09 it is reported that the objects were well-centered in the spectral direction after visual inspection during the spectroscopic run. However, we took into account the error related to the slit centering by assuming that it is a tenth of a pixel (∼3 km s⁻¹, which is the typical error found checking for the slit centering) and adding it in quadrature during the velocity estimation step.

As can be seen in Appendix A.2, there is good agreement for the measurements of the stars in common between this dataset (the P69 FORS2) and the P91 FORS2 (and for those cases where there is disagreement, the source of it can be traced back).

2.3. The FLAMES dataset

Gregory et al. (2019, hereinafter, G19) recently presented an additional sample of l.o.s. velocities (not metallicities) for individual stars in the direction of Tucana, taken with the FLAMES/GIRAFFE instrument at the VLT, as part of the ESO program 095.B-0133(A), PI: M. Collins. The spectrograph was used in MEDUSA mode, meaning in multi-fiber configuration, which allows for the simultaneous observation of up to 132 separate targets (sky fibers included). The instrument field of view (FoV) has a 25 arcmin diameter, and each fiber has an aperture on the sky of 1.2 arcsec. The grating used was the LR8, centered on 8817 Å and covering the CaT wavelength region, yielding a spectral resolution of R ∼ 6500.

The authors reported the detection of 36 probable member stars, out to very large distances from Tucana’s center, approximately up to 10 half-light radii. Given the larger spatial region probed by these data with respect to the FORS2 P91 and P69 datasets (compare Fig. 1 of this work to Fig. 1 in G19), it was interesting to explore whether the G19 catalog of l.o.s. velocities could be used together with our determinations from the P91 and P69 FORS2 data. The comparison of the six stars in common between the P91 and G19 samples yields a discrepancy of ∼30 km s⁻¹ for five of the six stars and of about −150 km s⁻¹ for the other object, which does not allow the direct combination of the measurements. The G19 paper also reports an offset of ∼23 km s⁻¹ between their velocities and those in the F09 catalog. We note that the offset of ∼30 km s⁻¹ is compatible with the 23 km s⁻¹ offset between G19 and F09 and the ∼7 km s⁻¹ offset we found when comparing our P91 velocities to the F09 catalog (see Appendix A.3 for further details). Given the above information, we proceeded to perform our own reduction of the FLAMES/GIRAFFE data, the characteristics of which we briefly describe below.

The observations were taken on six nights spread between June and September 2015, using two different fiber setups covering the same area: the first one (Tuc-1) had a total of 7 h of exposure time accumulated over seven OBs, and the second setup (Tuc-2) got 6 h of exposure time taken within six OBs, of which one was repeated twice. Each OB consisted of 3 × 1200 s exposures. The total number of individual targets was 164. The first setup had 14 fibers, and the second had 15 fibers allocated to empty sky regions distributed over the entire FoV. Only very few targets are spatially found inside the tidal radius of Tucana (∼30 targets), with the others scattered over an area much larger than the nominal extension of this galaxy. In fact, the two pointings are off-centered by ∼9′ from the optical center of Tucana.

The FLAMES data were downloaded from the science portal of the ESO archive as already processed spectra, meaning pre-reduced, wavelength-calibrated, extracted, and corrected to the barycentric velocity, but with no sky subtraction applied. For each OB, ESO delivers spectra stacked at the OB-level for the science targets and several auxiliary data, including the individual scientific exposures within an OB for both scientific targets and sky fibers³. The only remaining steps were to perform the sky subtraction, combine the spectra with repeated exposures, and, finally, calculate the radial velocities.

We first verified the quality of the wavelength calibration by cross-correlating the sky lines of the scientific exposures with a template sky spectrum, which was placed at the rest frame and obtained with a similar observational set-up (Battaglia et al. 2011)⁴. However, the delivered spectra are shifted to the heliocentric frame, including the sky lines, thus we had to remove this correction first. We then performed the cross-correlation for each spectrum of each OB dataset using the IRAF fxcor task, obtaining median offsets around 0.6 km s⁻¹ with a global scatter of 0.7 km s⁻¹. Therefore, the uncertainty related to the wavelength calibration resulted well below those from the velocity measurement (as shown later).

For the sky subtraction, we used the ESO skycorr tool (Noll et al. 2014). The idea behind this code is to adopt a physically motivated group scaling of the sky emission lines with respect to a reference sky spectrum according to their expected variability and given the date of the observations. We created the reference sky spectrum by first median combining the spectra of the fibers allocated to sky within the individual sub-exposures of each OB, and then by median combining the results for the individual sub-exposures. The optimized line groups in the reference sky spectrum are scaled to fit the emission lines in the science spectra and finally subtracted together with the sky continuum. Error spectra are also an output of the code. We used default input parameter while running skycorr. This method yielded satisfactory results for the majority of stars, particularly for spectra with low S/N.

The sky-subtracted science spectra were then normalized using a Chebyshev polynomial of order 3, together with their error spectra. Finally, repeated exposure of individual targets (including those the Tuc-1 and Tuc-2 set-ups have in common) were stacked together using a weighted average, with their error spectra combined accordingly. The typical S/N was ∼11 Å⁻¹, although in some cases it was as low as ∼2 Å⁻¹ (see also Table C.3 for the properties of each individual observed star).

The average S/N measured in our reduction of the FLAMES/GIRAFFE spectra is in good agreement with the S/N obtained from the GIRAFFE Exposure Time Calculator (ETC) considering typical values from the observed dataset: we used a black body template of T_eff = 4500 K, an I-band magnitude of 21, an airmass of 1.2, a moon illumination fraction of 0.2, a seeing of 1.0″, an object-fiber displacement of 0.3″ and a total exposure time of 10 h (36 000 s). With this setting we obtained a calculated S/N of 15 Å⁻¹, close to our typical S/N value.

2.4. Photometric data

Photometric data were used for all of the three datasets (FORS2 P91 and P69, and FLAMES/GIRAFFE) to exclude those objects whose magnitude and color were not compatible with being stars on Tucana’s RGB (see Sect. 4), and for the FORS2 datasets to determine metallicities from the equivalent width of the CaT lines using calibrations from the literature (see Sect. 3).

We used the pre-imaging FORS2 photometric catalog introduced in Sect. 2.1 to associate V- and I-band magnitudes to the target stars in the P91 and P69 spectroscopic datasets. The photometric catalog was astrometrized and the instrumental magnitudes were calibrated with the publicly available catalog from Holtzman et al. (2006) obtained with the Wide Field and Planetary Camera 2 of the Hubble Space Telescope (HST/WFPC2) being used as a reference. We used the aperture-photometry catalog provided in Johnson’s UBVRI-system to actually find the astrometric solution of the FORS2 catalog and to perform the photometric calibration using the suite of codes CataXcorr and CataComb, kindly provided to us by P. Montegriffo and M. Bellazzini (INAF-OAS).

The case of the FLAMES dataset, on the other hand, was different. Since we did not have a photometric catalog covering an area as wide as that of the spectroscopic targets, we instead used the publicly available photometry of individual point sources from the first data release of the Dark Energy Survey (DES-DR1, Abbott et al. 2018). We found a match for 154 out of 164 spectroscopic targets, considering a tolerance radius of 1 arcsec. To be conservative, we did not exclude the targets that did not have a match in the DES-DR1 photometry from the photometric selection. The DES-DR1 griz-photometry needed then to be converted first to the SDSS griz-system⁵ and finally to Johnson’s system (see Jordi et al. 2005). The FLAMES data showed magnitudes as high as V ∼ 23, which was also the limit of the DES-DR1 catalog, and, consequently, some targets had relatively large associated magnitude errors (δ_mag ∼ 0.2).

3. Line-of-sight velocity and metallicity measurements

The determination of l.o.s. velocities and metallicities ([Fe/H]) from the stacked spectra was done in the same way as in T18. While we were able to measure both l.o.s. velocities and metallicity for the FORS2 P91 and P69 data, only l.o.s. velocities could be determined for FLAMES/GIRAFFE data, due to the low S/N of those spectra.

Line-of-sight velocities were obtained using the fxcor task by cross-correlating with a synthetic spectrum resembling a low-metallicity RGB star convolved at the same spectral resolution of the dataset under consideration. For FORS2 we used the same template as in T18, and cross-correlated in the wavelength range 8400 − 8700 Å. For the FLAMES dataset we used a template from Zoccali et al. (2014), that is a synthetic spectrum of a star with T_eff = 4750 K, log(g) = 2.5 and [Fe/H] = −1.3 dex. In this case, we used the region around the two reddest lines of the CaT for the cross-correlation, since the first line often suffered from high residuals left by the subtraction of a sky emission line. This was not necessary for the FORS2 spectra, which suffer less from this problem due to their higher S/N.

In the following, we only keep objects whose stacked FORS2 and FLAMES spectra have a S/N ≳ 10 Å⁻¹, since this is the limit where the velocity errors provided by fxcor task appear reliable (see tests in Appendix A.1). Average velocity errors for the FORS2 (FLAMES) spectra were found to be ∼7 km s⁻¹ (∼9 km s⁻¹).

As can be seen from the histogram in the right panel Fig. 2, the three datasets show a clear peak in the l.o.s. velocities around 180 km s⁻¹, where we roughly expect to find the members of Tucana. The FLAMES dataset presents a higher fraction of contaminants due to the large area covered by the observations. With our homogeneous analysis, we show a clear detection of the stars with velocities compatible with Tucana in all three datasets and no significant offsets, unlike what we found in direct comparison with our velocity measurements and those published in F09 and G19, as shown in the left panel of Fig. 2. We refer the reader to Appendix A.3 and Sect. 5.2 for a detailed comparison with the studies taken from literature.

Fig. 2. Histogram of l.o.s. velocity measurements from the analyzed datasets. Left: comparing the velocities from the P91-FORS2 dataset with those from literature, i.e. Fraternali et al. (2009) and Gregory et al. (2019). Right: same comparison using instead the velocity measurements from our reduction of the P69-FORS2 and FLAMES datasets. Note that, in the left panel, the peaks of the histograms fall at different velocities, while in the right panel, where the datasets shown were analyzed homogeneously, these differences are absent.

To estimate the [Fe/H] values, we adopted the Starkenburg et al. (2010) relation, which is a function of the equivalent widths (EWs) of the two reddest CaT lines, and of the (V − V_HB) term, where V_HB is the mean magnitude of the galaxy’s horizontal branch (HB). We obtained the EWs from the continuum-normalized stacked spectra by fitting a Voigt profile over a window of 15 Å around the CaT lines of interest, and using the corresponding error-spectra as the flux uncertainty at each pixel in the fitting process. The errors on the EWs were then calculated from the covariance matrix of the fitting parameters. For the V_HB, we adopted the value of 25.32 from Bernard et al. (2009). Uncertainties on the [Fe/H] values were obtained by propagating the errors on the EWs, accordingly. Typical [Fe/H] errors were found to be around 0.15 − 0.25 dex.

4. Membership and kinematic analysis

Before proceeding with the analysis of Tucana’s kinematic and chemical properties, we needed to identify the stars that are probable members of Tucana and weed out possible contaminants (foreground MW stars and background galaxies). We followed the same steps for the membership selection in all the catalogs we analyzed.

We first selected targets located approximately along the RGB of Tucana, making our selection by magnitude and color. We used a set of isochrones (Girardi et al. 2000; Bressan et al. 2012) with age t_age  =  12.6 Gyr and [Fe/H]   ∼   −2.3 dex, and age t_age  =  8 Gyr and [Fe/H]   ∼   −0.4 dex to fix the blue and red color limits on the color-magnitude diagram (CMD). This color range was chosen to broadly cover the expected range of metallicities and stellar ages obtained from the SFH analysis of Tucana (Monelli et al. 2010; Savino et al. 2019). For the FLAMES dataset, due to the larger errors in the associated photometric dataset, we broadened the blue and red color limits applied to the FORS2 case by 0.2 mag. We further excluded the targets with spectra that showed high sky-residuals (and with S/N values < 10 Å⁻¹, as previously discussed) from all catalogs. Our P91 sample targets reduced from 50 to 43, the P69 objects reduced from 28 to 23, and the FLAMES one went from 164 to 58.

Given that the target selection of the three datasets was carried out completely independently, and might therefore include different biases, we proceeded by determining memberships and kinematic parameters by first considering our P91 FORS2 dataset on its own, a second time by combining it with the P69 sample, and, finally, while considering all the three sets. When we found common targets among the combined catalogs, we first kept those from the FORS2 datasets, with a higher priority for those of P91. Therefore, the combined FORS2 dataset has 63 targets, while including the FLAMES data gave a total of 109.

We then continued by assigning a membership probability to the individual targets in each considered dataset, applying a method based on the expectation maximization technique outlined in Walker et al. (2009b), but with the few modifications introduced by Cicuéndez et al. (2018). Briefly speaking, this approach makes it possible to carry out a Bayesian analysis to obtain the kinematic parameters of interest while assigning a probability of membership P_Mi to each i-star by maximizing the following log-likelihood equation:

$$\begin{array}{c}\hfill \ln L={\mathrm{\Sigma }}_i{P}_{{\mathrm{M}}_i}\ln \left[P_{\mathrm{mem}}{P}_{\mathrm{rad}}\right]+{\mathrm{\Sigma }}_i(1-P_{{\mathrm{M}}_i})\ln \left[P_{\mathrm{non}}(1-P_{\mathrm{rad}})\right],\end{array}$$(1)

where P_mem is the target’s probability distribution depending on its l.o.s. velocity, P_rad and P_non are the prior probabilities related to the surface density profile of the galaxy and the presence of possible contaminants, respectively, while P_Mi is defined as:

$$\begin{array}{c}\hfill P_{{\mathrm{M}}_i}=\frac{P_{\mathrm{mem}}{P}_{\mathrm{rad}}}{P_{\mathrm{mem}}{P}_{\mathrm{rad}}+P_{\mathrm{non}}(1-P_{\mathrm{rad}})}·\end{array}$$(2)

Both equations were adapted from Eqs. (3) and (4) of Walker et al. (2009b), respectively.

We run the Bayesian analysis using the MultiNest code (Feroz et al. 2009; Buchner et al. 2014), a multimodal nested sampling algorithm, in order to obtain the kinematic parameters and the membership probabilities all at once similarly to Cicuéndez et al. (2018). A further output of this code is the Bayesian evidence, which gives us the possibility to compare different kinematic models according to their statistical significance.

The spatial prior probability as a function of radius P_rad(R_i) accounts for the fact that we are more likely to observe a member star near the galaxy’s center than in its outer regions. To this aim, we assumed an exponentially decreasing surface-number density profile, which takes into account a uniform background surface-number density. The parameters of the profile were obtained from a VLT/VIMOS photometric dataset centered on Tucana, which was kindly provided by G. Beccari (ESO) and M. Bellazzini (INAF-OAS). This photometry was preferred to the FORS2 pre-imaging and the DES-DR1 photometry used for the CMD-selection since it is much deeper (almost 3 magnitudes in I-band), although it extends up to the tidal radius. To obtain the best-fitting structural parameters we applied a Bayesian Markov chain Monte Carlo (MCMC) analysis following the density profile of the RGB stars, as was done in Cicuéndez et al. (2018), while accounting for contamination. The assumed profile proved to be a good representation of the observed surface-number density profile for Tucana, and the best-fitting structural parameters (see Table B.1) are perfectly compatible with the parameters reported by Saviane et al. (1996)⁶.

The prior probability of contamination by foreground stars P_non was based on the Besançon model (Robin et al. 2003). The generated distribution of l.o.s. velocities was well-fit by a Gaussian profile (${\overline{\mathit{v}}}_{\mathrm{Bes}}=57$ km s⁻¹; σ_Bes = 99 km s⁻¹). The contamination model was generated in the direction of Tucana over an area equivalent to a FLAMES/GIRAFFE pointing selecting stars over the range of colors and magnitudes described above.

The l.o.s. velocity distribution of the probable member stars P_mem(v_i) was assumed to be Gaussian, as in T18, and accounted for the different kinematic models according to the following rotational term: v_rot(R_i) cos(θ − θ_i), with R_i being the angular distance from the galaxy’s center, θ_i the position angle (measured from north to east) of the ith target star, θ the position angle of the kinematic major axis (i.e., the direction of the velocity gradient, perpendicular to the axis of rotation), while v_rot(R_i) is the modeled rotational velocity term. We fit and compared three kinematic models: a dispersion-only model (i.e., with the velocity rotation term set to zero), a model with its rotational velocity linearly increasing with its radius (v_rot(R_i) = k R_i) and a flat one (v_rot(R_i) = v_c = constant).

The free kinematic parameters were as follows: the systemic velocity ${\overline{\mathit{v}}}_{\mathrm{hel}}$ and velocity dispersion σ_v common to the three models, the position angle θ of the kinematic major axis, the velocity gradient k of the linear rotation model, and the constant rotational velocity v_c of the flat model. In our definition, the position angle θ varies between 0° and 180°, which means that a rotation signal (either expressed as k or v_c) that has a negative sign implies a receding velocity on the west side of the galaxy (and would be equivalent to a positive gradient adding 180° to θ). The model evidences, Z, were combined through the Bayes factor, B_1, 2 = Z₁/Z₂, where the subscripts indicate a given model. To quantify the statistical significance of one model with respect to another, we made use of the Jeffreys scale, based on the natural logarithm of the Bayes factor: positive values of (0 − 1),(1 − 2.5),(2.5 − 5),(5 + ) correspond to inconclusive, weak, moderate, and strong evidence favoring one model over the other (T18; Hermosa Muñoz et al. 2020, but also Wheeler et al. 2017). In our case, we had the following Bayes factors: ln B_lin, flat, comparing the evidences of the two rotational models, and ln B_rot, disp, between the evidences of the best rotational model (choosing the one that has the largest Z) and that of the dispersion-only one.

We used the following priors for the kinematic parameters: $-50<({\overline{\mathit{v}}}_{\mathrm{hel}}-{\mathit{v}}_{\mathrm{g}})[\mathrm{km}{\mathrm{s}}^{-1}]<50$, where v_g is the initial mean value of the velocity distribution, 0 <  σ_v [km s⁻¹]< 50, −50 <  k [km s⁻¹ arcmin⁻¹]< 50, and −50 <  v_c [km s⁻¹]< 50. The prior over θ was set iteratively: we initially chose the prior range 0 <  θ <  π, ran the MultiNest code the first time in order to obtain the maximum value θ_m from the θ posterior distribution, and ran the MultiNest code again while updating the prior range to −π/2 <  θ − θ_m <  +π/2. This choice accounted for the limit case of θ near 0 or π.

Results of the recovered probability-weighted kinematic parameters and evidences for the three analyzed datasets, together with the effective number of probable members (defined as N_eff ≈ Σ_iP_Mi) are reported in Table 3. The P_Mi finally assigned to each target were those obtained from the most significant kinematic model.

The P91 FORS2 dataset yields 39 effective members, while the inclusion of the P69 and FLAMES data adds about 15 more effective members. These numbers already double the member stars reported by F09 (17) and G19 (36, although when accounting for their probability of membership they reduce to ∼20 effective members). We shall note that for each analyzed dataset, the effective number of members resulted to be approximately the same for each kinematic model. In all cases, the systemic velocity is stable around 180 km s⁻¹, while the velocity dispersion converges to 6 km s⁻¹ for the dispersion-only model, which proved to be the most statistically significant one. Both the systemic velocity and velocity dispersion values are significantly lower than what has been reported in literature studies. We refer to Sect. 5 for the full discussion of the results from the kinematic analysis.

We should also note that the adopted method does not take into account the different selection functions of the various datasets that were independently built. Therefore, to test whether this introduced a bias in our results, we repeated the analysis, relaxing the assumption of an exponentially declining surface-number density profile for Tucana’s stars, and simply required it to be monotonically decreasing (like in Walker et al. 2009b). No significant difference in the results was found.

5. Kinematic results

The analyses of the properties of the stellar component of Tucana carried out in the literature indicated a system with a relatively high l.o.s. velocity dispersion (σ_v  ∼  15 km s⁻¹) compared to other similarly luminous companions (see e.g., the compilation for LG dwarfs of Kirby et al. 2014; Wheeler et al. 2017). The works of F09 and G19 also reported the tentative presence of a velocity gradient likely due to internal rotation, since perspective effects related to transverse motion are negligible at the distance of Tucana. However, in a Bayesian analysis of the rotational support of the stellar component of LG galaxies using data from the literature, Wheeler et al. (2017) found no significant evidences for rotation and a l.o.s. velocity dispersion as high as 21 km s⁻¹ when analyzing the F09 catalog. Finally, there seems to be little agreement about the systemic velocity of Tucana among F09 and G19, showing an offset of ∼20 km s⁻¹ at more than 3-σ significance.

As shown in Table 3, our results are remarkably similar for all the cases we analyzed. When using the P91 FORS2 data alone, combining P91 and P69, or with all the datasets together, the systemic velocity and the velocity dispersion settled around ∼180 km s⁻¹ and 6 km s⁻¹, respectively. Furthermore, we found no significant evidence of rotation, with the dispersion-only model moderately favored in all of the cases. We note a slight increase in the velocity dispersion for the rotational models: this is caused by few targets acquiring a higher membership probability, due to the different fitted models. However, the effect is small, and all the velocity dispersion values we obtained are compatible within 1-σ.

We further performed a simpler analysis on the three datasets, considering only those targets with the highest membership probabilities (i.e., with P_M >  0.95) looking for the kinematic parameters of a dispersion-only model, whose results confirmed those of the probability-weighted analysis. Performing the same test by adding those stars with a lower membership probability (i.e., with 0.05 <  P_M <  0.95), would instead increase the velocity dispersion up to 8 km s⁻¹, which is still within 1-σ of the previous results taking into account the error bars. We also note that these stars are found further away from the center of Tucana compared to the more probable members (see Fig. 3), but they are still inside the tidal radius of the galaxy. However, it is hard to discern whether the increase in the velocity dispersion seen when including stars with lower membership probability is caused by a radially increasing velocity-dispersion profile, or simply because (all or part of) these stars are contaminants. Considering those stars with a lower membership probability that have metallicity measurements, they seem to be preferentially metal-poor (see Fig. 4). Therefore, it may also be that the increase in σ_v could be caused by the preferential inclusion of metal-poor stars with a hotter velocity dispersion than the more metal-rich stars (Tolstoy et al. 2004; Battaglia et al. 2006, 2008, 2011; Amorisco & Evans 2012a). We refer to Sect. 6.2 for more details on this point.

Fig. 3. Line-of-sight velocity distributions of the stars from the combined FORS2 + FLAMES dataset, having membership probabilities P >  0.95 (in black) and 0.05 <  P <  0.95 (in red). The red solid line indicates the systemic velocity obtained from the kinematic analysis. Left panel: distribution along the optical major axis; right panel: along the minor axis.

Fig. 4. Metallicity versus heliocentric velocity for the observed target of the combined FORS2 dataset. Black dots represent the targets with P >  0.95, while red dots show those with 0.05 <  P <  0.95.

We highlight the fact that underestimated (optimistic) or overestimated (pessimistic) velocity errors may impact on the measured velocity dispersion. We refer the reader to Fig. 1 in Koposov et al. (2011) for an analysis of the impact of underestimated/overestimated velocity errors as a function of the ratio between the true error and the true velocity dispersion. However, as reported in Appendix A, we conducted several consistency tests where we show that the velocity errors are well-determined.

Therefore, our results for Tucana point to a value for the velocity dispersion which is very unlikely to exceed 10 km s⁻¹. We assume ${\sigma }_{\mathrm{v}}=6.2_{-1.3}^{+1.6}$ km s⁻¹as our reference value, averaging between the results of the dispersion-only model from the analyzed datasets.

5.1. MultiNest mock tests

Although we have found that the rotation signal in our catalog is not statistically significant, we conducted a series of mock tests in order to explore which rotational properties can be detected according to the characteristics of our data. To this aim, we followed the approach already introduced in T18 and Hermosa Muñoz et al. (2020): we produced mock catalogs of l.o.s. velocities assuming the same number, spatial position, velocity distribution, and velocity uncertainties of the observed data. Our base catalog was the combined FORS2 P91 + P69 + FLAMES dataset after applying a P_M >  0.05 cut; the inclusion of less probable members was a compromise to have the highest number of targets (57) within the largest spatial area. To each target we assigned a mock velocity v_mock randomly extracted from a Gaussian distribution centered on zero and of standard deviation equal to the assumed velocity dispersion, fixed at σ_v, mock  =  6 km s⁻¹. This value of σ_v, mock was set according to the converging results from the probability-weighted analysis of our datasets. We further added a projected linear rotational component v_rot, mock, such that v_rot, mock/σ_v, mock = n = 1.5, 1.0, 0.75, 0.5, 0.25, 0 at the half-light radius. These correspond to velocity gradients of k = 11.2, 7.5, 5.6, 3.7, 1.9, 0 km s⁻¹ arcmin⁻¹, which we simulated at three different position angles, starting from the PA of the optical semi-major axis (97°) and then adding 45 and 90° (optical semi-minor axis). We chose an underlying linear rotation component since it resulted in higher evidence compared to a flat rotation model when analyzing our data. Each case was simulated N = 1000 times, in which we run our Bayesian kinematic analysis fitting just the linear rotation and the dispersion-only models in order to recover the related parameters and evidences. Results are shown in Table B.2.

Results from the tests indicate that the linear rotation would be spotted with high significance for velocity gradient values ≥5.6 km s⁻¹ arcmin⁻¹ aligned with the projected optical major axis. If the underlying rotation is instead milder, for example with velocity gradients ≤3.7 km s⁻¹ arcmin⁻¹, the recovered evidences for rotation are weak, inconclusive, or favoring the dispersion-only model, as we move through decreasing values of k and through different position angles. In any case, it is evident that if Tucana has a weak rotation signal (i.e., k ≤ 3.7 km s⁻¹ arcmin⁻¹), with the data we had, we could not detect it with high significance, particularly if it is not aligned with the optical major axis. If instead Tucana has a velocity gradient as reported in F09 and G19, meaning with a value close to k = 5.6 km s⁻¹ arcmin⁻¹ along the major axis, with our data we would have detected it with high significance, but this was not the case. Therefore, it seems that if rotation is actually present in Tucana, it is probably at a level of v_rot/σ_v ≲ 0.5, and to detect it with high significance a better sampling would be needed, in particular at radii around 3 ≲ R/R_e ≲ 5. We note that these results are conservative: in fact, if we had considered an input velocity dispersion as high as 10 km s⁻¹, the linear rotation signal would have been even more difficult to detect.

5.2. Comparison with other works

In Appendix A.3, we perform a comparative analysis of the velocity measurements for the individual stars derived in this work and those in F09 and G19. Here, we focus instead on the comparison of the recovered velocity parameters from the kinematic analyses of the different works.

Firstly, we ran our code on the F09 velocity measurements for member stars, and the same for G19, in order to see if we were able to recover their results. Considering a dispersion-only model, we found a 1-σ agreement between the recovered velocity parameters and those reported by these authors. Therefore, our procedure does not introduce a bias, and we can directly compare with the reported values in F09 and G19.

We found an offset of ∼15 (35) km s⁻¹ for the systemic velocity between the values reported by F09 (G19) and by us – ${\overline{\mathit{v}}}_{\mathrm{hel},\mathrm{F}09}=194.0\pm 4.3$ km s⁻¹ and ${\overline{\mathit{v}}}_{\mathrm{hel},\mathrm{G}19}=216.7_{-2.8}^{+2.9}$ km s⁻¹. These offsets are somewhat higher but still compatible with those reported in Appendix A.3, so we refer the reader to that section for an analysis of the possible causes. We stress that if we analyze our reduction of the P69 and the FLAMES data on their own, we obtain systemic velocities compatible at the 1-σ level with our value of ∼180 km s⁻¹; therefore, the differences encountered appear to be related to the treatment of the datasets. On the other hand, the velocity dispersion values (without accounting for the presence of possible gradients) reported by F09 and G19 – ${\sigma }_{\mathrm{v},\mathrm{F}09}=15.8_{-3.1}^{+4.1}$ km s⁻¹ and ${\sigma }_{\mathrm{v},\mathrm{G}19}=14.4_{-2.3}^{+2.8}$ km s⁻¹, differ from our reference σ_v value almost at the 3-σ level.

If we were to analyze our reduction of the P69 FORS2 dataset alone, we would obtain 17 effective members out of 23 input targets, showing a velocity dispersion value of ${\sigma }_{\mathrm{v},\mathrm{P}69}=11.1_{-2.7}^{+3.7}$ km s⁻¹, associated with a highly significant rotation signal. We noticed, however, that this gradient is driven by just two targets that have a very low membership probability when analyzing the combined FORS2 dataset. If we excluded them, the velocity dispersion would drop to ∼8 km s⁻¹ and the velocity gradient basically disappears. Therefore, it seems here that the low-number statistics strongly limit the conclusions we could get on the kinematic status of Tucana via the P69 dataset on its own.

The comparison with G19 results are even more puzzling. If we were to apply the same exercises as for the P69 dataset, analyzing the FLAMES data on their own, we would obtain 11 effective members out of 58 targets, which, however, could not resolve the σ_v value. This is probably due to the combination of a small number of effective members and the fact that the average velocity error of stars with a high probability of membership (δ_v ∼ 6 km s⁻¹) is comparable to the σ_v value we find in our main kinematic analysis.

Furthermore, we want to underline that we found several differences when comparing our velocity measurements to those of G19, as described in Appendix A.3. We suspect that the way the data were actually sky subtracted has led to an excess of stars with velocities around 220 km s⁻¹ in the G19 work, probably due to a combination of sky-line residuals and low S/N, which would have created fake CaT features.

5.3. Implications for Tucana’s dark matter halo properties

As previously discussed, the analysis of the internal kinematic properties of Tucana yields consistent results across the combination of datasets. Our best value for the velocity dispersion of ${\sigma }_{\mathrm{v}}=6.2_{-1.3}^{+1.6}$ km s⁻¹ is significantly lower (∼3-σ) than what was reported in the literature by both F09 and G19 (see also the discussion in the previous section), but closer now to the values observed for other similarly luminous dwarf galaxies of the Local Group (see e.g. Kirby et al. 2014; Revaz & Jablonka 2018).

We use the Wolf et al. (2010) mass-estimator, valid for dispersion-supported spherical systems, to calculate Tucana’s dynamical mass inside the half-light radius, $M_{1/2}=3G^{-1}{\sigma }_{\mathrm{v}}^2{r}_{1/2}\approx 4G^{-1}{\sigma }_{\mathrm{v}}^2{R}_{\mathrm{e}}$, where G is the gravitational constant and r_1/2 is the 3D de-projected half-light radius, which can be approximated to 4/3 R_e. Using the values from Table 1 and substituting for σ_v, we obtain $M_{1/2}=0.7_{-0.3}^{+0.4}\times {10}^7{M}_{\odot }$, which corresponds to a mass-to-light ratio within the half-light radius of $M_{1/2}/L_V={13}_{-7}^{+8}{M}_{\odot }/L_{\odot }$, assuming a luminosity of L_V = 5.5 ± 1.5 × 10⁵ L_⊙ (adapted from Saviane et al. 1996).

Using instead the velocity dispersion values from F09 and G19, applying the Wolf et al. (2010) mass-estimator we would obtain $M_{1/2,\mathrm{F}09}=4.8_{-2.0}^{+2.6}\times {10}^7{M}_{\odot }$ and $M_{1/2,\mathrm{G}19}=3.9_{-1.4}^{+1.6}\times {10}^7{M}_{\odot }$ respectively, which are more than four time as large as our own best estimation of Tucana’s mass.

Our measurements for velocity dispersion and dynamical mass of Tucana provide a new perspective to the discussion found in the literature. If instead of M_1/2 we use the value of the circular velocity $V_{\mathrm{circ}}(r_{1/2})=\sqrt{3}{\sigma }_{\mathrm{v}}$, we obtain a value of ${11}_{-2}^{+3}$ km s⁻¹ for Tucana, assuming our reference value for σ_v. This V_circ(r_1/2) is comparable to those of other similarly luminous dwarfs like Carina, Sextans, and Leo II, as it is possible to see from Fig. 10 in G19, but also from Fig. 1 in Boylan-Kolchin et al. (2012). From this last reference, we find that Tucana should live in a NFW dark matter halo with a maximum circular velocity V_max ≤ 24 km s⁻¹, as is the case for many dSphs of the LG.

In Fig. 5, we show the position of Tucana on the V_circ − r_1/2 plane with respect to the Local Group compilation of Kirby et al. (2014) and the updated value for Cetus from Taibi et al. (2018). Navarro-Frank-White halos sampled from the halo mass concentration relation of Dutton & Macciò (2014) are also shown (see Appendix B.2 for sampling details) and color-coded by V_max. The updated velocity dispersion measurements for Tucana and Cetus place these galaxies in a locus occupied by comparable luminosity dwarf galaxies. These results would imply not only that Tucana does not reside in a very centrally dense halo as predicted by G19, but also that this galaxy it is not an exception to the too-big-to-fail problem (see e.g., Kirby et al. 2014)

Fig. 5. Circular velocity at half light radius (under assumption of flat σ profile, isotropy, and sphericity) for Local Group galaxies from Kirby et al. (2014) and the isolated dwarfs Cetus and Tucana, from Taibi et al. (2018) and this work. Circular velocity profiles for NFW halos between 8.5 ≤ log(Mvir/M⊙)≤11 and following a mass concentration relation from Dutton & Macciò (2014) are shown color-coded by Vmax. Local Group galaxies represented as gray filled circles. Cetus and Tucana are highlighted with a black triangle and square, respectively; they occupy locations comparable to other dwarfs of similar stellar mass (see also Fig. B.2).

Cetus and Tucana’s isolation and well measured SFHs offer an intriguing leverage to potentially separate the multitude of solutions to the too-big-to-fail problem. In particular, what can be expected from any environmental or SFH dependence on the galaxy density profile in baryonic feedback scenarios (see discussion in Read et al. 2019), and how this could contrast with self-interacting dark matter solutions, which could predict more homogeneous behavior among all dwarfs. Further assessing these scenarios in light of the results for the isolated dSphs, and their spatial stellar population distributions will be the focus of a follow-up project.

6. Chemical analysis

The analysis of the [Fe/H] values of the FORS2 combined dataset led to the following results: considering the values with P_M >  0.95, we obtained median [Fe/H] = −1.58 dex, σ_MAD = 0.47 dex, σ_intrinsic = 0.39 dex; while adding the stars with 0.05 <  P <  0.95 we instead obtained a median of [Fe/H] = −1.61 dex, σ_MAD = 0.48 dex, σ_intrinsic = 0.39 dex. Therefore, Tucana is a metal-poor system with a significant spread in metallicity. The median [Fe/H] value measured for the likely members is in very good agreement with the integrated quantity derived from Tucana’s SFH: ⟨[Fe/H]⟩  =   − 1.52 ± 0.07 dex (Monelli et al. 2010). In addition, our average [Fe/H] value falls within the rms scatter of the stellar luminosity-metallicity relation for LG dwarf galaxies reported by Kirby et al. (2013a), while the intrinsic scatter agrees well with the values of other similarly luminous dwarf galaxies (Leaman et al. 2013).

The distribution of the [Fe/H] values as a function of the elliptical radius is shown in Fig. 6. The significant scatter in metallicity toward the inner part of the galaxy is evident. In addition, there is a bimodality in the metallicity histogram suggesting the presence of two subpopulations. This can be related to the results obtained from deep photometric data (Monelli et al. 2010), where it was shown that the splitting observed in the HB, in the RGB bump, and in the properties of the RR-Lyrae stars, imply that Tucana has been able to produce a second generation of metal-rich stars thanks to the self-enrichment from the first stars in a very short period of star formation (∼1 Gyr). In a more recent study, Savino et al. (2019) reanalyzed the photometric data from Monelli et al. (2010), by studying the main-sequence turn off and the HB of Tucana. They were able to obtain an SFH with a finer temporal resolution, showing that Tucana actually experienced two early phases of star formation (SF), followed by a third one ending between 6 and 8 Gyr ago, with the two initial episodes being the most intensive. According to the age-metallicity relation recovered by Savino et al. (2019), our [Fe/H] measurements could be related to the intermediate-old and intermediate-young age populations of the two last episodes of SF. This would explain the bimodality we observed in the metallicity histogram, while our most metal-poor stars could be related to the oldest episode of SF. However, despite the relatively high intensity of this SF period, we found very few stars with [Fe/H] < −2.25 dex. This is probably related to the fact that, for lower metallicity stars, due to weaker lines, a higher S/N is needed in order to get similar accuracy in [Fe/H] measurements. Furthermore, their spatial distribution tend to be more extended, which requires extra attention during the sampling phase.

Fig. 6. [Fe/H] as a function of the elliptical radius scaled with Re for Tucana’s probable member stars from the FORS2 combined dataset. Black dots represent the targets with membership probabilities P >  0.95, while the red dots show those with 0.05 <  P <  0.95. 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 gray band indicates the corresponding 1-σ confidence interval. The solid red line and the red band indicate the same, but using all targets with P >  0.05. The histogram on the right side represents the metallicity distribution of the stars with P >  0.95.

6.1. Looking for a metallicity gradient

Observations have shown that in many LG dwarf galaxies the young and metal-rich stars are more spatially concentrated than the old and metal-poor ones that display a more extended spatial distribution. Their overall radial distribution produces a decreasing metallicity gradient, e.g., Fornax (Battaglia et al. 2006; Leaman et al. 2013), Phoenix (Kacharov et al. 2017), which could eventually reach a plateau on the outside, e.g., Sculptor (Tolstoy et al. 2004), VV 124 (Kirby et al. 2013b), Cetus (T18). These results have also been reproduced by simulations (e.g., Schroyen et al. 2013; Revaz & Jablonka 2018), which have shown that the shape of these gradients strongly depends on the combination of the stellar mass, SFH, and dynamical history of the system under consideration, although their strength could be influenced by merger events (Benítez-Llambay et al. 2016) or environmental effects such as tidal stripping (Sales et al. 2010).

Therefore, we investigated the presence of a metallicity gradient as a function of radius, first focusing on the stars with P >  0.95, which extended up to ∼3 R_e (see Fig. 6). Performing an error-weighted linear least-square fit to the data, we obtained the value $m=\frac{\mathrm{d}[\mathrm{Fe}/\mathrm{H}]}{\mathrm{d}R}=-0.16\pm 0.09$ dex arcmin⁻¹ (= − 0.6 ± 0.4 dex kpc⁻¹ = −0.13 ± 0.07 dex $R_{\mathrm{e}}^{-1}$, using the values reported in Table 1 for the conversions). We also performed a Gaussian process regression (GPR) analysis, where we used a Gaussian kernel together with a noise component to account for the intrinsic metallicity scatter. The GPR has the advantage of being a kernel-based non-parametric probabilistic method that makes it possible to compute empirical confidence intervals. Since we are looking for a smooth function, it performs better than a least-square fit in finding the general trend in the data. In our case, we confirmed the decreasing trend, although the 1-σ confidence limits proved quite large due to the high intrinsic scatter of the data, making the presence of a metallicity gradient within ∼3 R_e dubious.

We further checked this result by performing a simple simulation. We assumed a double Gaussian metallicity distribution with parameters roughly fitting the observed one, but no spatial variation (μ_[Fe/H],1 = −2.0 dex, σ_[Fe/H],1 = 0.2, μ_[Fe/H],2 = −1.3 dex, σ_[Fe/H],2 = 0.2, assuming the same fraction of stars in the two Gaussians). We then randomly extracted [Fe/H] values at the radial positions of our data. We further reshuffled the [Fe/H] values according to the observed errors, and, finally, we performed a linear least-square fit looking for a spatial metallicity gradient. We repeated this process 1000 times. The obtained average gradient was compatible with zero, with the associated scatter large enough to include within 1-σ the observed value of m. Therefore, with the data we had, the observed gradient within R <  3 R_e is not statistically significant.

Adding the stars with 0.05 <  P <  0.95 would extend the spatial coverage up to R ∼ 6 R_e, thanks to the two outermost targets, but would lead to an even milder gradient: m = −0.07 ± 0.04 dex arcmin⁻¹ (= − 0.28 ± 0.16 dex kpc⁻¹ = −0.06 ± 0.03 dex $R_{\mathrm{e}}^{-1}$), by performing a linear least-square fit. The presence of a metallicity gradient in Tucana is expected from studies of deep-photometric data (Hidalgo et al. 2013; Savino et al. 2019). However it is probable that we are mainly targeting stars belonging to the more recent episodes of SF, whose populations share similar spatial extensions (see e.g., Fig. 11b in Savino et al. 2019). Therefore, the presence of a metallicity gradient in Tucana is very tentative, and we would need a better sampling, in particular, of the metal-poor component, around 3 ≲ R/R_e ≲ 5 to put our results on firmer ground.

We compared Tucana’s metallicity gradient (or rather the lack thereof) with those of some MW satellites with similar luminosities (L_V ≲ 5 × 10⁵ L_⊙) and short SFHs, for example, Draco (Aparicio et al. 2001), Ursa Minor (Carrera et al. 2002), and Sextans (Bettinelli et al. 2018). All of them formed the majority of their stars more than 10 Gyr ago within a short period of SF, which in some cases may have lasted no more than 1 Gyr (i.e., Sextans). It has been shown that such short SFHs may lead to mild metallicity gradients in these systems (Marcolini et al. 2008; Kirby et al. 2011; Revaz & Jablonka 2018). Indeed, both Draco and Ursa Minor show mild gradients with values of −0.05 dex $R_{\mathrm{e}}^{-1}$ and −0.03 dex $R_{\mathrm{e}}^{-1}$, as reported by Schroyen et al. (2013) and Kirby et al. (2011), respectively⁷. Sextans, on the other hand, seems to have a stronger gradient of −0.24 dex $R_{\mathrm{e}}^{-1}$, as reported by Schroyen et al. (2013) using the Battaglia et al. (2011) spectroscopic dataset. However, this value is somewhat overestimated since recent studies of the structural properties of Sextans (Roderick et al. 2016; Cicuéndez et al. 2018) show that this system is less extended than what was previously reported in the literature. Using the half-light radius value from Cicuéndez et al. (2018), we find a lower gradient of −0.18 dex $R_{\mathrm{e}}^{-1}$, which is still far from the other dwarf’s values and probably related to an early merger event that could have steepened its metallicity gradient (see Cicuéndez & Battaglia 2018, but also Benítez-Llambay et al. 2016). If the case of Tucana is similar to those of Draco and Ursa Minor, as it seems, we would expect it to host an equally mild metallicity gradient, but it would take a better sampling of the spatial extension of the metal-poor component in Tucana to confirm it.

6.2. Searching for two chemo-kinematically-distinct populations

Although we have not spotted a clear metallicity gradient in Tucana, the bimodality found in the metallicity distribution, may indicate the presence of two subpopulations which differ not only in their chemical properties, but also in their kinematics. Some of the dSphs satellites of the MW, such as Sculptor, Fornax, Carina, and Sextans, show such features where the metal-rich (usually more spatially concentrated) subpopulation has colder kinematics than the metal-poor (and more extended) one (e.g., Tolstoy et al. 2004; Battaglia et al. 2006, 2008, 2011; Koch et al. 2008; Amorisco & Evans 2012a). Determining the chemo-kinematic properties of dSphs is of great interest, not only to better understand their evolutionary path, but also to get an insight into their dark matter properties (e.g., Battaglia et al. 2008; Walker & Peñarrubia 2011; Amorisco & Evans 2012b; Strigari et al. 2018).

In the case of Tucana, we first analyzed the combined FORS2 dataset with the P >  0.95 cut applied (see Fig. 4). We took the median [Fe/H] value of −1.58 dex to split our sample into metal-rich (MR) and a metal-poor (MP) subsamples. We then ran our code to obtain the kinematic parameters of both samples (see Sect. 5). Using the dispersion-only model, we found ${\sigma }_{\mathrm{v},\mathrm{MR}}=6.0_{-1.6}^{+1.9}$ km s⁻¹ and ${\sigma }_{\mathrm{v},\mathrm{MP}}=9.0_{-2.4}^{+2.8}$ km s⁻¹, which are 1-σ from each other. Including instead the 0.05 <  P <  0.95 data, we obtained: ${\sigma }_{\mathrm{v},\mathrm{MR}}=5.9_{-1.6}^{+1.9}$ km s⁻¹ and ${\sigma }_{\mathrm{v},\mathrm{MP}}=10.9_{-2.2}^{+2.6}$ km s⁻¹, which are instead ∼2-σ from each other (see Table 4).

Therefore, there is weak evidence of two chemo-kinematically distinct subpopulations in Tucana. Additional data, in particular including external parts of the galaxy are necessary to reach firm conclusions, regarding the presence of a metallicity gradient as well as the possible distinct chemo-kinematic populations.

7. Summary and conclusions

In this paper, we present results regarding the internal kinematic and metallicity properties of the Tucana dwarf spheroidal galaxy, based on the analysis of multi-object spectroscopic samples of individual RGB stars.

This analysis is based on a novel set of 50 individual objects collected with the VLT/FORS2 instrument in MXU mode in P91, complemented by a re-reduction and reanalysis of two datasets from the literature, namely the VLT/FORS2-MXU dataset presented in Fraternali et al. (2009, F09) and the VLT/FLAMES-GIRAFFE one by Gregory et al. (2019, G19). Applying a probabilistic membership approach, we find 39 effective members in our P91 sample, which doubles the number of Tucana’s member stars found in F09 and G19.

A full re-reduction and analysis of the data presented in the literature was carried out because it became clear that the published catalogs could not be directly combined with the l.o.s. velocities we derived for the P91 sample: there are significant differences between the values of the systemic velocity reported in those studies with respect to that we derived from the P91 dataset (∼195 km s⁻¹ for F09, ∼215 km s⁻¹ for G19, while ∼180 km s⁻¹ in our case); and the comparison of the individual l.o.s. velocities for the stars in common were supporting the presence of shifts with respect to F09, but were not sufficient to fully quantify whether that was the only source of difference, or if it was even more unfavorable for the comparison with G19.

Following our homogeneous data reduction, we find an excellent agreement between velocity measurements of the three datasets both for stars in common (Fig. A.1 top row panels) as well as for systemic velocity, which is stable around 180 km s⁻¹ for the three datasets (see Fig. 2, right panel, and Table 3). We proceeded to analyze the P91 dataset alone and also in combination with our treatment of the P69 and FLAMES data. The resulting values of the intrinsic l.o.s. velocity dispersion are consistently around 6 km s⁻¹ when considering the three combinations of datasets (P91, P91 + P69, all the three combined) and the highly probable members (probability of membership P_M >  0.95); when including lower probability members, the velocity dispersion increases, but it is unlikely that σ_v is higher than 10 km s⁻¹. Therefore, our analysis leads to the conclusion that the l.o.s. velocity dispersion of Tucana’s stellar component is much lower than the values reported by F09 and G19 – ${\sigma }_{\mathrm{v},\mathrm{F}09}=15.8_{-3.1}^{+4.1}$ km s⁻¹ and ${\sigma }_{\mathrm{v},\mathrm{G}19}=14.4_{-2.3}^{+2.8}$ km s⁻¹.

Furthermore, we find no significant signs of internal rotation. Mock tests suggest that if Tucana had had a maximum rotational velocity of ∼10−15 km s⁻¹ along the projected major axis (as previously reported in literature) with the data available to us, we should have detected it with high significance. On the other hand, lower levels of rotation are not completely ruled out, but a larger sample would be needed to quantify their presence. Nevertheless, it seems improbable that Tucana is a fast rotator (v_rot/σ_v ≳ 1).

Assuming an average ${\sigma }_{\mathrm{v}}=6.2_{-1.3}^{+1.6}$ km s⁻¹for our data, we obtain a dynamical mass within the half-light radius of $M_{1/2}=0.7_{-0.3}^{+0.4}\times {10}^7{M}_{\odot }$. This translates into a circular velocity at the half-light radius of $V_{\mathrm{circ}}(r_{1/2})={11}_{-2}^{+3}$ km s⁻¹, which implies that, if Tucana inhabits a NFW dark matter halo, it should have a similar density to those of other MW dSphs (Boylan-Kolchin et al. 2012). Therefore, Tucana is not an exception to the too-big-to-fail problem and not “a massive failure”, as it had gained fame for being.

The analysis of Tucana’s chemical properties was only carried out on the combined P91 and P69 FORS2 data, due to their higher S/N. We establish that the galaxy is mainly metal-poor with a significant scatter in metallicity (having a median [Fe/H] = −1.58 dex, σ_MAD = 0.47 dex and σ_intrinsic = 0.39 dex when considering only highly likely members, and median [Fe/H] = −1.61 dex, σ_MAD = 0.48 dex, and σ_intrinsic = 0.39 dex when including the less likely members). The derived values agree very well with SFH studies (Monelli et al. 2010; Savino et al. 2019). In addition, the average [Fe/H] falls between the rms scatter of the stellar luminosity-metallicity relation for LG dwarf galaxies (Kirby et al. 2013a).

Looking at the distribution of the [Fe/H] values as a function of radius, we find a mild metallicity gradient. However, the size and spatial distribution of the current datasets do not lead to a statistically significant detection. The presence of a gradient in Tucana would be expected from the age gradients inferred from deep photometric studies (Hidalgo et al. 2013; Savino et al. 2019), but also from both observations and simulations of similarly luminous dwarfs (see e.g., Leaman et al. 2013; Schroyen et al. 2013; Revaz & Jablonka 2018), which indeed host mild metallicity gradients. Therefore, the presence of an underling gradient in Tucana is not excluded, but it would need a better sampling of the metal-poor component in Tucana (particularly at R ≳ 3 R_e) to confirm it. Finally, we find a hint of the presence of multiple stellar populations with distinct chemo-kinematical properties, although in this case, the addition of new data would also help to put this result on firmer ground.

Acknowledgments

We would like to thank the anonymous referee for the helpful and constructive comments. We thank Giacomo Beccari and Michele Bellazzini for sharing with us their VLT/VIMOS photometric data. We also thank F. Fraternali, E. Tolstoy, A. Gregory and M. Collins for useful discussions in the different phases of this project. S.T. thankfully acknowledges ESO for a one month visit funded by an SSDF grant. G.B. and S.T. acknowledge financial support through the grants (AEI/FEDER, UE) AYA2017-89076-P, AYA2014-56795-P, as well as by the Ministerio de Ciencia, Innovación y Universidades (MCIU), through the State Budget and by the Consejería de Economía, Industria, Comercio y Conocimiento of the Canary Islands Autonomous Community, through the Regional Budget. S.T. acknowledges an FPI fellowship BES-2015-074765, while G.B. acknowledges financial support through the grant RYC-2012-11537. This research has made use of NASA’s Astrophysics Data System and extensive use of IRAF, Python, Numpy, Scipy and Astropy ecosystems.

References

Appendix A: Line-of-sight velocity and metallicity measurements – consistency checks and comparisons

A.1. Consistency checks

We performed consistency checks in order to unveil the presence of possible systematic errors in our velocity determinations using fxcor, or other issues. For all the datasets analyzed, we searched for systematics by shifting the template spectrum at several velocities (from −50 km s⁻¹ to 500 km s⁻¹ at step of 50 km s⁻¹) and cross-correlating them with the templates themselves at rest, through the full wavelength range of data, and also only around the CaT. We did not find any significant systematic shift introduced by the cross-correlation procedure as a function of the velocity shift nor the cross-correlation range.

We also verified that the random errors were well-treated by the fxcor task and if they led to reliable velocity errors. In this case, we performed a Monte Carlo analysis, randomly adding the error-spectra from targets having measured S/N of 10, 15, 20, 25, and 30 Å⁻¹ for the FORS2 datasets, and of 5, 10, 20, 50, and 100 Å⁻¹ for the FLAMES one, to the shifted templates of the previous step. We performed 250 realizations for each case and found that the procedure tends to produce velocity errors that are underestimated with respect to the velocity scatter from the individual Monte Carlo runs for S/N lower than 10 Å⁻¹, for both datasets. Therefore, we take these values as the limit S/N above which we will trust the velocity error estimations. At this S/N, the velocity errors were found to be ∼ ± 10 km s⁻¹ in both cases.

We further checked the choice of using the two reddest lines of the CaT for the cross-correlation measurements of the FLAMES dataset (see previous section), creating two more catalogs of velocity estimations using both the entire CaT range and the two bluest lines. Cross-correlating between catalogs and selecting those targets with an S/N >  10 Å⁻¹, which reduced our sample to 76 objects out of 164, we found very good agreement, excluding eight targets that had velocity measurement differences more than 3-σ from each other. We visually checked the corresponding spectra for these targets and found that for the outliers the first and/or third line of the CaT resulted affected by sky residuals or were hidden by noise, while in general, the central line resulted more clearly visible. For fiber 118, the opposite happened, and we considered the measurement obtained using the full CaT range for this case.

Considering once again the velocities from targets having an S/N > 10 Å⁻¹, we looked for cases were the cross-correlation procedure did not get a velocity solution, finding only two cases (fibers 18 and 55). Their visual inspection did not show obvious CaT lines, and we thus decided to discard these two targets.

Finally, we checked the internal accuracy of our P91 FORS2 dataset by inspecting the four targets in common between the two pointings. Of these, three had reliable measurements of velocity and metallicity that were compatible within 1-σ, while the other one did not have reliable measurements. In fact, it did not show a visible CaT, and its I-magnitude and color were also not compatible with the RGB of Tucana (see also Sect. 4). We then decided to discard this target from the sample and to average the measurements from the other three targets together.

We further excluded two more targets from the P91 FORS2 dataset since their measurements were not reliable due to high residuals in their spectra: these are slit 14, chip-2 of Tuc0 field, and slit 17, chip-1 of Tuc1. The high residuals in these two spectra were due to the stellar trace falling on a bad CCD-row and a poor sky-subtraction that could not be improved, respectively.

A.2. Comparison between measurements of different datasets obtained with our reduction

We compared the measurements obtained from the P91 dataset with those from our reduction of the P69 and FLAMES datasets, using the targets they had in common. This was an important step to ensure we could combine the velocity catalogs.

The P91 and P69 FORS2 datasets have five targets in common with measured velocities and metallicities, while the FLAMES dataset has seven targets in common (four taking those with S/N >  10 Å⁻¹) with the P91 FORS2 dataset and 18 (15) with the P69 one (see the top row of Fig. A.1 for the velocity comparison).

Fig. A.1. Velocity comparison for common targets between datasets – vhel vs. Δvhel. Top: between the datasets from our internal reduction – left: FORS2 datasets; center: flames vs. P91-FORS2; right: flames vs. P69-FORS2. Red circles indicate the problematic targets highlighted in the main text. Bottom: between our datasets and those from the literature – left: F09 vs. P69-FORS2; central: G19 vs. P69-FORS2; right: G19 vs. flames.

The agreement between these common targets is excellent, mostly within 1-σ, except for three targets in common between the FLAMES dataset and the P69 FORS2 one (marked as object 1, 2 and 3 in the top row of Fig. A.1) with velocities over 3-σ from each other. Of these targets, two are in common with the P91 FORS2 dataset (objects 1 and 3), which are in good agreement with the FLAMES measurements (Fig. A.1 top row, central panel), but not with those of the P69 FORS2 dataset (Fig. A.1 top row, right panel). We then found that the problem may reside with the measurements from the P69 FORS2 catalog. Therefore, we inspected these outliers, namely slits 9 and 10 of chip-1 and slit 13 of chip-2 (respectively object 1, 2, and 3 in the cited figures), finding that slits 9 and 13 have particularly noisy sky-line residuals around the first two CaT lines that may have compromised their measurements, while the spectrum of slit 10 was clean. However, this target was marked as a double star and with slit mis-centering problems in F09. We then decided to exclude these values from the P69 FORS2 catalog, but to maintain them in the P91 FORS2 and FLAMES datasets. In conclusion, after a homogeneous analysis, the measurements for stars in common, which were taken during several epochs as well as with different instruments and spectral resolution in the case of FORS2 and FLAMES, are found to be in very good agreement.

A.3. Comparison with other works

Comparison with F09: The l.o.s. velocity and metallicity distributions of the P91 and F09 datasets show a difference of ∼10 km s⁻¹ and ∼0.25 dex in the mean values, respectively. We investigated whether clear offsets could be found by comparing the measurements of the targets in common between the two studies. There are only three such stars, for which the metallicities agree at the 1-σ level, but the l.o.s. velocities differ at 3-σ level for two stars, and 1.5-σ for the other one, respectively. The stars discrepant at the 3-σ level were slits 9 and 13 of chip-1 and -2, respectively (targets 7 and 18 in F09), whose P91 spectra do not show any particular issue.

Since we could not pinpoint the sources of these differences in the line-of-sight velocities, we decided to re-reduce and reanalyze the F09 FORS2 dataset, so as to be homogeneous with the treatment of the P91 data.

The comparison of the l.o.s. velocities and metallicities derived from our reduction of the P69 FORS2 dataset and those published by F09 showed average offsets of 7 km s⁻¹ and 0.2 dex in the velocity and metallicity measurements, of the same signs and comparable to those previously found between the targets in common between the P91 FORS2 dataset and the F09 catalog (see Fig. A.1 bottom row, left panel for the velocity comparison). It is unclear what leads to the measured difference in velocity. For the metallicities, the [Fe/H]-CaT EW calibration in F09 is equivalent to that used here in the metallicity and magnitude regime of Tucana’s stars, and the 0.1 mag difference in the adopted V_HB value does not have a significant impact on the final [Fe/H] values. The main source of difference is likely to reside in how the EWs were calculated: in F09, Gaussian profiles were used to fit the CaT lines, while we adopt Voigt profiles. We checked that if we were to adopt Gaussian profiles, we would obtain lower metallicities of a difference comparable to the systematic shift we previously found. We want to stress, however, that the Voigt profiles represent a better fit to the CaT lines, in particular along the line wings (see e.g., Starkenburg et al. 2010).

Comparison with G19: We compared our own reduced datasets with the full catalog of observed targets by G19. There was not a good agreement in all cases. The P91 FORS2 dataset was the one with the lowest number of common targets (6), which showed mainly a systematic offset of ∼30 km s⁻¹. One target, however, showed a velocity difference of ∼130 km s⁻¹. We visually checked for it in our dataset (target 13 of chip-1 in Tuc0), finding a clean, high S/N spectrum, with no features that could have biased the velocity measurement.

The P69 FORS2 dataset proved to have 16 targets in common with G19 (see Fig. A.1 bottom row, central panel). While for a few targets there is good agreement, overall, there is a large scatter in the distribution of velocity differences with one significant outlier (target 1 of chip-2). An inspection of its corresponding spectrum in the P69 data showed the CaT lines shifted with respect to the rest frame by ∼200 km s⁻¹, compatible with what we found in the cross-correlation measurements, but significantly different from the ∼25 ± 14 km s⁻¹ reported by G19.

The comparison between the velocities we derived from our treatment of the FLAMES dataset was even more puzzling, considering that the values obtained come from the same sample. As can be seen from the bottom-right panel of Fig. A.1, there is a significant scatter among the 59 targets in common (we compared only those measurements from our dataset with a, S/N >  10 Å⁻¹), in particular at high velocities (> 100 km s⁻¹), and there might even be a dependence as a function of velocity between the two datasets. We also confirmed that these problems remain when considering only those targets marked as likely members of G19.

A systemic velocity for Tucana around 200 km s⁻¹ was found by G19. At v_hel, G19 ∼ 200 km s⁻¹, the comparison in Fig. A.1 (bottom row, right panel) gives velocities around 170 − 190 km s⁻¹ for our FLAMES measurements of the same stars (i.e., velocities close to the systemic velocity we find for Tucana), as well as values between 0 and 50 km s⁻¹, so, typical of foreground stars. Our suspicion is that the G19 catalog contains an excess of stars with velocities around 200 km s⁻¹.

We investigate this possibility in two ways: by comparing the radial distribution of the G19 members with that obtained from photometric observations of Tucana’s stars (test A) and comparing the G19 l.o.s. velocity distribution for the stars they considered as nonmembers with the expectations from a Galactic foreground model (test B). In both cases, we took into account the displacement of the FLAMES pointings from the center of Tucana.

Test A: We looked at the normalized cumulative radial distribution of G19’s member stars, and compared it to that obtained from the observed surface density profile of Tucana that we measured from the VLT/VIMOS photometric catalog introduced in Sect. 4. Although the VLT/VIMOS catalog does not cover the entire area of the FLAMES pointings, it was sufficient for following the surface-density profile of the RGB stars of Tucana up to its nominal tidal radius.

We found that the number of G19’s member stars in the outer parts of Tucana, up to where VLT/VIMOS photometry extends, tends to be overestimated with respect to that expected from the photometry for Tucana’s RGB stars in the same area. It could, however, be argued that the FLAMES/GIRAFFE fiber set-up might imprint a different distribution than that expected from the photometry. Therefore, we performed the following “Test B”, which is instead free from this possible issue.

Test B: We now concentrate on the l.o.s. velocity distribution of the stars marked as nonmembers in G19, which should only contain Galactic contaminants: this distribution shows two clear peaks in the velocity histogram, one around 0 km s⁻¹, and, unexpectedly, the other at ∼200 km s⁻¹ (G19, Fig. 6). Making the reasonable assumption that the nonmembers are mainly foreground Galactic contaminants, we compared their velocity distribution with that obtained from the Besançon model (Robin et al. 2003) generated in the direction of Tucana over an area equivalent to that of a FLAMES-GIRAFFE pointing, and by opting for the Besançon model stars to have similar positions on the CMD as on the FLAMES targets. The velocity distribution of the Besançon model stars showed just a single peak around 0 km s⁻¹, with a smooth decline towards negative and positive velocities, with a tail extending to 300 km s⁻¹. We checked if the peak at ∼200 km s⁻¹ for G19’s nonmembers could be explained from the distribution expected from the Galactic model. We randomly chose, from the synthetic dataset, a number of stars equal to that of the G19 nonmembers, and, over 1000 trials, calculated the number of objects extracted from the Besançon model that would have velocities > 150 km s⁻¹: we never got a number of contaminants as high as that of G19’s nonmembers over the same velocity range. This indicates that the number of stars with velocities ∼200 km s⁻¹ in the G19 catalog of nonmembers could likely be overestimated. Performing the same exercise for our FLAMES targets that proved to be nonmembers (i.e., with P <  0.05) considering the targets with velocities > 120 km s⁻¹ (accounting for the observed shift between the datasets), we found that the observed number of nonmembers is within the 87% (1.5-σ) of the distribution obtained for the contaminants.

We speculate that this excess of velocities around ∼200 km s⁻¹ in the G19 dataset can be attributed to a sky subtraction problem around the CaT region. In fact, the sky-lines at 8504 Å and 8548 Å, which are around the first and second lines of the CaT, if badly subtracted could lead to strong absorption residuals. These features, in a low-S/N regime and during the cross-correlation procedure, could be mistaken for the first two CaT lines shifted at ∼200 km s⁻¹, so around the value of the systemic velocity reported for Tucana by G19.

Appendix B: Supplementary material from the kinematic analysis

B.1. Supplementary tables

We provide here supplementary tables with results from the structural analysis of Tucana using the VLT/VIMOS photometric dataset, as reported in Sect. 4, along with the results from the mock tests presented in Sect. 5.1.

In all tables, the reported values of the parameters obtained from a Bayesian analysis represent the median of the corresponding marginalized posterior distributions, with 1-σ errors set as the confidence intervals around the central value enclosing 68% of each distributions.

B.2. CDM halo profiles

The NFW circular velocity curves shown in Fig. 5 are generated for 300 mock halos, which are uniformly sampled in virial mass between 8.5 ≤ log M_vir ≤ 11, and follow the redshift zero halo mass concentration relation from Dutton & Macciò (2014) with scatter of σ_lnc = 0.25. These halos are shown in Fig. B.1, color-coded by their V_max.

Fig. B.1. Distribution of NFW halos in virial mass – concentration plane used to construct circular velocity profiles in Fig. 5. Halos are color-coded by Vmax.

In Fig. B.2, to be comprehensive, we reproduce plots where the mock halos are color-coded by concentration, or virial mass. In the top panel, the Local Group dwarf galaxies are color-coded according to their stellar mass; a probable correlation is visible between the circular velocity and the stellar mass. In the bottom panel, observed dwarf galaxies with virial masses estimated from dynamical modeling (Read & Erkal 2019; Leaman et al. 2012) are also color-coded, illustrating the tension between the predicted and observed density profiles.

Fig. B.2. Same as Fig. 5, but with halo circular velocity profiles color-coded by concentration (top) and virial mass (bottom). Cetus and Tucana are shown as a triangle and a square, respectively. Top panel: the Local Group dwarf galaxies are color-coded according to their stellar mass; a probable correlation is visible between the circular velocity and the stellar mass. Bottom panel: the color-coded dots that follow the same color schema as the halo virial mass profiles are Local Group dwarf galaxies, for which virial masses were obtained from dynamical modeling (Read & Erkal 2019; Leaman et al. 2012).

Appendix C: Additional tables

All Tables

All Figures

Fig. 1. Spatial distribution (left) and color-magnitude diagram (right) of stars along the line of sight to the Tucana dSph. Black points represent the objects classified as stars in the FORS2 photometric data (see main text); red and blue dots indicate the P91-FORS2 MXU targets classified as probable members (i.e., with P >  0.05) and non-members, respectively. Yellow triangles and green squares represent the probable member stars from the P69-FORS2 and FLAMES datasets, respectively, which were added to the P91-FORS2 and analyzed through the text. The two observed FORS2 pointings are represented as large squares, while the ellipses denote the galaxy half-light radius and the spatial extension of the dataset (i.e., up to 5 × Re ∼ Rtidal). We note that the photometric data are not corrected for reddening.

In the text

Fig. 2. Histogram of l.o.s. velocity measurements from the analyzed datasets. Left: comparing the velocities from the P91-FORS2 dataset with those from literature, i.e. Fraternali et al. (2009) and Gregory et al. (2019). Right: same comparison using instead the velocity measurements from our reduction of the P69-FORS2 and FLAMES datasets. Note that, in the left panel, the peaks of the histograms fall at different velocities, while in the right panel, where the datasets shown were analyzed homogeneously, these differences are absent.

In the text

	Fig. 3. Line-of-sight velocity distributions of the stars from the combined FORS2 + FLAMES dataset, having membership probabilities P >  0.95 (in black) and 0.05 <  P <  0.95 (in red). The red solid line indicates the systemic velocity obtained from the kinematic analysis. Left panel: distribution along the optical major axis; right panel: along the minor axis.
In the text

	Fig. 4. Metallicity versus heliocentric velocity for the observed target of the combined FORS2 dataset. Black dots represent the targets with P >  0.95, while red dots show those with 0.05 <  P <  0.95.
In the text

Fig. 5. Circular velocity at half light radius (under assumption of flat σ profile, isotropy, and sphericity) for Local Group galaxies from Kirby et al. (2014) and the isolated dwarfs Cetus and Tucana, from Taibi et al. (2018) and this work. Circular velocity profiles for NFW halos between 8.5 ≤ log(Mvir/M⊙)≤11 and following a mass concentration relation from Dutton & Macciò (2014) are shown color-coded by Vmax. Local Group galaxies represented as gray filled circles. Cetus and Tucana are highlighted with a black triangle and square, respectively; they occupy locations comparable to other dwarfs of similar stellar mass (see also Fig. B.2).

In the text

Fig. 6. [Fe/H] as a function of the elliptical radius scaled with Re for Tucana’s probable member stars from the FORS2 combined dataset. Black dots represent the targets with membership probabilities P >  0.95, while the red dots show those with 0.05 <  P <  0.95. 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 gray band indicates the corresponding 1-σ confidence interval. The solid red line and the red band indicate the same, but using all targets with P >  0.05. The histogram on the right side represents the metallicity distribution of the stars with P >  0.95.

In the text

Fig. A.1. Velocity comparison for common targets between datasets – vhel vs. Δvhel. Top: between the datasets from our internal reduction – left: FORS2 datasets; center: flames vs. P91-FORS2; right: flames vs. P69-FORS2. Red circles indicate the problematic targets highlighted in the main text. Bottom: between our datasets and those from the literature – left: F09 vs. P69-FORS2; central: G19 vs. P69-FORS2; right: G19 vs. flames.

In the text

	Fig. B.1. Distribution of NFW halos in virial mass – concentration plane used to construct circular velocity profiles in Fig. 5. Halos are color-coded by Vmax.
In the text

Fig. B.2. Same as Fig. 5, but with halo circular velocity profiles color-coded by concentration (top) and virial mass (bottom). Cetus and Tucana are shown as a triangle and a square, respectively. Top panel: the Local Group dwarf galaxies are color-coded according to their stellar mass; a probable correlation is visible between the circular velocity and the stellar mass. Bottom panel: the color-coded dots that follow the same color schema as the halo virial mass profiles are Local Group dwarf galaxies, for which virial masses were obtained from dynamical modeling (Read & Erkal 2019; Leaman et al. 2012).

In the text