Probing the Galactic Halo with RR Lyrae Stars – V. Chemistry, Kinematics, and Dynamically Tagged Groups

We employ a sample of 135,873 RR Lyrae stars (RRLs) with precise photometric-metallicity and distance estimates from the newly calibrated P – φ 31 – R 21 –[Fe/H] and Gaia G -band P – R 21 –[Fe/H] absolute magnitude-metallicity relations of Li et al., combined with available proper motions from Gaia EDR3, and 6955 systemic radial velocities from Gaia DR3 and other sources, in order to explore the chemistry and kinematics of the halo of the Milky Way (MW). This sample is ideally suited for characterization of the inner-and outer-halo populations of the stellar halo, free from the bias associated with spectroscopically selected probes, and for estimation of their relative contributions as a function of Galactocentric distance. The results of a Gaussian Mixture-Model analysis of these contributions are broadly consistent with other observational studies of the halo, and with expectations from recent MW simulation studies. We apply the HDBSCAN clustering method to the speciﬁc energies and cylindrical actions ( E , J r , J φ , J z ), identifying 97 Dynamically Tagged Groups (DTGs) of RRLs, and explore their associations with recognized substructures of the MW. The precise photometric-distance determinations ( δ d/d < 5 %), and the resulting high-quality determination of dynamical parameters, yield highly statistically signiﬁcant (low) dispersions of [Fe/H] for the stellar members of the DTGs compared to random draws from the full sample, indicating that they share common star-formation and chemical histories, inﬂuenced by their birth environments.


INTRODUCTION
RR Lyrae stars (RRLs), which are old (>10 Gyr), low-mass (< 1 M ), generally metal-poor periodic pulsating variables stars located on the instability strip of the horizontal branch, have been used for decades to explore the nature of the Galactic halo of the Milky Way (MW).Their utility as halo tracers is primarily driven by their relatively high luminosities (MV ∼ 0.65; e.g., Catelan & Cortés 2008;Muraveva et al. 2018), and the fact that they are "self selecting"; identified by their variability alone, and not based on spectroscopic follow-up of alternative probes, such as main-sequence turn-off and giant-branch stars, which often require complex targetselection criteria.
The availability of well-sampled light curves from multiple pho-JSPS Research Fellow tometric observations by the Gaia project (Gaia Collaboration et al. 2022b) has enabled large samples of RRLs to be assembled, from which accurate determinations of metallicity have been obtained (Li et al. 2023), based on refined calibrations of the relationships between [Fe/H] and parameters derived from light curves (e.g., the period and Fourier-decomposition parameters φ31 and R21) for type RRab (fundamental-mode pulsating stars) and type RRc (first-overtone pulsating stars), with a scatter of about σ[Fe/H] = 0.2 dex.When combined with the newly calibrated G-band absolute magnitude-metallicity relations of Li et al. (2023), such a sample provides halo tracers from the local neighbourhood out to Galactocentric distances on the order of 100 kpc.Once large-scale substructures (such as the Sagittarius Core and Stream), dwarf galaxy satellites of the MW, and globular clusters (GCs) are identified and removed, the remaining RRLs can be used to infer the relative con-tributions of the inner-and outer-halo populations as a function of distance.
With the addition of proper motions and systemic radial-velocity measurements from Gaia DR3 (Gaia Collaboration et al. 2022a), RRLs also provide a powerful tool for examination of the nature of the halo of the MW based on their derived space motions.An analysis of the energy-action space, derived from positions and velocities, can reveal stars with similar dynamical properties, known as Dynamically Tagged Groups (DTGs), which are stars with common birthplaces that have been stripped from their parent structures (e.g., dwarf satellites and GCs).Previous analyses of DTGs (e.g., Limberg et al. 2021a;Shank et al. 2022a,b) and Chemo-Dynamically Tagged Groups (CDTGs; e.g., Gudin et al. 2021;Shank et al. 2023;Zepeda et al. 2023) have provided insight into the star-formation and chemical-evolution histories of their natal environments.For example, the inner-halo region of the Galaxy is thought to comprise a large amount of stellar debris from at least one ancient merger, the substructure Gaia-Sausage-Enceladus (GSE) (Belokurov et al. 2018;Helmi et al. 2018), and possibly others.Clustering RRLs based on their energies and actions provides the opportunity to associate them with mergers and accreted components that occurred over the full assembly history of the MW (Helmi et al. 1999), making the importance of identifying DTGs crucial for this study.
In Paper I (Liu et al. 2020) of this series, we constructed a catalogue of over 6000 RRLs with metallicity and systemic radial velocity measured from SDSS/SEGUE (York et al. 2000;Yanny et al. 2009;Rockosi et al. 2022) and LAMOST (Cui et al. 2012) spectra.In Paper II (Wang et al. 2022), we performed comprehensive studies of the global structure and substructures of the Galactic stellar halo, based on a sample of about 3000 RRLs.Paper III (Liu et al. 2022) addressed the chemical and kinematical properties of the halo, based on a final sample of about 5000 RRLs.Li et al. (2023) used the spectroscopic catalogue of Paper I to calibrate the photometric metallicities and distances of RRLs from Gaia light curves.This sample, together with RRLs with radial velocity (RV) measurements from Paper I and other sources, was therefore adopted to further halo studies in this series.In Paper IV (Zhang et al., submitted), the origin of the Oosterhoff dichotomy of RRLs is investigated based on samples of stars from large-scale surveys.
We now continue this work based on the much larger initial sample of over 135,000 RRLs from Li et al. (2023), spanning from close to the Galactic Centre to over 100 kpc, derive changes in the metallicity distribution function (MDF) from the inner-halo to the outerhalo region, and analyse DTGs of RRLs in the halo and disk systems of the MW.The large numbers of RRLs in this catalogue with homogeneous metallicity and precise distance estimates can be used to draw a clearer multi-dimensional picture of the Galaxy.We carry out two approaches: (1) analysis of the change in the metallicity distribution function (MDF) of RRLs as functions of radial distance from the Galactic Centre and vertical distance from the Galactic plane, to quantify the relative contributions of the inner-and outer-halo populations of the Galactic halo, and (2) an unsupervised dynamical clustering approach for the subset of RRLs with available full space motions, to identify DTGs that were once members of previous mergers and accreted components.
This paper is organized as follows.In Section 2, we briefly describe the adopted data, the observed properties of the full sample (magnitudes, distances, and metallicities), and compare with the subset of stars with available RVs.We also describe the identification and removal of stars associated with large-scale structures, such as dwarf galaxies (including the Large and Small Magellanic Clouds), the Sagittarius Core (Sgr Core) and Sagittarius Stream (Sgr Stream, including the leading and trailing arms), GCs, and the removal of possible artefacts and potentially problematic stars (e.g., those with possibly compromised astrometry or blended sources, or with very high interstellar reddening), assembling a Cleaned Sample.Section 3 examines the distribution of derived [Fe/H] for stars in the Cleaned Sample, as functions of Galactocentric distance and distance from the Galactic plane.We perform a Gaussian Mixture-Model (GMM) analysis to determine the relative contributions of the inner-and outer-halo populations in each distance bin.Section 4 briefly summarizes our determinations of dynamical parameters for the stars with available RVs. Section 5 describes the dynamical properties for this sub-sample, the identification of DTGs, and their association with known substructures in the halo and disk systems, GCs, and with previously identified dynamical groups.Section 6 describes the chemical structure of the DTGs, and an evaluation of their statistical significance.Section 7 presents a summary of our results, along with future prospects.

DATA
Our analysis is based on the catalogue of refined photometric estimates of metallicity ([Fe/H]) and distances of RRLs recently reported by Li et al. (2023), which also includes proper motions from Gaia EDR3 (Gaia Collaboration et al. 2021), and systemic RVs for 6955 stars, where the vast majority of them are taken from Gaia DR3 (Gaia Collaboration et al. 2022a), and the rest are taken from Liu et al. (2020) and Clementini et al. (2022).The Initial Sample of 135,873 RRLs consists of a total of 115,410 RRab and 20,463 RRc stars.Precise photometric-metallicity and distance estimates were derived by calibrating P − φ31 − R21 − [Fe/H], P − R21 − [Fe/H], and G−band absolute magnitude relations, respectively; for details, see Li et al. (2023).
The global properties of the Initial Sample are presented in Fig. 1, showing the distribution of the G-band magnitude (left panel), heliocentric distance D (middle panel), and photometric metallicity [Fe/H] (right panel) in this sample.Stars lacking available RVs are indicated with black histograms; those with available RVs are shown in blue.In the middle panel, the vertical dashed lines indicate the distances of recognized structures, including the Galactic Bulge, the Sgr Core, the Large Magellanic Cloud (LMC) and Small Magellanic Cloud (SMC), and a number of dwarf galaxies.Note that each panel only includes RRLs in the metallicity range −4.0 <[Fe/H] +0.5.
Fig. 2 presents maps of the Galactic coordinates (l, b) for our RRL stars.The panels show the Initial Sample (top left), the RRLs identified as belonging to recognized structures, including GCs, dwarf galaxies (the LMC, the SMC, and others) and the Sgr Stream and Sgr Core (top right), the RRLs that are removed as part of the trimming procedure described below (bottom left), and the Cleaned Sample of RRLs (bottom right) after these are excised, as described below.
Heliocentric distances and Galactic coordinates are used to derive positions of RRLs in the Cartesian Galactocentric right-hand frame coordinate system, (X,Y ,Z).Fig. 3 shows a map of the Initial Sample of RRLs within 75 kpc from the Galactic Centre.The left panel labels the recognized structures.For the first part of our analysis below, we discard potential members of GCs, the LMC, the SMC, and other dwarf galaxies, the Sgr Stream and Sgr Core, and potential artefacts, as described below, in order to focus on the inner-and outer-halo components of the Galaxy.We retain these structures (but remove the artefacts) for the dynamical analysis.

Removal of globular clusters and dwarf galaxies
Using the catalogue of GCs from Harris (2010) with known equatorial positions (α, δ) and their corresponding half-light radii r h , RRLs are removed that are located within ten times the half-light radius for each GC.As a result, 3085 RRLs are discarded from the Initial Sample.For the dwarf spheroidal satellites, RRLs are removed within an angular distance of 0.5 • from the center of Carina I, Carina II, and Ursa Minor, 0.3 • from Leo I and Leo II, and 1 • from Sculp- tor, Fornax, Draco, and Sextans.The centers of these dwarf galaxies are quoted in Table C.2 of Gaia Collaboration et al. (2018).In the process, a total of 530 RRLs are removed from the Initial Sample.

Removal of the Magellanic Clouds
The LMC and SMC are the largest nearby satellites of the Galaxy, and can be seen clearly in Fig. 3. To remove the likely RRL members of the Magallenic Clouds, we first consider stars that fall within 16 • and 12 • from the LMC and SMC, respectively.Secondly, the subset of those stars with proper motions relative to the LMC or SMC less than ±5 mas yr −1 with respect to the values reported in Gaia Collaboration et al. ( 2018) are kept.Next, we require that RRLs have G > 18.5, based on the expected apparent magnitude of the RRL stars at the distances of the LMC and SMC (Belokurov et al. 2017).Finally, stars are selected within 5 • from the LMC and SMC centers and with distance modulii µLMC > 17.5 or µSMC > 18.0 to treat the dense central region of the Magallenic Clouds.In this process, a total of 22,869 RRLs are removed from the Initial Sample.

Removal of the Sagittarius Stream and Sagittarius Core
Discovered in 1994 (Ibata et al. 1994), the Sgr dwarf galaxy is an ancient relic, crucial for exploring the gravitational potential of the MW and for understanding how its leading and trailing stream arms are influenced by tidal disruption.As seen in Fig. 3

Removal of artefacts
After the removal of potential members of substructures, artefacts are also considered in the trimming process.These artefacts are identified based on the following three criterion: where RUWE is the renormalized unit weight error, BRE is the phot bp rp excess factor, representing the ratio between the combined flux in Gaia BP -and RP -bands (GBP and GRP) and the total flux in the G-band, and E(B − V ) is the reddening.Criterion (i) is applied to remove sources whose astrometric solutions are not wellrepresented by the single-star five-parameter model, as described in Lindegren et al. (2018).Criterion (ii) is applied to identify sources that are considered blends due to their high BRE values (see Evans et al. 2018).Criterion (iii) is applied to remove stars in regions with very high reddening.The number of artefacts in the Initial Sample that meet one or more of the three criteria is 24,619, which are excised from the sample.

The Cleaned Sample
We define the Cleaned Sample of RRLs by application of the trimming procedure described above, resulting in a total of 78,898 stars (see the bottom-right panel of Fig. 2 and the right panel of Fig. 3).We considered whether we could recover some of the trimmed stars by relaxing one or more of the criteria used.For example, criterion (i), which guards against poor astrometry, is not a concern for RRLs with photometric-distance estimates, and for which proper motions do not enter the calculation, and rejection on the basis of criterion (ii), which may be impacted by the nature of the flux distribution within the Gaia photometric bands for strongly variable sources such as RRLs.However, as described below, when we attempt to better isolate the halo populations of the MW by imposing a selection on distance from the Galactic plane of |Z| > 3 kpc, we find that most of the "recovered" stars are eliminated by this additional cut, so this attempt was abandoned.
For our present analysis, we also decided to only include stars with photometric-metallicity estimates in the range −4.0 < [Fe/H] +0.5, on the grounds that the estimated [Fe/H] outside this range may be spurious, and require spectroscopic follow-up to be used with confidence.This cut only removed a relatively small number of stars, resulting in a total of 78,740 RRLs in the Cleaned Sample.

[FE/H] DISTRIBUTION OF RR LYRAES IN THE STELLAR HALO
Fig. 4 shows the MDFs for different radial regions in Galactocentric distance in the stellar halo within the metallicity range −4.0 < [Fe/H] +0.5, after applying the additional |Z| > 3 kpc cut to the Cleaned Sample, which left a total of 42,195 stars.Superposed on each panel are two black dashed vertical lines, which demarcate the metallicity peaks of the inner-halo and outer-halo stellar populations, as derived by Carollo et al. (2007Carollo et al. ( , 2010) ) and Beers et al. (2012) on the basis of kinematic analyses of metal-poor stars from SDSS/SEGUE in a limited local volume.The red dashed line in each panel indicates the position of the median metallicity in each radial bin for comparison.From inspection of this figure, it is clear that there is a general trend for decreasing [Fe/H] with increasing Galactocentric distance, as has been suggested by numerous previous analyses (e.g., Carollo et al. 2007Carollo et al. , 2010;;dee Jong et al. 2010;Beers et al. 2012;An & Beers 2020;Dietz et al. 2020;An & Beers 2021).It is also apparent that the fractions of stars near the black dashed line marking the peak metallicity of the inner-halo population ([Fe/H] = −1.6),while initially roughly constant, begins to drop beyond 20-25 kpc, while the fractions of stars near the peak metallicity of the outer-halo population ([Fe/H] = −2.2) increases beyond this radius.
Fig. 5 shows the MDFs for different vertical regions from the Galactic plane in the stellar halo for this same sample.The black dashed lines and red dashed line are as in Fig. 4. From inspection of this figure, a general trend for decreasing [Fe/H] with increasing distance from the Galactic plane is clearly observed, as inferred by Carollo et al. (2010) based on the kinematics of local halo stars.The fractions of stars near the black dashed line marking the peak metallicity of the inner-halo population, while initially roughly constant, begins to drop beyond 10-15 kpc, while the fractions of stars near the peak metallicity of the outer-halo population increases beyond this distance.
In order to gain more insight about the nature of the stellar populations in the Galactic halo, and their relative contributions as a function of Galactocentric distance and distance from the Galactic plane, we have carried out a GMM analysis, as described below.

Gaussian Mixture-Model Analysis
Similar to Liu et al. (2022), we employ two Gaussian components to model the contribution from the inner-halo (IH) and outer-halo (OH) populations, respectively.The metallicities are assumed to follow a Gaussian distribution with means [Fe/H]IH / [Fe/H]OH and dispersions σIH /σOH.The fractions of IH and OH stars are represented by fIH and fOH, respectively.In this way, for the model parameters Θ = {[Fe/H]IH, σIH, fIH, [Fe/H]OH, σOH}, the likelihood of observing the ith sample star with metallicity [Fe/H] i is given by: The probability PIH([Fe/H] i ) of a star with measured metallicity [Fe/H]i belonging to the inner-halo population can be calculated from: (2) The probability POH([Fe/H]i) of a star with measured metallicity [Fe/H] i belonging to the outer-halo population can be obtained from: The likelihood of a specific bin j is calculated by multiplying the function of equation 1 for a total of Nj stars located within bin j: Li. (4) The posterior distribution of the model parameters is obtained by:  1 and 2 for the radial and vertical bins, respectively.

The Dual-Halo Interpretation
From inspection of the dual-halo model fits shown in Figs. 4 and  5, this simple picture does a credible job capturing the change in the MDF for the radial and vertical regions considered.The initial parameter guesses (the peaks of the MDFs for these components obtained by Carollo et al. 2010), based on a kinematic analysis of local metal-poor stars, leads to fitted mean metallicities that agree with these to within 0.1 to 0.2 dex across all regions.There remain deviations between the observations and the mixture model in some cases, but none sufficiently large that introduction of additional components appears to be required.It should be kept in mind that our  procedure for identifying and removing clear substructure from the data set described above may well be imperfect, so some deviations between the mixture-model predictions and the observations might be associated with this "residual substructure." There originally existed some controversy about whether or not target-selection biases or other systematic errors (such as in distance or proper-motion estimates) in the Carollo et al. (2007) and Carollo et al. (2010) papers that originally presented kinematic evidence for the existence of a dual halo, which were largely refuted by Beers et al. (2012).In Paper III (Liu et al. 2022), on order of 5000 RRLs Numerical simulations of MW-like galaxies have also provided support for a dual-halo interpretation (e.g., Zolotov et al. 2009;Font et al. 2011Font et al. , 2020;;McCarthy et al. 2012;Tissera et al. 2014;Cooper et al. 2015;Pillepich et al. 2015;Rey & Starkenburg 2022).For example, Zolotov et al. (2009) found that the inner regions (r < 20 kpc) of their simulated stellar haloes are composed of both in-situ and accreted components, while stars in the outer region originated from accreted satellites.More recently, the stellar haloes in the ARTEMIS simulations have shown that the stellar-density profiles can be fitted with broken power laws (Font et al. 2020).Their results suggest that the break radii (typically ≈ 20-40 kpc) correspond to the transition between in-situ formation and accreted origin.With a semi-empirical approach, Rey & Starkenburg (2022) pointed out that a prominent last major merger results in a steeper radial profile for the accreted components.These simulations support that stellar haloes in MW-like galaxies tend to have dual haloes.These dual haloes are composed of inner haloes formed by in-situ and accreted components and outer haloes dominated by accreted components.

DYNAMICAL PARAMETERS DETERMINED WITH AGAMA
From the Initial Sample (135,873 RRLs), there are 6955 RRLs with RVs, and 6932 RRLs that also are within the metallicity range −4.0 < [Fe/H] +0.5.After application of the artefact-removal conditions described in Section 2.1.4above, our Final Sample of RRLs with RVs contains a total of 5355 stars.We employ this subset of stars in order to estimate dynamical parameters as summarized here.
The orbital characteristics of the stars with available RVs are determined using their distances, proper motions, and RVs and corresponding errors as inputs to the Action-based GAlaxy Modelling Architecture1 (AGAMA) package (Vasiliev 2019), adopting the Solar positions and peculiar motions described in Shank et al. (2022b) 2 , and the gravitational potential MW2017 (McMillan 2017).Similar to Shank et al. (2022b), the input quantities and their errors are run 1000 times through the orbital-integration process in AGAMA to calculate the cylindrical velocities (vr, v φ , and vz), cylindrical actions (Jr, J φ , Jz), orbital specific energy (E), and eccentricity (ecc), along with their associated errors.We also determine estimates of the maximum distances from the Galactic plane Zmax, the pericentric distance rperi, and apocentric distance rapo.Refer to Shank et al. (2022b) for definitions of these orbital parameters and for the Monte Carlo error calculation treatment.Note that, according to these cal-  Note-This table is a stub; the full table is available in the electronic edition.culations, all 5355 RRLs in the Final Sample are bound within the Galaxy (E < 0 km 2 s −2 ).
A Toomre diagram of the the Final Sample is shown in Fig. 6.For this sample, 58% of the RRLs are on prograde orbits, while 42% are on retrograde orbits.A small fraction of these stars (7%) fall within a radius of 100 km s −1 from the Local Standard of Rest indicated by the red circle in this figure .Of the 5355 RRLs, there are 4324 RRab stars and 1031 RRc stars.For the subset of RRab stars, 59% follow prograde orbits, while 41% follow retrograde orbits, and 6.4% are within the red circle.For the subset of RRc stars, 56% follow prograde orbits, while 44% follow retrograde orbits, and 4.5% are within the red circle.Thus, the fractions of stars on prograde and retograde orbits for both RRab and RRc stars are roughly the justifying our choice to consider them together in our dynamical analysis.Fig. 7 presents histograms of the derived estimates of Zmax, rperi, and rapo.The majority of the stars have their Zmax and rapo within ∼ 200 kpc.The rperi distances are all within ∼ 50 kpc.There are RRLs with Zmax exceeding 200 kpc and rapo with values beyond 200 kpc, which are likely due to errors in these derived quantities primarily arising from the larger input proper-motion errors.

DYNAMICALLY TAGGED GROUPS OF RR LYRAES IN THE MILKY WAY
The Final Sample of 5355 RRLs with orbital parameters are grouped with the clustering algorithm HDBSCAN, which identifies Dynamically Tagged Groups (DTGs) of stars based on their similar dynamical properties (E, Jr, J φ , Jz).It is important to note that the metallicities are not used in this procedure.This algorithm takes into consideration the phase space of the orbital energies and cylindrical actions and uses the density of the data in this space to identify clusters in a hierarchical fashion.For a more detailed description about the HDBSCAN clustering algorithm, we refer the interested reader to Shank et al. (2022b).
The following inputs are set to initiate HDBSCAN: min cluster size = 5, min samples = 5, cluster selection method = 'leaf', prediction data = 'True', Monte Carlo runs set to 1000, and the minimum confidence level is set to 20% for a cluster to pass the confidence test.The min cluster size sets the minimum number of stars allowed in a cluster, and the min samples is set to the min cluster size, which is used to account for noise levels in the data.Setting cluster selection method = 'leaf' permits tighter clustering, and prediction data = 'True' allows HDBSCAN to store memory of nominal clusters for the subsequent Monte Carlo run.We found that setting the min cluster size input to 5, as opposed to a higher number such as 10, has little impact on the mean confidence levels associated with the identified DTGs, while providing a superior resolution of the fine-scale dynamical substructure.
The HDBSCAN clustering routine found 97 DTGs.Table 3 lists the identified DTGs and their corresponding number of stellar members, the confidence levels, the means and dispersions of their metallicities, and associations made with MW substructures, stellar associations, previous group associations, GC associations, and dwarf galaxy associations.Note that the nomenclature for previous group associations for DTGs is adopted from Yuan et al. (2020).If a DTG is not associated with any of the types of association, it is labeled with the word 'new'.The average confidence level for the 97 DTGs is 48.6%.The average confidence level for the DTGs with 10 or more members is 50.1%; for the DTGs with fewer than 10 members, the average confidence level is 45.1%.
Table 7 lists the DTGs identified by HDBSCAN with their associations and respective references, along with the names (using Gaia IDs) and [Fe/H] for each member of the particular DTG.
Table 4 lists the mean cluster dynamical parameters for each DTG, their corresponding number of stellar members, along with mean cylindrical velocities, cylindrical actions, energies, eccentricities, and their corresponding dispersions.
We emphasize that the DTGs of RRLs we identify can also be used to target non-variable stars with similar dynamical parameters for which more complete elemental-abundance studies can be carried out in the future.

Milky Way Substructures
The DTGs are associated with known MW substructures.These associations are determined by the orbital parameters and [Fe/H] for the average of the individual RRLs present in each DTG.If the DTG meets the criteria of a particular substructure, then they are associated with that substructure.We generally adopt the criteria presented in Naidu et al. (2020) to assign associations with MW substructures.As a result, we identify DTGs associated with MW substructures such as the Gaia-Sausage-Enceladus (GSE), Metal-Weak Thick Disk (MWTD), Helmi Stream, Sequoia, Sagittarius Stream, and LMS-1 (Wukong).However, since the current final sample of RRLs with orbital dynamical parameters do not have [α/Fe] ratios available, we remove the Naidu criteria associated with [α/Fe] ratios.This influences how we treat the criteria for the MWTD, so we employ the global rotational velocity range quoted in Carollo et al. (2010) as a criterion.
Table 5 lists the MW substructures that are identified as associations, their corresponding number of stellar members, their means and dispersions of metallicity, and the mean cylindrical velocities, cylindrical actions, energies, and eccentricities.
Table 6 lists the DTGs associated with each identified substructure mentioned in Table 5.
The most populated substructure is GSE with 447 members.The second-most populated substructure is the MWTD with 70 members.The third-most populated substructure is the Helmi Stream with 65 members.The fourth-most populated substructure is Sequoia with 56 members.The fifth-most populated substructure is the Sgr Stream with 34 members.LMS-1 (Wukong) is the leastpopulated substructure, with 8 members.
Fig. 8 shows the Lindblad Diagram (top panel) and projectedaction plot (bottom panel) for the MW substructures listed in Table 6.The Lindblad Diagram describes the energy and azimuthal-action positions (E, J φ ) of the respective DTG members.The projectedaction plot captures the contributions of the cylindrical radial, azimuthal, and vertical components (Jr, J φ , Jz) of the actions.A detailed description of the meaning of null and non-zero actions (Ji = 0 , Ji = 0) is provided in Shank et al. (2022b).
The GSE substructure exhibits the widest range in energies and the most extended radial component compared to the other MW substructures.GSE members also show a null azmimuthal component, with the exception of a handful of DTGs, with DTG-40 overlapping in action space with Sequoia.This DTG is associated with GSE instead of the Helmi Stream due to its high average eccentricity ( ecc > 0.7) and average [Fe/H] not meeting one of the criteria for the Helmi Stream ( [Fe/H] < −1.6).Also, a few GSE members are on planar orbits.Helmi Stream members are all prograde, have a stronger vertical component compared to their radial components, and exhibit more circular orbits.The LMS-1 (Wukong) substructure exhibits a prograde behaviour with a significant vertical component.The Sgr Stream has three distinct positions, one coinciding in the action-space plot with Helmi Stream members; the two others are isolated from the rest of the DTGs in the Lindblad Diagram with the highest energies.They also have more vertical contributions in their orbits, with a handful close to circular in the projected-action space.The MWTD exhibits prograde orbits with roughly equal vertical and radial components, with the exception of some DTGs being more planar in their orbits than circular and vice versa.Finally, Sequoia members show retrograde orbits, and overlap with one DTG of GSE (DTG-40); two of its four DTGs are more planar, while the remaining two have more circular orbits.
It is interesting to note that the large number of DTGs ( 43) as-  In the previous work of Shank et al. (2023), the MWTD and Splashed Disk substructures shared CDTGs overlapping in the energy-action space due to the criteria selected quoted in Naidu et al. (2020), which not only utilizes a metallicity range, but also employs a range in [α/Fe] ratios.Due to the absence of alpha-toiron ratios in our RRL sample, we replace the [α/Fe] criterion in Naidu et al. (2020) with a global rotational velocity criterion 100 < V φ < 150 km s −1 for the MWTD quoted in Carollo et al. (2010).For the Splashed Disk, we employed a velocity range 25 < V φ < 175 km s −1 and a metallicity range −1.0 < [Fe/H] < −0.2 , adopted from An & Beers (2021).The limited numbers of of metalrich RRLs in our sample precluded association with the Splashed Disk.

Statistical Framework
Following Gudin et al. (2021), we perform a statistical analysis of our DTGs to determine how probable the observed abundance dispersions of [Fe/H] would be if their member stars were selected at random from the full set of RRL stars in the Final Sample.
We calculate the statistical significance at three thresholds of the CDF values for [Fe/H] (α ∈ {0.50, 0.33, 0.25}), as well the significance across all α values, resulting in two probabilities: (i) The Individual Elemental-Abundance Dispersion (IEAD) probability, which is the individual binomial probability for specific values of α and [Fe/H], and (ii) The Global Element Abundance Dispersion (GEAD) probability, which is the multinomial probability for [Fe/H], grouped over all values of α.This represents the overall statistical significance.
For a more detailed discussion of the above probabilities, and their use, the interested reader is referred to Gudin et al. (2021).

Results
The number of the 97 DTGs satisfying the α levels of 0.50, 0.33, and 0.25 obtained by our calculations are 66, 54, and 47, respectively, which are individually highly statistically significant, as their IEAD probabilities are each p << 0.001.We underscore that almost half of the sample of DTGs lies below the α = 0.25 level.Unsurprisingly, the GEAD probability (which considers all three levels of α) is also highly significant, p << 0.001.We conclude that the stars in the DTGs, which we identify on the basis of their dynamical parameters alone, share common star-formation and chemical histories, influenced by their birth environments.

SUMMARY AND FUTURE PROSPECTS
We have analysed a large sample of 135,873 RR Lyrae stars (RRLs), identified on the basis of variations in light curves from Gaia DR3, and with precise photometric-metallicity and distance estimates from the newly calibrated P -φ31-R21-[Fe/H] and Gaia Gband P -R21-[Fe/H] absolute magnitude-metallicity relations of Li et al. (2023), combined with available proper motions from Gaia EDR3, and 6955 systemic radial velocities from Gaia DR3 and other sources, in order to explore the chemistry and kinematics of the halo of the MW.
After removal of stars that are potential members of globular clusters and dwarf galaxies, the Magellanic Clouds, and the Sagittarius Stream, and excision of possible artefacts, we obtain a Cleaned Sample of 78,898 RRLs, 78,740 of which fall in the metallicity range −4.0 < [Fe/H] +0.5.
We use the subset of 42,195 stars in this sample located at |Z| > 3 kpc to consider the nature of the MDFs for MW halo RRLs as func-tions of Galactocentric distance and distance from the Galactic plane (centered on the Sun).By application of a Gaussian Mixture Model, under the assumption that the halo system comprises two stellar populations (the dual-halo interpretation advocated by Carollo et al. 2007, 2010, and Beers et al. 2012), we demonstrate that this simple model can account for the observed in-situ distribution of [Fe/H] remarkably well.Some small deviations remain between the model prediction and the observations, which may arise from imperfect removal of recognized substructures and artefacts.
The Lindblad Diagram and the projected-action space of the MW substructures reveal that GSE indeed has the most spread in its energy profile and a strong vertical component.The substructures of Helmi Stream, LMS-1 (Wukong), MWTD, and Sgr Stream are all prograde and have more circular orbits.The MWTD is the most prograde and most metal-rich substructure identified.There are no DTGs associated with the Splashed Disk, as expected, due to the dearth of metal-rich RRLs.Sequoia members are all on retrograde orbits, while there is an interesting split between members that have stronger vertical components compared to their horizontal components, and vice versa, a behaviour not seen in the other substructures.
An analysis of the dispersions of [Fe/H] for the RRL DTGs yields highly statistically significant (low) dispersions of [Fe/H] for the stellar members of the DTGs compared to random draws from the full sample, indicating that they share common star-formation and chemical histories, influenced by their birth environments.
To our knowledge, the MW RRL sample we investigate is the largest to date.It is already yielding important constraints on the nature of the halo system, its substructures, and its formation and evolution.In Paper VI of this series, we plan to further explore the nature of the MW halo system, and make use of the sub-sample of RRLs with available RVs in order to construct measures of its global structure, providing strong constraints on the shape of the dark matter halo.
We anticipate that a much larger sample of RRLs will become available with the next data release from Gaia, anticipated within the next two years.This will not only enable improved estimates of photometric metallicities and distances, as well as provide additional systemic radial velocities, but will expand the number and reach of the sample suitable for in-situ studies of the halo system.

Figure 1 .
Figure 1.From left to right, the three panels display the Gaia G magnitude, the heliocentric distance D, and photometric-metallicity distributions of RRLs with and without RVs, respectively.The blue histograms correspond to the 6933 RRLs with available RVs, while the black histograms are for the 128,489 RRLs without RVs.The vertical dashed lines in the middle panel highlight known distances for Galactic Bulge, Sgr Core, the Magallenic Clouds (LMC and SMC), and dwarf galaxies: Carina I and II, Ursa Minor, Draco, Sculptor, Sextans, and Fornax.These distributions are based on the sub-sample of 135,422 RRLs with metallicities in the range −4 < [Fe/H] +0.5.

Figure 2 .
Figure 2. Top Left: Galactic coordinate map of the Initial Sample of 135,873 RRLs.Top Right: The sample of 32,356 RRLs identified as members of GCs (black), dwarf galaxies (blue), Sgr Stream and Sgr Core members (orange), and LMC or SMC members (red).Bottom Left: The sample of 56,975 RRLs that were removed from the Initial Sample if members are identified as MW substructures as indicated in the top-right panel or possible artefacts (gray), RRLs that might have compromised photometry or astrometry, as described in the text.Bottom right: The Cleaned Sample of 78,898 stars, after removing the sample of RRLs shown in the bottom-left panel.Note that the artefacts only correspond to the gray points in the bottom-left panel.

Figure 3 .
Figure 3. Galactocentric Cartesian map of the RRLs in the Z vs. X plane within 75 kpc.The left panel displays the sub-sample of 135,422 RRLs with metallicities in the range −4 [Fe/H] +0.5, along with labeled dwarf galaxies: Draco and Ursa Minor, the LMC, the SMC, and the Sgr Core, the Sgr leading arm (Sgr LA), and the Sgr trailing arm (Sgr TA).Note that there are more stars outside the ranges shown in the figure.The right panel shows the sample after trimming these labeled structures, the Cleaned Sample.The colour bar corresponds to the number of stars in each pixel, where the pixel size is 0.4 kpc by 0.4 kpc.The black star and red star represent the positions of the center of the Galaxy and Sun, respectively.
(left panel), the Sgr Core is a relatively dense region located at around (X = 20, Z = −5) kpc.Using the sample of RRL Sgr Stream and Sgr Core candidates from Ramos et al. (2020), we cross-match with our Initial Sample via Gaia DR3 source IDs.A total of 5872 RRLs are identified as potential Sgr Stream and Sgr Core members, and removed from the Initial Sample.
where O represents the observables, i.e.,[Fe/H]i, and I(Θ) represents the priors of the model parameters.In this study, a Bayesian Markov Chain Monte Carlo (MCMC) technique is used to obtain the posterior probability distributions of the model parameters.The best-fit values of the model parameters are then delivered by the marginalized posterior probability distributions shown in Figs. 4 and 5 and summarized in Tables

Figure 4 .
Figure 4.The Cleaned Sample of 42,195 RRLs with −4.0 < [Fe/H] +0.5 (x-axis) and |Z| > 3 kpc, divided into seven radial regions in Galactocentric distance within the halo.The y-axis displays the fraction of RRLs in a particular [Fe/H] bin with widths of 0.1 dex.The upper-left corner of each panel lists the region in the stellar halo considered, and the upper-right corner shows the number of RRLs in that distance bin.The two black dashed vertical lines are set at [Fe/H] = −2.2 and [Fe/H] = −1.6, as determined for the inner-halo and outer-halo populations by Carollo et al. (2007, 2010) and Beers et al. (2012) on the basis of kinematic analyses of metal-poor stars from SDSS/SEGUE in a limited local volume.The red dashed line indicates the median metallicity of stars in each radial bin.The transition in median metallicity with increasing radial distance is clearly evident, in particular beyond 20 to 25 kpc.The solid red and green lines indicate the best fits to presumed inner-halo and outer-halo components obtained by a Gaussian mixture-model analysis; the solid black line indicates the sum of these components.The lower-right corner shows the fractions of each component.

Figure 5 .
Figure 5. Same as Fig. 4, but divided into six vertical regions of distance from the Galactic plane within the halo.

Figure 6 .
Figure 6.Toomre Diagram of the Final Sample.The axes are v ⊥ = v 2x + v 2 z vs. vy.The red circle represents stars that are within a radius of 100 km s −1 from the Local Standard of Rest (238.5 km s −1 ), while the vertical blue dashed line represents the division between prograde (vy > 0 km s −1 ) and retrograde (vy < 0 km s −1 ) stellar orbits.For the Final Sample of 5355 stars, 58% and 42% of RRLs follow prograde and retrograde orbits respectively.There are 4324 RRab stars and 1031 RRc stars that make up the Final Sample.For the RRab stars, 59% and 41% follow prograde and retrograde orbits, respectively, while the RRc stars follow 56% and 44% follow prograde and retrograde orbits, respectively.The fractions of stars within the disk region (red circle) out of number of RRLs in the Final Sample that are RRab and RRc stars are 6.4%, and 4.5%, respectively.

Figure 7 .
Figure 7. Logarithmic histograms for Zmax (left), r peri (middle), and rapo (right) for the Final Sample of 5355 RRLs.Note that all stars have r peri that are within ∼ 50 kpc.Most stars fall within ∼ 200 kpc for both Zmax and rapo.For both these derived quantities, stars beyond ∼ 200 kpc primarily arise due to large proper-motion errors.

Figure 8 .
Figure 8. Top panel: Lindblad Diagram of the identified MW substructures.The different substructures are associated with the different colours in the legend.Bottom Panel: The projected-action plot of the same substructures.This space is represented by J φ /J Tot on the horizontal axis and (Jz − Jr)/J Tot on the vertical axis with J Tot = Jr + |J φ | + Jz.For more details on the projected-action space diagram, see Figure 3.25 in Binney & Tremaine (2008).

Table 2 .
1. GMM fitting results for radial bins in Galactocentric distance.[Fe/H] IH , f IH , [Fe/H] OH , and f OH indicate the fitted mean values of metallicity and fractions of the inner-halo and outer-halo components, respectively.GMM fitting results for vertical bins in height above the Galactic plane (centered on the Sun).[Fe/H] IH , f IH , [Fe/H] OH , and f OH indicate the fitted mean values of metallicity and fractions of the inner-halo and outer-halo components, respectively.

Table 3 .
(Roederer et al. 2018ers of stars, confidence, mean abundance and dispersion, and associations.We adopt the nomenclature for previously identified DTGs fromYuan et al. (2020).For example, IR18:E is the first and last initial of the first author (IR)(Roederer et al. 2018), the year the paper was published (2018).Following the colon is the name of the group obtained in that paper (E).