© The Authors 2024

1 Introduction

The characterization of the internal kinematic properties of the stellar component of galaxies is an important step for unraveling the processes that have shaped their evolution. It is also propaedeutic to the application of dynamical modeling techniques aimed at inferring the galaxy’s mass distribution.

In the Local Group (LG), we have the possibility of carrying out this type of characterization in detail, on a star-by-star basis, for a large sample of dwarf galaxies, reaching down to a regime in stellar mass that is otherwise inaccessible to current facilities. This regime is particularly interesting as it is where the galaxies with the largest dynamical mass-to-light ratios are found (Mateo 1998; Walker & Peñarrubia 2011; Simon 2019; Battaglia & Nipoti 2022, and references therein), making them very valuable probes to test cosmological models.

For LG galaxies, we can also rely on a wealth of additional, detailed information, such as the stars’ chemistry and relative ages, which allow one to explore whether there are causal links between observed properties (e.g., the impact of rotational versus velocity dispersion support on age and metallicity gradients or the potential imprints of star formation episodes on the chemokinematic properties of the stellar component). In addition, the impact of environmental effects can be assessed at the level of the individual dwarf galaxies, at least to some extent. For example, from the 3D distance from the large LG spirals, it is possible to identify what dwarf galaxies might have had a negligible level of interaction with these much more massive systems. In particular for Milky Way (MW) satellites, the availability of accurate and precise systemic proper motions derived from Gaia and Hubble Space Telescope (HST) has opened the door to the possibility of relating their observed properties to their orbital history (Sohn et al. 2017; Li et al. 2021; Pace et al. 2022; Battaglia et al. 2022), even though uncertainties in both the satellites’ systemic motions and in the mass distribution of the MW, and how it is assembled over time, imply that such considerations are most robust for the last one or two pericentric passages (e.g., D’Souza & Bell 2022; Vasiliev 2024; Santistevan et al. 2024).

The direction and amplitude of velocity gradients, combined with the information on the galaxy’s orbit, can give several useful indications of a galaxy’s past life. As a result of various observing campaigns from multiple groups (e.g., Tolstoy et al. 2004; Walker et al. 2008), it has become clear that the stellar component of LG dwarf galaxies is mainly supported by random motions (see e.g., Wheeler et al. 2017, who find that the measured ratio between ordered and random motions V/σ < 1 for three-fourths of the sample). This is also true for systems that are unlikely to have ever interacted with the large LG spirals, or that might have had experienced one pericentric passage around the MW or M31 at most (Taibi et al. 2018, 2022). Such findings disfavor, or at least call for a revision of, those models suggesting that the pressure-supported dwarf spheroidal galaxies form from tidally stirred, rotationally supported dwarf galaxies (e.g., Mayer et al. 2001, 2006; Kazantzidis et al. 2011). On the other hand, the presence of rotation along the minor-axis (prolate rotation) has been unveiled in And II (Ho et al. 2012) and Phoenix (Kacharov et al. 2017) and linked to merger events (Lokas et al. 2014; Cardona-Barrero et al. 2021). In other systems, velocity gradients aligning with the projection of the orbital path on the sky are revealing tidal stripping (e.g., AntliaII, TucanaIII) (Mutlu-Pakdil et al. 2018; Torrealba et al. 2019), and relations have been made between the velocity gradients’ amplitude and orbital phase (e.g., approaching versus leaving pericenter; Martínez-García et al. 2023).

However, robustly measuring velocity gradients of just a few km/s/kpc entails a significant challenge, depending on factors such as the spatial coverage, velocity uncertainties, and sample sizes. Therefore, they are worth revisiting as new datasets become available, in particular at a time when more precise determinations of the systemic proper motions allow for “perspective” gradients to be corrected (e.g., Kaplinghat & Strigari 2008; Walker et al. 2008).

Another main outcome of the spectroscopic campaigns has been the finding that the stellar component of several LG dwarf galaxies can be described as the superposition of at least two components with a different metallicity, spatial distribution, and kinematics (dubbed “chemo-kinematic components” or “multiple chemo-dynamic populations,” e.g., Tolstoy et al. 2004; Battaglia et al. 2006; Walker & Peñarrubia 2011; Fabrizio et al. 2016; Kordopatis et al. 2016; Pace et al. 2020). This has been possible through the combined analysis of line-of-sight (l.o.s.) velocities and metallicities (or indicators of relative metallicity) for large samples of individual stars. Trying to understand whether these are indeed separate populations can give information on formation scenarios. For example, multiple populations can form from early mergers (that would puff up a preexisting old and metal-poor stellar component), followed by later gas accretion, which would form the more metal-rich component (e.g., Benítez-Llambay et al. 2016). Or, possibly, these components could be due to star formation bursts, or prominent phases of star formation at given times, rather than a more continuous star formation activity that would favor the formation of metallicity gradients (Revaz & Jablonka 2018; Mercado et al. 2021). From a dynamical modeling perspective, the importance of multiple chemo-kinematic populations lies in the fact that they can be used to place tighter constraints on the DM halo density distribution (e.g., Battaglia et al. 2008a; Walker & Peñarrubia 2011; Amorisco & Evans 2012; Zhu et al. 2016; Strigari et al. 2017) by relieving the well-known mass-anisotropy-(stellar density) degeneracy (Binney & Mamon 1982).

The Sculptor dSph, found at an heliocentric distance of 83.9 ± 1.5 kpc (Martínez-Vázquez et al. 2015), was the first LG dwarf galaxy in which the presence of multiple chemokinematic components was unveiled (Tolstoy et al. 2004) and used as independent tracers of the galaxy’s gravitational potential (Battaglia et al. 2008a). The rationale for considering these as separate components came mainly from the presence of a well-defined blue- and red-horizontal branch in the color-magnitude diagram of this galaxy, which also exhibited a distinct spatial distribution. Additionally, Schwarzschild modeling of this galaxy naturally led to the identification of two peaks in the energy-angular momentum space, directly related to the two previously known populations (Breddels & Helmi 2014), giving support to the hypothesis of two distinct stellar components. Dynamical modeling using both Sculptor’s chemo-kinematic components individually, typically has resulted in a preference for a cored DM halo density distribution (Battaglia et al. 2008a; Walker & Peñarrubia 2011; Amorisco & Evans 2012; Zhu et al. 2016; Strigari et al. 2018), while single component modeling results are more diverse, finding both cores in some cases and cusps in others (Breddels et al. 2013; Read et al. 2019; Pascale et al. 2019; Kaplinghat et al. 2019; Pascale 2020).

Battaglia et al. (2008a) also found a l.o.s. velocity gradient of ${7.6}_{-2.2}^{+3.0}$ km s⁻¹ deg⁻¹ along the projected major axis. Nonetheless, the evidence for a significant velocity gradient in this galaxy has weakened, both in studies that have used Gaia data release 1 (DR1)+HST or Gaia early data release 3 (eDR3) systemic proper motions to correct the earlier spectroscopic samples for perspective velocity gradients (e.g., Massari et al. 2018; Martínez-García et al. 2021) as well as in new larger datasets with improved, Gaia-aided memberships (Tolstoy et al. 2023). Additionally, in the initial analysis of Tolstoy et al. (2023), the two previous chemokinematic components are no longer as easily identifiable as in earlier studies; while, on the other hand, the most distinct kinematics is detected for stars with [Fe/H] <−2.5, whose average velocity is shifted by a few km/s from that of the rest of the stars.

Our goal is to perform a comprehensive and quantitative chemo-dynamical analysis of the stellar populations in Sculptor, using the new dataset of Tolstoy et al. (2023) and simultaneously considering the information on the stars’s spatial location, l.o.s. velocity, and metallicity. We aim to untangle whether Sculptor is better described by a single population with a negative metallicity gradient or whether it is composed of distinct stellar components, and to revisit the presence of a possible v_los gradient. The results of this analysis will inform future dynamical modeling (Arroyo-Polonio et al., in prep).

In Sect. 2 we briefly describe the main characteristics of the dataset used. In Sect. 3 we characterize Sculptor’s stellar populations. In Sect. 4, we present the results of the search for 1.o.s. velocity gradients. In Sect. 5 we discuss the implications and possible explanations of our results in the context of all the findings in the literature. Finally, in Sect. 6 we summarize our work, draw the main conclusions, and discuss future prospects.

2 Data

In this work we make use of the v_los and metallicity ([Fe/H]) measurements of individual stars presented in Tolstoy et al. (2023). We refer the reader to this article for a detailed discussion of the characteristics of the dataset, the data reduction process, and the methodology to determine v_los and [Fe/H].

Here suffices to say that the spectroscopic dataset, on which the measurements are based, was acquired with the GIRAFFE spectrograph in Medusa mode, on the FLAMES instrument at the Very Large Telescope (VLT). The grating used was LR8, covering between 820.6 and 940.0 nm, with a spectral resolving power of 6500 (Pasquini et al. 2002). This wavelength region includes the near-IR Ca II triplet (CaT) lines, a well-known and well-tested indicator of the metallicity of red giant branch (RGB) stars also in composite stellar populations (Battaglia et al. 2008b; Starkenburg et al. 2010; Carrera et al. 2013; Vásquez et al. 2015).

The dataset consists of 67 independent observations of 44 pointings, from ten different observing programs acquired over the period from 2003 to 2018; these were homogeneously reduced and analyzed, using the most recent ESO pipeline, and then specific software developed over many years by M. Irwin. A zero-point calibration, based on a global shift applied to the spectra for each individual field, was applied to ensure that there were no velocity offsets between different exposures. Parallaxes and proper motions from Gaia Data Release 3 (DR3) (Gaia Collaboration 2021), as well as the FLAMES v_los determinations, were jointly used in order to effectively select stars members of Sculptor and discard contaminants.

This resulted in a sample of 1604 individual RGB stars likely members of Sculptor for which the spectrum signal-to-noise ratio (S/N) was large enough (>5) to yield precise v_los; among these, 1339 have also reliable metallicities. Since in our analysis we will use both the metallicity and v_los information, we will mainly focus on the second set, unless explicitly said otherwise. The spectra of the selected stars have a mean S/N of 53 with a minimum of 13. The velocities have a mean uncertainty of ±0.6 km s⁻¹, while the metallicities have mean errors around ±0.1 dex. In Appendix A we test the robustness of the velocity uncertainties provided in the dataset and propose corrections for them.

The large angular size of Sculptor causes the 3D systemic velocity of the galaxy to project differently over the lines-of-sight to different stars (Kaplinghat & Strigari 2008; Walker et al. 2008). This gives rise to a perspective velocity gradient that can be mistaken for rotation. Fortunately, one can correct the l.o.s. velocities for this effect as long as the galaxy’s systemic proper motion is known. To this aim, we adopted the formulae in the Appendix of Walker et al. (2008) and the systemic proper motion for Sculptor listed in Table 1. These values are already corrected for the zero-points offset measured with quasars located within 7 degrees from Sculptor (Battaglia et al. 2022). The perspective velocity gradient is visible in Fig. 1: it amounts to about 1 km s⁻¹ in the outermost region (around 1 degree from the center) at a position angle of 37°. From now on, we will work with the l.o.s. velocities corrected for the perspective velocity gradient.

3 Chemo-dynamical characterization

In Sect. 3.1, we present the methodology used for the identification of different stellar populations. In Sect. 3.2, we show the results of the chemo-dynamical characterization of Sculptor’s stellar populations, and analyze the performance of the method on a mock galaxy.

Fig. 1 Projection on the sky of the spatial distribution of Sculptor’s member stars, color-coded by the value of the l.o.s induced velocity because of the perspective gradient (see the color bar). The purple arrow shows the direction of the systemic proper motion of Sculptor. The black ellipse has the ellipticity and position angle in Table 1 and semi-major axis equal to twice the semi-major projected half-light radius.

3.1 Methodology

According to the analysis provided by Tolstoy et al. (2023), the existence of Sculptor’s two distinct chemo-dynamical populations is not as clear as before. Rather, Sculptor seems to be possibly well described by only one population with a negative metallicity gradient. However, since the metallicity, position, and l.o.s. velocity distributions of Sculptor were analyzed separately, here we revisit the modeling jointly analyzing these quantities to properly account for possible dependencies in Sculptor’s chemodynamical properties. In practice, we compare three different models: a single population model, where the stellar component is described by one surface number density profile and a constant l.o.s. velocity dispersion, but exhibits a metallicity gradient (1-pop); a two-population model, where the stellar component is described via the super-position of two components, allowing for differences in their spatial distribution, global l.o.s. velocity dispersion and metallicity distribution (2-pop); and finally, a three-population model (3-pop), which we used to search for statistical evidence of some other population or residual contamination from MW stars (although it is important to point out that in our case all the stars are astrometric and l.o.s. velocity members of Sculptor, hence if present, residual contamination should be minimal).

We characterize the chemo-dynamical properties of Sculptor’s stellar population broadly following the methodology presented in Walker & Peñarrubia (2011); Pace et al. (2020), with some modifications.

In Sect. 3.1.1, we quantify the selection function of the VLT/FLAMES dataset. In Sect. 3.1.2, we define the likelihood function used to compare different models to the observed data. In Sect. 3.1.3, we explain how the parameter space is explored and how our results from different models are compared. Finally, in Sect. 3.1.4 we explain the method used to compute the probabilities of membership to each population.

3.1.1 Selection function

The number of spectroscopically observed member stars is only a subset of the overall population of Sculptor’s stars. Therefore, we include a selection function describing the ratio between spectroscopically observed member stars and overall members to Sculptor. This ratio is not uniform across the galaxy nor as a function of the stars’ magnitudes, as a consequence of the placement of the fibers, of the magnitude limit and other observational conditions.

In order to correct any selection bias, we compare our dataset with the confirmed astrometric members of Tolstoy et al. (2023), derived using only Gaia DR3 data, considered complete down to a magnitude limit of 20 (Cantat-Gaudin et al. 2023).

Therefore, we define the selection function as: $$\omega (R,G)=\frac{d{N}_{obs}(R,G)}{d{N}_{mem}(R,G)},$$(1)

where dN_obs(R, G) and dN_mem(R, G) are the number of spectroscopically observed and astrometric member stars with semi-major axis radius¹ between R and R + dR and Gaia G-magnitude between G and G + dG, respectively.

Even though our dataset is discrete, instead of binning the data to build a discrete selection function, we adopt a smoothing Gaussian-Kernel to model it as a continuous function in the following way: $$\omega (R,G)=\frac{\sum _i^{N_{obs}}\exp -[\frac{{(R_i^{obs}-R)}^2}{k_R^2}+\frac{{(G_i^{obs}-G)}^2}{k_G^2}]}{\sum _j^{N_{mem}}\exp -[\frac{{(R_j^{mem}-R)}^2}{k_R^2}+\frac{{(G_j^{men}-G)}^2}{k_G^2}]},$$(2)

where the subscripts i and j range over all the spectroscopically observed member stars and all the astrometric members, respectively. $R_i^{obs}\left(G_i^{obs}\right)$ and $R_j^{mem}\left(G_j^{mem}\right)$ are the semi-major axis radii (magnitudes) of the spectroscopically observed and astrometric member i and j stars, respectively; $k_R^2$ and $k_G^2$ are the smoothing parameters for the semi-major axis radius and the magnitude, respectively. The two smoothing parameters must be chosen taking into account that the density of stars in each point of the (R, G) plane is different; for example, the maximum density in the number of stars is found close to the center and at G~19.

We adopt the following parametrizations of the softening parameters to describe our dataset²: $k_R^2=0.008+R/12+R^2/10$ where R is expressed in deg (so $k_R^2$ ranges from 0.5 arcmin in the inner parts to 16 arcmin in the outer part) and $k_G^2=0.25+|G-19|/6$ (from 0.25 mag in the region around G = 19 to 0.75 for the brighter magnitudes). As a comparison, panels a) and b) of Fig. 2 compare the selection function computed using eq. (2) along with the aforementioned parametrization for $k_R^2$ and $k_G^2$ (left panel), with the selection function computed binning the dataset on a 15×15 grid.

In order to simplify the terminology, we also define the integral of the selection function over all the magnitudes as: $$\mathrm{\Omega }(R)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\omega (R,G)dG.$$(3)

In panel c) of Fig. 2 we can see Ω(R) for the chosen softening parameters, compared with the observational data for 100 different arbitrary binning grids: the general good agreement can be visually appreciated here. One can notice that there is a small discrepancy between the data and the modeled selection function in the inner parts and this might warrant a different softening parameter for those regions; however, this is difficult to determine with the data at hand, since we want to model the 2D map G vs R, and we do not have many stars in the inner part for all the values of G.

Fig. 2 Selection function of the dataset. Panel a) shows the continuous selection function used in this work and defined in Eq. (2); panel b) shows the ratio between spectroscopically observed and astrometric member stars for a 15×15 gird in the (R, G) plane. Panel c) shows the selection function integrated over all magnitudes (Eq. (3)) as a solid line and the ratio between observed and member stars for 100 different grids of binning in R.

3.1.2 Likelihood

We work within the framework of Bayesian inference. Recalling Bayes theorem, we have that $$P(Y\mid \zeta ,M)=\frac{P(\zeta \mid Y,M)P(Y\mid M)}{P(\zeta \mid M)}.$$(4)

Above, ζ indicates the observables that we are going to fit, M the model that we are going to use to reproduce our observed dataset, and Y the parameters used to describe that model. P(ζ|Y, M) is the likelihood function (hereafter, ℒ), which represents the probability that our observational data ζ is described by the parameters Y using the model M. P(Y|M) is the prior probability of having a set of parameters Y favored in the model M. P(Y|ζ, M) is the posterior probability distribution, the probability of a set of parameters Y in the model M that describes our data. Finally, P(ζ|M) is the normalization factor.

Each star of our sample is characterized by a set of observables, {ζ_i} = {R_i,[Fe/H]_i, v_los,i} with i = 1,.....,N_obs, that is, R_i the semi-major axis radius, [Fe/H]_i the metallicity, v_los,i the l.o.s. velocity of the i star and N_obs the total number of stars in the sample, respectively; when dealing with the selection function, we will also consider the star’s G magnitude.

For a N_pop-population model, being N_pop the total number of stellar populations in the model, the likelihood can be expressed as³: $$\log (\mathcal{L})={\displaystyle \sum _{i=1}^{N_{obs}}}\log \left({\displaystyle \sum _{p=1}^{N_{pop}}}\frac{f_p\omega (R_i,G_i){\mathcal{L}}^p\left({\zeta }_i\right)}{\int \int \int \int \omega (R,G){\mathcal{L}}^p(\mathit{\zeta })d^3\zeta dG}\right).$$(5)

Above, the index p runs over the populations and the index i over all the stars in our sample; f_p is the fraction of stars that belong to the population p with respect to the total number of stars in our sample, so that ${\sum }_{p=1}^{Npop}{f}_p=1;\omega (R_i,G_i)$ is the value of the selection function corresponding to the R and G of the i star; and finally, ℒ^p (ζ_i) is the individual likelihood of the star i belonging to the population p. In the denominator d³ ζ indicates that we are integrating over R, [Fe/H] and v_los, and dG over the magnitude as well. As ℒ^p does not depend on the magnitude G, we can integrate the selection function over the magnitude and the denominator reduces to ∫ ∫ ∫ Ω(R) ℒ^p (ζ) d³ ζ.

The individual likelihood ℒ^p of a star in a certain population is: $${\mathcal{L}}^p\left({\zeta }_i\right)=P_r^p\left(R_i\right)P_{\mathcal{M}}^p([Fe/H{]}_i,R_i)P_{\mathcal{V}}^p\left(v_{\text{los },i}\right),$$(6)

where $P_r^p\left(R_i\right),P_{\mathcal{M}}^p([Fe/H{]}_i,R_i)$ and $P_{\mathcal{V}}^p\left(v_{\text{los },i}\right)$ are the semi-major axis radius, metallicity and velocity probability distribution functions respectively. Hereafter, we will drop the super index p from the notation.

For the single population model with a metallicity gradient, 1-pop, the metallicity distribution depends on R and only the v_los term is separable. If we integrate over it, the denominator of eq. (5) becomes ∫ ∫ Ω(R) P_R P_ℳ dR dℳ.

For the multiple-populations models, 2-pop and 3-pop, we assume that P_ℳ does not depend on R, then ℒ(ζ_i) is fully separable and the normalization term of the denominator in Eq. (5) can be directly integrated on [Fe/H] and v_los; only the radial term remains, and the denominator of Eq. (5) becomes ${\int }_0^{\mathrm{\infty }}\mathrm{\Omega }(R)P_{R}dR$.

The stellar density profile is assumed to be a spherically symmetric Plummer distribution (Plummer 1911) for simplicity. Then the probability distribution of projected radii is: $P_r\left(R_i\right)=\frac{2R_i/R_h^2}{{(1+R_i^2/R_h^2)}^2},$(7)

where R_h is the projected half-light radius.

For the metallicity distribution, we assume that it follows a Gaussian: $$P_{\mathcal{M}}([Fe/H{]}_i)=\frac{1}{\sqrt{2\pi ({\sigma }_{\mathcal{M}}^2+e_{i,\mathcal{M}}^2)}}\exp (-\frac{1}{2}\frac{{([Fe/H{]}_i-\mathcal{M})}^2}{{\sigma }_{\mathcal{M}}^2+e_{i,\mathcal{M}}^2}),$$(8)

where, for the multiple-populations model, ℳ is the mean metallicity, σ_ℳ is the standard deviation of the metallicity distribution and e_i,ℳ the errors in the i-th metallicity measurement. In the one-population model with the metallicity gradient, ℳ = ℳ₀ + m R, where ℳ₀ is the metallicity at the center of the galaxy and m is the metallicity gradient.

Finally, we also assume that the v_los follow a Gaussian distribution: $$P_{\mathcal{V}}\left(v_{los,i}\right)=\frac{1}{\sqrt{2\pi ({\sigma }_{\mathcal{V}}^2+e_{i,\mathcal{V}}^2)}}\exp (-\frac{1}{2}\frac{{(\mathcal{V}-v_{\text{los},i})}^2}{{\sigma }_{\mathcal{V}}^2+e_{i,\mathcal{V}}^2}),$$(9)

where 𝒱 is the mean l.o.s. velocity, σ_𝒱 is the standard deviation of the v_los distribution and e_i,𝒱 the errors in the l.o.s velocity measurements.

3.1.3 Model comparison and MCMC

The three models that we are exploring have a different number of free parameters: 1-pop has six free parameters {R_h, ℳ₀, m, σ_ℳ, 𝒱, σ_𝒱}; 2-pop has 11 free parameters, that is, a set of five parameters ({R_h, ℳ, σ_ℳ, 𝒱, σ_𝒱}) per population, and one extra parameter f which is the fraction of stars belonging to the population 1 over the total. Finally, 3-pop has 17 free parameters, that is, one set of five parameters per population, and two extra parameters f₁ and f₂, that are the fraction of stars belonging to population 1 and 2 over the total, respectively.

To explore the parameter space and sample from the posterior distribution, we use the EMCEE package (Foreman-Mackey et al. 2013), a python implementation of the Affine Invariant Markov Chain Monte Carlo (MCMC) ensemble sampler (Goodman & Weare 2010), further details about the MCMC runs are explained in Appendix B. However, since the likelihoods for different models are different, and they have a different number of parameters, we cannot directly use this method to choose the best model. For that, one can use the Akaike information criteria AIC = k − 2 log(ℒ) or the Bayesian information criteria BIC = k log(n) − 2 log (L), where k is the number of free parameters, n the number of data points that we are using to compare (1339 stars in our case), and ℒ the value of the median of the likelihoods for the posterior probability distribution. We use the second one, as it penalizes more strongly the extra free parameters when comparing numerous points, as in our case. A ΔBIC = BIC(A) − BIC(B) larger than 10 corresponds to very strong evidence of the A model being favored over the B one (Bauldry 2015). Finally, to characterize the parameters that fit better our data, we define the best-fitting parameters as the median and the 16th and 18th percentiles of the posterior probability distribution for each parameter.

3.1.4 Probability of membership

The method used allows us to attribute a probability of membership of each star to a given population, we define: $$L_i^p=\frac{f_p\omega (R_i,G_i){\mathcal{L}}^p\left({\zeta }_i\right)}{\int \int \int \mathrm{\Omega }(R){\mathcal{L}}^p(\zeta )d^3\zeta }.$$(10)

So that, for each star i, the probability of belonging to the population p is: $$P_i^p=\frac{L_i^p}{\sum _k^{N_{pop}}{L}_i^k}.$$(11)

In our Bayesian approach, we compute the probability of membership taking into account all the posterior probability distribution. We select 5000 models according to it and compute the probability of membership for each one of them. Then we use its median, scaled so that the sum of probabilities for each star of belonging to one of the populations is 1, as the final indicator for the probability of membership.

3.2 Results

In Sect. 3.2.1, we present the results of our methodology applied to the dataset from Tolstoy et al. (2023). As explained in Appendix A, the authors provided two datasets in which the l.o.s. velocities and corresponding uncertainties were derived using different methods. After comparing the velocity variation with the uncertainties, we found the uncertainties not to be well estimated, and therefore, we provided some corrections. We decided to use all sets of velocities in this work, the corrected and uncorrected versions; see Appendix A for the details. Since the results produced in all the cases are the same, here we quote only those from the same velocities and uncertainties set used in Tolstoy et al. (2023); the rest of them are shown in Appendix C. In Sect. 3.2.2, we apply the method to a mock galaxy to validate our methodology.

3.2.1 Sculptor

In Table 2, we report the best-fitting parameters for the three different models, along with the ΔBIC (see Sect. 3.1.3) with respect to 1-pop (BIC_1−pop − BIC_n−pop).

In 1-pop, we recover the results by Tolstoy et al. (2023), where they found a systemic velocity of 111.2 ± 0.2 km s⁻¹ and a metallicity gradient around −0.7 dex/deg. Those values are in excellent agreement with our results, systemic velocity of 111.2 ± 0.3 km s⁻¹ and a metallicity gradient of 0.69 ± 0.06 dex/deg.

With a ΔBIC = −40, there is very strong evidence that 2-pop provides a much better description of the data than 1-pop. In this case, the distribution of metallicities are centered at $ℳ=-{1.44}_{-0.03}^{+0.03}$ for the metal-rich (hereafter, MR) and $-{2.02}_{-0.02}^{+0.02}$ for the metal-poor (MP) components, with a scatter of ${\sigma }_{ℳ}={0.27}_{-0.02}^{+0.02}$ and ${\sigma }_{ℳ}={0.34}_{-0.01}^{+0.01}$ respectively. The metallicity distribution function of both populations overlaps within 1σ, so they are not well differentiated when considering only this observable. The MR component has a lower l.o.s. velocity dispersion than the MP one, ${\sigma }_{𝒱}={6.4}_{-0.4}^{+0.3}$ km s⁻¹ and ${11.1}_{-0.3}^{+0.3}$ km s⁻¹, respectively, and it is more spatially concentrated, the 2D semi-major axis half-light radii are ${0.129}_{-0.008}^{+0.008}$ deg and ${0.27}_{-0.01}^{+0.01}$ deg respectively.

In terms of kinematic properties, we recover the results of Tolstoy et al. (2004), Battaglia et al. (2008a) and Walker & Peñarrubia (2011) always within 1σ, but with considerably lower errors thanks to the quality of the dataset. For the metallicity distribution we cannot directly compare with the literature as in Tolstoy et al. (2004) and Battaglia et al. (2008a) the authors made a hard-cut in the metallicity to distinguish between populations and Walker & Peñarrubia (2011) used the Magnesium-index rather than [Fe/H]. Regarding the extension in the sky of both components, Walker & Peñarrubia (2011) found $R_h={0.113}_{-0.006}^{+0.007}$ deg for the MR and $R_h={0.24}_{-0.02}^{+0.02}$ deg for the MP using circular projected radii, that are systematically smaller than the semi-major axis projected radius. Therefore, the differences in the results are attributed to this effect. Moreover, our R_h are within 1σ compared with the ones from Battaglia et al. (2008a), obtained using the semi-major axis projected radius.

Finally, in 3-pop, the first two components have very similar model parameters to those already found in the 2-pop case for the MR and the MP. However, the presence of a third population (pop 3, hereafeter) with a few stars is inferred⁴. This population is very metal-poor with a metallicity around $ℳ=-{2.90}_{-0.2}^{+0.3}$; it is very extended and displaced with respect to Sculptor’ center; and it displays a shift of 15 km s⁻¹ in the mean v_los (${125.5}_{-2.6}^{+2.4}$ km s⁻¹) with respect to the systemic one in Sculptor (${111.2}_{-0.3}^{+0.3}$ km s⁻¹). This is in line with the results from Tolstoy et al. (2023), who found an increase in the mean v_los as a function of radius. Other works have also found extended stellar components in this dwarf galaxy (Jensen et al. 2023) by analyzing the number surface density profiles derived from Gaia DR3 probable members. 3-pop has a ΔBIC = −69 and ΔBIC = −29 with respect to 1-pop and 2-pop, respectively. Therefore, there is statistically significant evidence to choose 3-pop as the model that describes best the data.

In Fig. 3 we show the projection of the stars in the v_los vs [Fe/H] and R vs v_los planes, color-coded according to the population they have the higher probability of belonging, according to Eq. (11), as well as the overall observed and predicted metallicity and velocity distributions. In this figure, the difference in the velocity dispersion between the MR and MP populations can be seen clearly by looking at the different spreads in the l.o.s. velocity axis. It also shows that pop 3 is more metal-poor and has a higher v_los compared with the rest of the galaxy. Also, this component is not concentrated in the center but rather uniform in R. In the panel a) of Fig. 3 we show the observed metallicity distribution and the one predicted by the best-fit 3 population model. We can see that we are describing the two known populations, but also the extended metal-poor tail of the distribution with the third component. Here, we are showing the underlying metallicity distribution obtained with the models taking into account the selection function (dashed lines and bands), and comparing it to the observed metallicity (solid histogram). We experimented with running the code without taking into account the selection function, to get an idea of how much the underlying and observed metallicity distribution differ. The result remained the same for the metallicity distribution and only the R_h of the different populations changed. Therefore, we consider the one shown here to be a sufficiently good approximation of the underlined one.

In Fig. 4 we show the projection of the Sculptor’s stars in the sky, with the stars that belong to a specific population are color-coded by velocities. For this figure, we assigned each star to the population that yields the highest probability of membership. We can clearly see again the MR population with low velocity dispersion in the center of Sculptor and the MP population more extended and with higher velocity dispersion (left and middle panels, respectively. The third population is located in the outskirts of Sculptor, it is very extended, and its center is displaced compared to the one of Sculptor. In this figure, it is clear again that the systemic velocity of this population is shifted.

Finally, we computed the stellar surface number density and l.o.s. velocity dispersion profiles for both populations according to the membership criteria established in the previous paragraph (Fig. 5). Here, stars belonging to the third population were excluded due to their distinct characteristics, and off-centered spatial distribution. The velocity dispersion, and corresponding uncertainties, were determined by fitting a Gaussian distribution to the l.o.s. velocities of the stars in each bin and accounting for velocity uncertainties (see Sect. 5.1 of Taibi et al. 2018, for details). The number density profile has been computed by counting the number of stars in each bin, and we note that these are not the “true” underlying profiles, but rather the number of observed MR and MP stars. The selection function corrected ones are shown with bands, the main difference is that it rises the outer slope of the MP component, and it also increases slightly the density in the inner part for both components. Here we can see that the transition from the MR dominance to the MP occurs at around 0.15 deg.

For the velocity dispersion profile, the bin size was chosen as to maintain the number of stars per bin around 80. The velocity dispersion profiles of both populations differ of ~4 km s⁻¹ across the whole radial extend. They are slightly flatter than those obtained by Battaglia et al. (2008a). Compared with the σ_los profiles computed by Strigari et al. (2017) using the data from Walker & Peñarrubia (2011), ours show more difference in velocity dispersion between both populations, mostly in the inner part. This could be because the populations are better differentiated in [Fe/H] than in the magnesium-index used in that work These changes in the shape of the velocity dispersion profiles raise the urge to perform new dynamical modeling to infer how the dynamical mass and the dark-to-luminous mass ratio change with this new dataset.

Fig. 3 Chemo-dynamical properties of Sculptor. (a) / (d) panels: Normalized metallicity and l.o.s. velocity distribution of the spectroscopically observed member stars of Sculptor fitted by the 3-pop model. The dashed-lines indicate the predicted distribution of our models: red for the MR component, blue for the MP, green for the third population, and black for the sum of the three components as indicated in the label. Bands indicate the 1σ confidence intervals. (b) / (c) panels: Projection of the stars in the vlos vs [Fe/H] / R vs vlos planes. In each panel, the stars are color-coded by the population to which they have the highest probability of belonging, red for the MR, blue for the MP and green for the third population as indicated in the label. The systemic velocity of Sculptor is indicated with straight dashed lines.

Fig. 4 Projection on the sky of all stars in the LR8 sample (gray). In each panel we assign stars to each population according to the criteria explained in the text (from left to right population 1,2 and 3) and we color-code them by vlos, as indicated in the color bar.

Fig. 5 Upper panel: projected surface number density profile for the MP and MR components, obtained directly by the observations and corrected by the selection function, as indicated in the legend. Lower panel: vlos dispersion profile for the MP and MR components, as indicated in the legend.

3.2.2 Reliability test

In order to test the performance of our methodology, we also analyzed a mock galaxy. We chose to use one from the Gaia challenge wiki⁵ (Walker & Peñarrubia 2011), so that we analyze a galaxy generated by distribution functions not related to the distributions that we used for the likelihood. This mock is a N-body realization extracted from a distribution function assuming two distinct stellar populations produced by Hernquist density profiles in a Hernquist dark matter halo. The orbits of the stars have been integrated 100 crossing times to ensure stability. For more details, we refer the reader to the documentation in the wiki, the specific file used is “SplitCompCore.data”. We generate a random sample of 1339 stars of that mock galaxy. For the velocities, we use the same errors we have in our dataset for each specific stars by adding Gaussian errors to the velocities. The [Fe/H] is not provided in the mock, but the reduced magnesium index (Mg-index), and populations are not as well differentiated in Mg-index as they are in [Fe/H] for Sculptor. That is why we decided to increase the distance between the Mg-index distributions of both populations by 0.3 dex so that they are more alike to the [Fe/H] distinction we have in Sculptor. The errors are the same as the ones provided in the mock. Moreover, we added an artificial component of 25 stars resembling the third population we find in Sculptor. The parameters of the stellar populations in the mock are listed in Table 3. The metallicities and velocities have been computed by selecting all the stars in the mock belonging to populations 1 and 2, and computing their mean value and dispersion deconvolving from errors by using Gaussian distributions.

We fit the mock using the same three different models used in the previous section. We find that 3-pop is the one that describes better the data, showing a ΔBIC of −322 and −122 with respect to 1-pop (BIC_1−pop − BIC_3−pop) and 2-pop (BIC_2−pop − BIC_3−pop), respectively. In Table 3, we compare the mock parameters with thones inferred from our analysis. We can clearly see that all the parameters are within 1σ but the systemic velocity of the MR component and the velocity dispersion of the MP component which are within 2σ, this serves to validate our methodology.

Now we use this result to analyze how well the membership probability method performs. First, we test the same criteria used for assigning the membership in Fig. 4, that is, we assign a star to the population it has the highest probability of belonging. In this case, the purity⁶ of the sample is 87, 95, 100% for the stars assigned to populations 1, 2, 3, respectively. If, in instead of choosing the population with higher probability to belong, we put a threshold of a membership probability >80%, the purity of the samples goes up to 95, 98, 100% for populations 1, 2 and 3, respectively. But at the cost of not having all stars classified into a certain population.

We performed a final reliability test, as we want to be sure that our method is not biased toward finding new populations. We run the same code on the mock galaxy but without the third population of 25 stars in the dataset. In this case, the favored model is 2-pop, showing a ΔBIC of -106 and -57 with respect to 1-pop (BIC_1−pop − BIC_2−pop) and 3-pop (BIC_3−pop − BIC_2−pop), respectively.

4 Velocity gradient

Here we analyze the sample of 1339 probable member stars to look for signals of a velocity gradient in Sculptor. Apart from perspective effects, which we have corrected for (see Sect. 2), a velocity gradient could be due to internal rotation or tidal effects⁷. However, we point out that we consider the latter hypothesis unlikely, as by re-analysis the Iorio et al. (2019) N-body simulations, tailored to reproduce a Sculptor-like galaxy on a Sculptor-like Gaia motivated orbits, show that we should not be detecting a tidally induced velocity gradient in the region we are studying. Therefore, in the following, we will use the terms “rotation” and “velocity gradient” interchangeably.

To search for possible rotation in Sculptor, we use the same methodology as the one presented in Taibi et al. (2018), which we summarize here. We look for evidence of rotation by comparing two different models: a radial velocity gradient model, in which apart from the velocity dispersion in the galaxy, there is a v_los gradient signal. It is parametrized as v_los = v_sys + nR_i cos(θ_i − θ), where v_sys is the systemic velocity of the galaxy, n the magnitude of the gradient in km s⁻¹ deg⁻¹, R_i and θ_i are the distance from center of the galaxy⁸ and the position angle of the i star, and θ is the position angle of the maximum rotation gradient, perpendicular with respect to the rotation axis; and a dispersion-only model, in which we assume that the galaxy is fully pressure-supported. We use again a Bayesian approach to explore the posterior probability distribution, we assume Gaussian velocity distributions centered in v_los and with a velocity dispersion σ_los to build our likelihoods, so in total there are 2 free parameters for the dispersion-only model (v_sys, σ_los) and 4 for the rotational model (v_sys, σ_los, n, θ). For this analysis, we can compute the Bayesian evidence as the space parameter is explored using nested sampling. Therefore, in order to determine which one better describes the data, we use the Bayes factor B₁₂ = Z₁/Z₂, the ratio between the Bayesian evidences of two different models, 1 and 2. The Bayes factor inherently penalizes the introduction of additional free parameters in the rotational models. To compute the Bayesian evidences, we use the Multinest code (Feroz et al. 2009). The significance of one model with respect to another can be based on Jeffrey’s scale, computing the natural logarithm of the Bayes factor: values from (0–1), (1–2.5), (2.5–5), (5+) correspond to inconclusive, weak, moderate and strong evidence favoring the first model over the second one (Kass & Raftery 1995).

In the upper part of Table 4, we report the best fit parameters of the two models. The Bayes factor between the dispersion-only model and the radial model is 0.6 favoring the radial model. However, according to Jeffry’s scale, this is still inconclusive evidence of the model with rotation being favored over the dispersion-only one. In Battaglia et al. (2008a) the authors found a velocity gradient along the major axis analyzing a subset of this FLAMES dataset. The angles of the rotation axes are different by 56 degrees, and the velocity gradient they found was $-{7.6}_{-2.2}^{+3.0}$ km s⁻¹ deg⁻¹, in agreement at 2σ level with ours. It is important to note that the authors do not correct by the perspective induced velocity gradient (see Fig. 1) but for the motion of the Local Standard of Rest around the Milky Way center. In Martínez-García et al. (2023) the authors also looked for l.o.s. velocity gradients in Sculptor using the dataset from Walker et al. (2009), finding, however, a gradient consistent with 0 km s⁻¹ deg⁻¹ within 1σ. In this case the differences in the gradient are likely to be caused by the different datasets analyzed.

We also look for velocity gradients independently in the MR and MP populations. We find similar results as for the general analysis: gradients of −4.4 ± 2 km s⁻¹ kpc⁻¹ and −3.0 ± 3 km s⁻¹ kpc⁻¹ with θ = 32 degrees ± 20 and θ = 25 ± 40 for the MP and MR components, respectively. However, according to the Bayes factor, the dispersion-only model is favored with inconclusive evidence for the MP and weak evidence for the MR, showing a B₁₂ of −0.9 and −2, respectively. So, there is no evidence of rotation in any of the populations independently, neither. This is in contrast with the results from Zhu et al. (2016), who found rotation in the MR and the MP independently around different axes. They report amplitudes of the rotational signals of 1.1 ± 0.1 km s⁻¹ and 0.9 ± 0.1 km s⁻¹ with the maximum gradient at position angles around 30 deg (changing the sign of the amplitude to match our sign) and 140 deg, for the MR and MP components respectively. The amplitudes are not directly comparable, as they assume radially constant rotation in their models, but the angles of rotation are within 1σ for the MR component and completely off for the MP one. However, these signals could be artificially introduced by the rotational models they assume to assign the membership of the stars. This is why they test models without rotation in their Appendix B, in this case, only the MR component shows rotation, which is the signal similar to the one we found. Moreover, they do not compare the results with a simpler dispersion-only model to evaluate the evidence of the signals. It is in this test of comparing models where we find that the rotation models are not a better description of the observational data for Sculptor compared to dispersion-only models.

Fig. 6, shows, in the top panel, the l.o.s. velocity field map computed using Voronoi bins with, overplotted with the rotational axis. The bottom panel shows, instead, the binned rotational gradient compared to that of the model. In the upper panel we can see a ring of low velocity stars at around 15–20 arcmins from the center, more prominent in the upper part. However, given the lack of statistics, we cannot confirm this is a physical structure. In order to see if this is an artifact of the binning procedure or not, we show velocity maps computed following different approaches in Appendix D. Also, in the velocity map of Fig. 6, there is a high-velocity region in the outer galaxy parts, more intense in the right-bottom part. This is the region where the stars belonging to the third population are located (see Fig. 4). To test whether those high-velocity stars are responsible for the rotation found, we repeat the same analysis, but removing the stars with larger probability of belonging to population 3 than to the MR or the MP ones. The results are listed in the lower part of Table 4. ⟨v_los⟩ and σ_los are very similar to the model with all the stars and within 1σ. Second, the velocity gradient found is now slightly larger but with larger uncertainties too. Third and more important, the Bayes factor between the dispersion-only model and the rotational model drops by 1.5 points to −0.9 favoring the dispersion-only model, meaning that there is even less evidence of a velocity gradient. In general, removing or not the third population, our results are in agreement with the latest results in the literature, in which they find no evidence of a v_los gradient in Sculptor (Martínez-García et al. 2023; Tolstoy et al. 2023).

Fig. 6 L.o.s. velocity gradient of Sculptor. Upper panel: L.o.s. velocity map of Sculptor centered in its systemic velocity and generated with Voronoi bins; it is color-coded by velocities as indicated in the color bar. The black dots indicate the position of the stars in our sample; the purple arrow indicates the systemic proper motion of the galaxy; finally, the green line is the rotation axis, with the bands indicating the 1σ error range. Lower panel: velocity profile of Sculptor along the direction perpendicular to the rotation axis (maximum velocity gradient). The dots indicate the value of the mean velocity with errors in certain bins at different distances to the rotation axis, in black using all stars and in gray removing those from pop 3; the solid line shows the rotational gradient of 4.0 km s−1 deg−1, with the bands showing the 1σ error range; finally, the dashed line indicates the systemic vlos of Sculptor.

5 Discussion

In this section, we discuss the results in the context of the dSphs evolution. In Sect. 5.1, we consider the different implications of multiple populations in the context of Sculptor SFH. In Sect. 5.2, we search for the possible origin of the new third population found in Sculptor. Finally, in Sect. 5.3 we discuss whether there is enough evidence of a v_los velocity gradient and the implications it has in the evolution history of Sculptor.

5.1 Multiple populations in the context of Sculptor star formation history

We found that a multiple populations model describes Sculptor better than a single population model with a metallicity gradient. Then there are at least two main distinct stellar components in Sculptor, a MP extended one and a MR one more concentrated and less dynamically hot. There are two different scenarios that can lead to this system: first, the metal-rich population being formed internally after the metal-poor one, or second, by external causes such as gas infall or merger events. In the former scenario, after a strong initial star-forming epoch giving rise to the MP component, the star formation was partially quenched and the gas enriched. Then, the star formation started again in a more metal-rich environment, forming the MR component. This scenario has been shown to successfully reproduce the metallicity distribution of the Sculptor dSph, according to hydrodynamical simulations of isolated dwarf galaxies (Marcolini et al. 2008; Revaz et al. 2009). The overlap in metallicity occurs as most of the gas is retained, and the enrichment can be inhomogeneous. This process of self-enrichment would take some Gyr to occur for Sculptor (Babusiaux et al. 2005; Revaz et al. 2009) and would favor the extended SFHs proposed in some works in the literature (de Boer et al. 2012). The second scenario is that in-falling gas or an early merger with another dwarf galaxy produced a second burst of star formation. In this case, the process could happen faster, being able to explain the short SFHs proposed in some other works in the literature (Bettinelli et al. 2019; de los Reyes et al. 2022). In Benítez-Llambay et al. (2016) the authors analyze simulations of dSphs with two epochs of star formation, with the second one caused by interactions with another dSph. They are able to reproduce the observed characteristics of the stellar populations of dSphs like Sculptor. However, they choose galaxies with two distinct star-forming epochs separated by at least 3 Gyr in time. Given the latest reported short SFHs of Sculptor, simulations with closer star-forming epochs in time are needed to account for the possibility of short SFHs.

One could argue that any of those scenarios would have left some signatures imprinted in elemental abundances. Sculptor is known to display a “knee” in alpha elements at a metallicity [Fe/H] ~ −1.8 (Hill et al. 2019; de los Reyes et al. 2022); this is where we start to see a significant fraction of MR members, as shown in panel c of Fig. 3. This supports the idea that the MP component formed quickly, therefore presenting a constant abundance of α-elements. Subsequently, the MR population formed in an environment progressively enriched by type Ia supernovae. This could result from any of the two scenarios we presented above.

One could hypothesize that, since they are two well differentiated populations, there should be a change in the [α / Fe] scatter between the region where both populations coexist (−2 < [Fe/H] < −1) and the region dominated by the MR population ([Fe/H] > −1). We computed the scatter in [Mg/Fe] at different bins of [Fe/H] for the datasets of Hill et al. (2019) and de los Reyes et al. (2022). However, we can not draw any conclusion because of the large uncertainties. In the future, thanks to surveys like WEAVE and 4MOST we expect to increase the statistics, so that we can test this scenario.

5.2 Possible origins for the third population

In the following, we examine possible origins of the third population found in Sect. 3.2.1, with a very metal-poor mean metallicity, mean velocity shifted by 15 km s⁻¹ with respect to the other two populations and a position offset with respect to Sculptor’s center in the sky.

5.2.1 Contamination from disk-halo stars

We consider the possibility that the third, very metal-poor population could be due to contamination from MW stars, belonging to the disks or to the halo. There are several reasons to discard this possibility. If the third population is made by contaminant stars, we should find more of them with velocities lower than the systemic velocity of Sculptor. Also, if they were contaminants we would expect their spatial distribution in the sky to be uniform across the area sampled. Moreover, if we compare the characteristics of these stars and the ones that have been cataloged as contaminants by Tolstoy et al. (2023) they do not follow the same trends in [Fe/H], v_los and R. Furthermore, according to the Besançon model (Robin et al. 2003), not a single star with that v_los and [Fe/H] is expected in the direction of Sculptor that surpasses the probability membership criteria in terms of proper motions and parallaxes used to select the possible members. Finally, we have used the methodology presented by Battaglia et al. (2022) to assign probabilities of membership based on Gaia photometry and astrometry, position on the sky and augmented by the use of l.o.s. velocities, and the stars of this component appear as high probability members as well. We also cross-matched with a recent catalog of Sculptor members from Jensen et al. (2023), and all the stars in common from population 3 appear to be members in their selection too.

5.2.2 Contamination from streams

Approximately a hundred streams have been identified in the vicinity of the MW (Mateu 2023). Distributed across all the sky, these streams can reach very low metallicity and exhibit distinctive velocities, suggesting a potential explanation for the characteristics observed in the third population. There are two known streams of stars that pass close enough to Sculptor to be considered as possible sources of contamination. First, the Cetus-Palca stream (Newberg et al. 2009; Thomas & Battaglia 2022; Yuan et al. 2022). However, even though the position in the sky coincides, in the region around Sculptor, the v_los of the stream is around 40 km s⁻¹ and the proper motion is in the same direction but considerably larger than that of Sculptor. Thus, we discard Cetus-palca as a possible source of contamination. Second, there is the Murrumbidgee-G17 stream (Grillmair 2017). In this case the v_los coincides and its proper motion is roughly comparable to that of Sculptor, but it is so narrow in the sky that it could be responsible for at most 2 stars in the third component.

5.2.3 Tides

Another potential explanation is tidal effects due to the pericentric passages of Sculptor around the MW.

We analyzed the N-body simulations by Iorio et al. (2019) as these were tailored to reproduce many of Sculptor’s properties (3D position, 3D motion, surface number density profile, l.o.s. velocity dispersion profile) and in a set-up that would maximize the effects of tidal disturbances, in terms of MW mass, DM density profile for the simulated dwarf galaxy and Gyr of evolution.

The analysis suggests that tidal effects are an unlikely explanation for the origin of this population. First, the streaming motions caused by tides should appear at a completely different angle, that is, following the direction of Sculptor’s systemic proper motion. Second, tidal effects would shift the velocity just a couple of km s⁻¹ off the systemic velocity, at the location of the stars of the third population. Finally, the effect should be symmetric, that is, we should find a similar group of stars in the N-E direction at the same distance from the center and with a shift similar to the one of population 3, but mirrored to negative velocities.

5.2.4 Recent minor merger

The final physical explanation we propose is a recent minor merger. This scenario could explain the shift in position and in the mean v_los of this population. Some candidates for the minor merger, like ultra-faint dwarfs, can easily reach the metallicity of the third population (Simon 2019). On the other hand, the metallicities of pop 3 do not appear consistent with those of a stellar cluster, as the spread in metallicity of the most metal-poor clusters (Beasley et al. 2019; Martin et al. 2022) is not large enough to explain the metallicity dispersion we see in here.

In order to more directly compare the abundances of the most metal-poor stars of Sculptor and those of UFDs, and see whether pop 3 showed some specific signature, we performed a crossmatch of all our very metal-poor stars with medium-high resolution spectroscopic observations from the literature (Tafelmeyer et al. 2010; Starkenburg et al. 2013; Simon et al. 2015; Jablonka et al. 2015; Skúladóttir et al. 2024), and we compared their abundances with those in UFDs (Frebel et al. 2014; Roederer et al. 2016; Marshall et al. 2019; Ji et al. 2019; Chiti et al. 2018, 2023). We found literature elemental abundance measurements for 10 stars belonging to the pop 3 and 5 to the MP component. They fall within the scatter and trends for UFD stars at the same metallicity, as we show in Fig. 7. Therefore, at present, the comparison of elemental abundances does not exclude this scenario; although it should be emphasized that the number statistics is small and scatter exists between the abundances in various UFDs. The 4MOST/4DWARFS data, with the spectra calibrated and homogenized in the same way, will be transformative by reducing the scatter in the abundance measures and increasing the statistics (Skúladóttir et al. 2023).

Fig. 7 Chemical abundances of stars belonging to 3-pop and UFDs. From top to bottom, we present the [Sr/Fe], [Ba/Fe], and [Mg/Fe] abundances as a function of [Fe/H]. The green points represent pop 3 stars with available abundances from the literature, while the gray points indicate abundances from UFDs, as indicated in the label. Downward arrows indicate upper limits. Iron abundances are from the studies in the literature.

5.2.5 Limitations of the method

Finally, the limitations of our assumptions could be artificially producing this component. We assume that both the velocity and metallicity distributions have Gaussian shapes, however, it is common to see metal-poor tails in the metallicity distribution of dwarfs (e.g., Kirby et al. 2011; Leaman 2012; Taibi et al. 2022). Moreover, the models assume constant l.o.s. velocity dispersion profiles, but in Fig. 5 it can be clearly seen that the velocity dispersion profiles drop at least 2 km s⁻¹ for both, the MR and the MP populations through the radial extent. So, even though these are general assumptions and they work for distinguishing the main properties of large populations, they might not be general enough to properly characterize and differentiate such a small component as the one found here. We will address this problem by using general distribution functions to characterize the populations in Arroyo-Polonio et al. (in prep). Nonetheless, we point out that it is not a systemic difference in a single parameter what is characterizing this population. It is the combination of systematic differences in metallicities, l.o.s. velocities, and positions in the sky, what constrains statistically the existence of this component.

5.3 Velocity gradient

Even though it is true that we find a model with a 3σ detection for a velocity gradient, it does not have any statistical evidence if we compare it with a simpler dispersion-only model. Moreover, this evidence fades even more if we remove the third population of stars, whose origin is unclear. Therefore, we find no evidence of a detection of a l.o.s. velocity gradient in Sculptor.

The non-detection of a velocity gradient in Sculptor, up to this level, is also important to learn about the evolution of this galaxy. Recent literature reports that Sculptor recently passed its pericenter at around 47.7 kpc half gigayear ago (Battaglia et al. 2022)⁹. Taking into account both, the MW and the LMC to integrate its orbit, it would be its first pericenter passage, considering only the MW, it would be its second one (Battaglia et al. 2022). If we believe the scenario in which dSphs have become pressure-supported because of the interaction with the MW (Mayer et al. 2001), our results indicate that these effects have been strong enough to remove almost all the rotation in Sculptor in a single passage through the pericenter, or two if we do not consider the model with the LMC. However, the simulations from Kazantzidis et al. (2011) indicate that, for the pericentric distance of Sculptor, it is unlikely that one or two passages are enough to completely shut down the rotation of a dwarf galaxy. Moreover, other works in which they analyze the rotational state of more dwarf galaxies did not find clear evidence of a relation between the pericentric distance and the rotation support (Wheeler et al. 2017; Martínez-García et al. 2021). Notably, even in the core Sagittarius, a highly perturbed dSph suffering the stripping of the MW, substantial rotation has been found (del Pino et al. 2021), suggesting that tidal stripping alone may not be able to entirely eliminate the rotational motion.

6 Conclusions

In this work, we have reanalyzed the internal chemo-kinematic properties of the stellar population of the Sculptor dSph by using the new homogeneous VLT/FLAMES dataset from Tolstoy et al. (2023), with an unprecedented sample size, purity in the selection of member stars, and velocity and metallicity precision. We find that a model with two stellar populations with different half-light radii, mean metallicity, and l.o.s. velocity dispersion describes Sculptor much better than one with a single component with a negative metallicity gradient toward the outskirts. Our best-fitting models, however, also include a third component that accounts for few very metal-poor stars, which have a clearly shifted v_los with respect to the rest of the galaxy and a spatial distribution slightly displaced from the center. With this new dataset, the shape of the l.o.s. velocity dispersion for the MR and MP components that have been used in the literature for dynamical modeling of this galaxy have seen some mild variations. We have tested the robustness of the methodology used by applying it to mock data, and we find that it can correctly recover the parameters for the populations given in input with a high accuracy and that it can distinguish between 2-pop and 3-pop systems.

In the 3-pop model, the MR component has a mean metallicity [Fe/H] = $-{1.41}_{-0.03}^{+0.03}$, half-light radius R_h = ${0.126}_{-0.007}^{+0.008}$ deg, and a velocity dispersion σ_los = ${6.5}_{-0.3}^{+0.3}$ km s⁻¹. The MP component has a mean metallicity [Fe/H] = $-{2.00}_{-0.02}^{+0.02}$, halflight radius of R_h = ${0.26}_{-0.01}^{+0.01}$ deg, and a velocity dispersion σ_los = ${10.8}_{-0.3}^{+0.3}$ km s⁻¹. The third component has a metallicity of [Fe/H] = $-{2.9}_{-0.2}^{+0.3}$ and a mean velocity v_los = ${125.5}_{-2.6}^{+2.4}$ km s⁻¹.

As for the group of very metal-poor stars classified in the third population, it does not seem likely that we can explain it as contamination by known MW streams, or field disk or halo MW stars, and it is unlikely to be explained by tidal effects as well. Therefore, the only physical explanation we find is a recent minor merger, in which the outer stars of the galaxy are not dynamically thermalized yet. Still, around 20 stars are too few to be considered as strong evidence for its presence and it cannot be excluded that limitations in the modeling assumptions are driving the detection of this component. In order to confirm the existence of this third component, more spectroscopic followups on the outskirts of Sculptor will be needed. If, with more data, we keep finding more very metal-poor stars with the offset in velocity, we could characterize this new population more accurately.

We also revised the issue of whether Sculptor exhibits a 1.o.s. velocity gradient, taking advantage of the possibility, afforded by the knowledge of Sculptor‘s proper motion from Gaia eDR3 data, of correcting for spurious velocity gradients first due to perspective effects. We found inconclusive evidence of a l.o.s. velocity gradient of $-{4.0}_{-1.5}^{+1.5}$ km s⁻¹ deg⁻¹, after correcting for the perspective gradient; this is statistically indistinguishable from a model without a l.o.s. velocity gradient. This holds both when including or when excluding stars from the third component with a shifted l.o.s. velocity distribution. Therefore, this new dataset does not suggest the presence of statistically significant rotation in this galaxy.

Acknowledgements

J. M. Arroyo acknowledges support from the Agencia Estatal de Investigación del Ministerio de Ciencia en Innovación (AEI-MICIN) and the European Social Fund (ESF+) under grant PRE2021-100638. J. M. Arroyo, G. Battaglia, G. Thomas acknowledge support from the Agencia Estatal de Investigación del Ministerio de Ciencia en Innovación (AEI-MICIN) and the European Regional Development Fund (ERDF) under grant number AYA2017-89076-P, the AEI under grant number CEX2019-000920-S and the AEIMICIN under grant number PID2020-118778GB-I00/10.13039/501100011033. The authors acknowledge the referee for the constructive report, which enhanced the quality of the discussion in the manuscript.

Appendix A Analysis of the velocity uncertanties

The set of l.o.s. velocities on which the analysis in Tolstoy et al. (2023) is based were obtained via a maximum-likelihood method (MLM) where 3 parameters (v_los, Δ v_los, and the best-fit Gaussian width for blurring the delta function of the CaT template) are explored to obtain the most likely combination. However, the authors also determined an additional set of l.o.s. velocities through a cross-correlation (CC) with a template with a zero continuum and three Gaussians mimicking the CaT lines in absorption with the uncertainties being derived from the scatter of the velocities returned by fitting a Gaussian to each of the 3 CaT lines separately as in Battaglia et al. (2008b).

Since these methodologies produce different uncertainties, this prompted us to test whether the choice of one set over the other could have an impact on our analysis.

From the Tolstoy et al. (2023) sample, Arroyo-Polonio et al. (2023) had already analyzed the l.o.s. velocity measurements from spectra of individual exposures with high S/N and for stars with >8 exposures over a range of years, for the purpose of monitoring velocity variations. For that subset of measurements, the authors found that the uncertainties for the MLM method were underestimated while those from the CC method were slightly overestimated. In that case, the observations compared were taken over several months to years in order to characterize the motion of the binaries. Thereby introducing some velocity variation, even though it may be small, it is still present. Here we perform a similar analysis, but comparing pairs of observations of the same stars taken on the same night and focusing on a larger regime of S/N, to cover the overall range of the full Tolstoy et al. (2023) sample (180 stars with spectra with 13 < S/N < 76). Even though we are using observations obtained on the same night, there is still jitter introduced by stellar granulation that can produce velocity variation on shorter time scales (e.g., Yu et al. 2018). However, since we cannot model this velocity variation due to the randomness of the periodicity of these variation, this is essentially included within the uncertainties.

We compute the velocity difference normalized by the velocity uncertainties: $${\alpha }_i=\frac{v_1-v_2}{\sqrt{{\left(\mathrm{\Delta }{v}_1\right)}^2+{\left(\mathrm{\Delta }{v}_2\right)}^2}},$$(A.1)

where v₁ and v₂ are the l.o.s. velocities of a given star i in the two observations 1 and 2 in the same night, and Δ v₁ and Δ v₂ their respectively uncertainties. If the statistical uncertainties are well estimated, the scaled median absolute deviation (sMAD)¹⁰ of the distribution of values of the parameter α should be 1. The distribution of α as a function of S/N ratio and the value of sMAD at different S/N are shown in red in Fig. A.1.

In order to evaluate how robust our uncertainties are, we test them assuming a correction of the following form: $$\mathrm{\Delta }{v}^{\mathrm{\prime }}=max[(a-b\cdot S/N)\mathrm{\Delta }v,c],$$(A.2)

so that the new uncertainties Δv′ are proportional to the old ones Δv but corrected by a factor that scales with the signal to noise ratio (a − b · S/N). The parameter c is a minimum uncertainty floor for the highest S/Ns. We divide our sample into four different regimes of S/N delimited by the values 0, 19, 38, 57 and 76. Next, we minimize the function, $f={\sum }_1^4(sMAD-1)$, that characterizes how close the sMAD is to 1 in each of the S/N bins, looking for the best a, b and c parameters. For the set of MLM velocities, we find the correction: $$(\mathrm{\Delta }v{)}_{mlm}^{\mathrm{\prime }}=max[(1.448-0.021\cdot S/N)(\mathrm{\Delta }v{)}_{mlm},0.65\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}],$$(A.3)

and for the set of CC velocities: $$(\mathrm{\Delta }v{)}_{cc}^{\mathrm{\prime }}=max[(1.253-0.011\cdot S/N)(\mathrm{\Delta }v{)}_{cc},0.53\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}].$$(A.4)

Both correction methods yield remarkably similar results in terms of sMAD and final uncertainties across different S/N regimes (see Fig. A.1). It is important to consider that the uncertainties (Δv)_mlm are considerably lower compared to (Δv)_cc, but the correction is larger for them. Despite this, the corrected MLM uncertainties $(\mathrm{\Delta }v{)}_{mlm}^{\mathrm{\prime }}$ are still lower compared with the corrected CC ones $(\mathrm{\Delta }v{)}_{cc}^{\mathrm{\prime }}$, meaning that the MLM has more precision for v_los determination. Regarding the floor found, we can test it if we look at our spectral resolution. With the LR8 grating, we are resolving points in the spectra within 7 km s⁻¹ (in the wavelength regime of the CaT), and each CaT line is sampled with ~60 points. So this means we are fitting 180 points when we use the 3 lines. Therefore, assuming that most of the information of the v_los comes from the CaT, the minimum theoretical error given by our spectral resolution should be around $7/\sqrt{180}=0.5\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}$; a value consistent with the floor we are finding. The little offset in both methods can be attributed to jittering.

In the upper panels of Fig. A.1 we show the values of the α parameter for all the pairs of observations taken in the same night, for both the CC (left) and the MLM (right) sets of velocities, before and after applying the corrections in Eq. A.4 and A.3. In the lower panels, we show the values of the sMAD of these distributions of α values. In both cases, at S/N ~25 for the CC and ~45 for the MLM methods, the correction factor is very close to 1. However, when we move from that specific S/N the values of sMAD before the correction start to go away from 1, indicating a bad calibration of the uncertainties. For the corrected uncertainties, we obtain a sMAD close to 1 in all the regimes of S/N in both cases, CC and MLM. Therefore, these new uncertainties reproduce the velocity variation that we see from different exposures more robustly for all the S/N regimes. Since the corrections found in here are for the individual exposures, we combined them weighting by the errors to compute the averaged one.

Appendix B MCMC fitting

We use uniform priors for all the parameters, the range for each of them is shown in Tab. B.1. For the multiple component populations, the priors listed are the ones we used for the first run. This serves to find the main characteristics of each population, however the walkers mix between the different populations in the chains, that is why we had to do a second run restricting the range of metallicities to those we obtained to be the likely ones for each of the components after the first run. Also for the third component the prior for the mean metallicity was changed to [-4,0].

We sample the posterior probability distribution with a number of chains equal to 2 times the number of free parameters, each of them evolved 10000 steps. We burn 2000 of those steps, ensuring model convergence. We also applied a thinin of 32 on the order of the autocorrelation length for the chains.

Fig. A.1 Upper panels: distribution of values of α as a function of the minimum S/N between the pair of velocities compared. In red, we show the values before the uncertainty correction and in blue after it, as indicated in the label. Left panel: CC velocities and uncertainties. Right panel: MLM velocities and uncertainties. Lower panels: sMAD of the α parameter distribution at different bins of S/N. In red, we show the values before the uncertainty correction and in blue after it, as indicated in the label. Left panel: CC velocities and uncertainties. Right panel: MLM velocities and uncertainties.

Appendix C Results with corrected uncertainties

In Tabs. C.1, C.2 and C.3 we show the results of the chemodynamical analysis applied to the different datasets: uncorrected CC, corrected CC and corrected MLM after applying the uncertainty corrections derived in Appendix A. In all the cases the correction does not play any role in the results. This is probably because the process of averaging v_los over different exposures for the same star mitigates its effect.

Appendix D Velocity maps

In order to characterize better whether the low velocity ring structure shown in Fig. 6 is a physical property or whether it is an artifact caused by low-number statistics, we plotted in Fig. D.1 three different v_los maps in which the velocity of each point was computed differently as indicated in the caption. In the left panel, we associate each point in the sky with the mean l.o.s. velocity of the stars closer than 12 arcmins; in the middle panel, for each point we compute the mean for the 25 closest neighbors; finally, in the right panel, we group the stars in Voronoi bins. We can see that the structure remains similar in all the cases. However, as explained in Sec. 4, there are not enough stars in that region to study it in detail.

Fig. D.1 L.o.s. velocity maps of Sculptor. They were generated: for each point, computing the mean velocity in a radius of 12 arcmins (left); for each point, computing the mean velocity of the 25 closest neighbors (middle); binning the stars in groups using Vonoroi bins (right).

References

All Tables

All Figures

Fig. 1 Projection on the sky of the spatial distribution of Sculptor’s member stars, color-coded by the value of the l.o.s induced velocity because of the perspective gradient (see the color bar). The purple arrow shows the direction of the systemic proper motion of Sculptor. The black ellipse has the ellipticity and position angle in Table 1 and semi-major axis equal to twice the semi-major projected half-light radius.

In the text

Fig. 2 Selection function of the dataset. Panel a) shows the continuous selection function used in this work and defined in Eq. (2); panel b) shows the ratio between spectroscopically observed and astrometric member stars for a 15×15 gird in the (R, G) plane. Panel c) shows the selection function integrated over all magnitudes (Eq. (3)) as a solid line and the ratio between observed and member stars for 100 different grids of binning in R.

In the text

Fig. 3 Chemo-dynamical properties of Sculptor. (a) / (d) panels: Normalized metallicity and l.o.s. velocity distribution of the spectroscopically observed member stars of Sculptor fitted by the 3-pop model. The dashed-lines indicate the predicted distribution of our models: red for the MR component, blue for the MP, green for the third population, and black for the sum of the three components as indicated in the label. Bands indicate the 1σ confidence intervals. (b) / (c) panels: Projection of the stars in the vlos vs [Fe/H] / R vs vlos planes. In each panel, the stars are color-coded by the population to which they have the highest probability of belonging, red for the MR, blue for the MP and green for the third population as indicated in the label. The systemic velocity of Sculptor is indicated with straight dashed lines.

In the text

	Fig. 4 Projection on the sky of all stars in the LR8 sample (gray). In each panel we assign stars to each population according to the criteria explained in the text (from left to right population 1,2 and 3) and we color-code them by vlos, as indicated in the color bar.
In the text

	Fig. 5 Upper panel: projected surface number density profile for the MP and MR components, obtained directly by the observations and corrected by the selection function, as indicated in the legend. Lower panel: vlos dispersion profile for the MP and MR components, as indicated in the legend.
In the text

Fig. 6 L.o.s. velocity gradient of Sculptor. Upper panel: L.o.s. velocity map of Sculptor centered in its systemic velocity and generated with Voronoi bins; it is color-coded by velocities as indicated in the color bar. The black dots indicate the position of the stars in our sample; the purple arrow indicates the systemic proper motion of the galaxy; finally, the green line is the rotation axis, with the bands indicating the 1σ error range. Lower panel: velocity profile of Sculptor along the direction perpendicular to the rotation axis (maximum velocity gradient). The dots indicate the value of the mean velocity with errors in certain bins at different distances to the rotation axis, in black using all stars and in gray removing those from pop 3; the solid line shows the rotational gradient of 4.0 km s−1 deg−1, with the bands showing the 1σ error range; finally, the dashed line indicates the systemic vlos of Sculptor.

In the text

Fig. 7 Chemical abundances of stars belonging to 3-pop and UFDs. From top to bottom, we present the [Sr/Fe], [Ba/Fe], and [Mg/Fe] abundances as a function of [Fe/H]. The green points represent pop 3 stars with available abundances from the literature, while the gray points indicate abundances from UFDs, as indicated in the label. Downward arrows indicate upper limits. Iron abundances are from the studies in the literature.

In the text

Fig. A.1 Upper panels: distribution of values of α as a function of the minimum S/N between the pair of velocities compared. In red, we show the values before the uncertainty correction and in blue after it, as indicated in the label. Left panel: CC velocities and uncertainties. Right panel: MLM velocities and uncertainties. Lower panels: sMAD of the α parameter distribution at different bins of S/N. In red, we show the values before the uncertainty correction and in blue after it, as indicated in the label. Left panel: CC velocities and uncertainties. Right panel: MLM velocities and uncertainties.

In the text

	Fig. D.1 L.o.s. velocity maps of Sculptor. They were generated: for each point, computing the mean velocity in a radius of 12 arcmins (left); for each point, computing the mean velocity of the 25 closest neighbors (middle); binning the stars in groups using Vonoroi bins (right).
In the text