A Galactic Eclipse: The Small Magellanic Cloud Is Forming Stars in Two Superimposed Systems

The structure and dynamics of the star-forming disk of the Small Magellanic Cloud (SMC) have long confounded us. The SMC is widely used as a prototype for galactic physics at low metallicity, and yet we fundamentally lack an understanding of the structure of its interstellar medium (ISM). In this work, we present a new model for the SMC by comparing the kinematics of young, massive stars with the structure of the ISM traced by high-resolution observations of neutral atomic hydrogen (H i) from the Galactic Australian Square Kilometre Array Pathfinder survey. Specifically, we identify thousands of young, massive stars with precise radial velocity constraints from the Gaia and APOGEE surveys and match these stars to the ISM structures in which they likely formed. By comparing the average dust extinction toward these stars, we find evidence that the SMC is composed of two structures with distinct stellar and gaseous chemical compositions. We construct a simple model that successfully reproduces the observations and shows that the ISM of the SMC is arranged into two superimposed, star-forming systems with similar gas mass separated by ∼5 kpc along the line of sight.


INTRODUCTION
Corresponding author: Claire Murray cmurray1@stsci.eduAs one of our nearest neighbors, the Small Magellanic Cloud (SMC) is among the most highly scrutinized galaxies in the Universe.Located at a distance of ∼ 62 kpc (see compliation of historical estimates in de Grijs & Bono 2015), the SMC is close enough that individual stars and interstellar medium (ISM) structures can be resolved.
What makes the SMC especially interesting is its markedly different interstellar conditions to the Milky Way.In particular, the SMC has a low metallicity (∼ 20% Solar; Russell & Dopita 1992), making it an excellent laboratory for deciphering the physics of the ISM at higher redshift (e.g., cosmic noon), where similar pristine conditions are expected to be prevalent.
However, the line of sight structure of the SMC is not well constrained, due in large part to the fact that different observational tracers paint disparate pictures for its geometry.For example, the oldest stellar populations in the SMC are distributed spherically within a radius of 7 − 12 kpc (Subramanian & Subramaniam 2012), and do not show significant signs of rotation (Harris & Zaritsky 2006;Gaia Collaboration et al. 2018a;Niederhofer et al. 2018;Zivick et al. 2018;Niederhofer et al. 2021), except for in the very central region (within ±1 kpc from the SMC optical center; Zivick et al. 2021).Stars with measurable distances (red clump stars, Cepheid variables, RR Lyrae) are greatly dispersed along the line of sight (∼ 20 − 30 kpc depth; Mathewson et al. 1988;Nidever et al. 2013;Scowcroft et al. 2016;Jacyszyn-Dobrzeniecka et al. 2016;Ripepi et al. 2017).In contrast, young main sequence and red giant branch (RGB) stars display a radial velocity gradient indicative of rotation (Evans & Howarth 2008;Dobbie et al. 2014;El Youssoufi et al. 2023).In addition, the velocity distribution of stars in the central SMC exhibits signs of tidal disruption at the hands of the Large Magellanic Cloud (LMC) (Niederhofer et al. 2018;Zivick et al. 2019;De Leo et al. 2020;Zivick et al. 2021;Niederhofer et al. 2021).There is mounting evidence for the presence of distinct substructures along the line of sight (e.g., Hatzidimitriou & Hawkins 1989;Martínez-Delgado et al. 2019;Cullinane et al. 2023).In addition, there is a well-known "bimodality" in distance observed in red clump stars on the Eastern side of the system (e.g., Nidever et al. 2013) and recently traced throughout the SMC (Subramanian et al. 2017;Tatton et al. 2021;Omkumar et al. 2021;Almeida et al. 2023).Furthermore, recent chemical analysis of RGB stars in the SMC exhibit complex distributions in both radial velocity and metallicity ([Fe/H]), suggesting that chemical enrichment in the system "has not been uniform", pointing to either strong influence of star formation (SF) bursts (Hasselquist et al. 2021;Massana et al. 2022) or chemically distinct substructures (Mucciarelli et al. 2023).
In addition to its diversely structured stellar populations, the ISM of the SMC is highly disturbed (Stanimirović et al. 1999).Neutral hydrogen (Hi) emission in the SMC is organized into multiple, distinct peaks in radial velocity.This structure was originally interpreted as the presence of two, distinct sub-systems separated in space (Kerr et al. 1954;Johnson 1961;Hindman 1964).Mathewson & Ford (1984) further showed that the velocities of stars, HII regions and planetary nebulae matched that of the Hi, and concluded that the SMC has been "badly torn" by its interactions with the LMC, in agreement with the leading theoretical models at that time (Murai & Fujimoto 1980).The widely held interpretation is that the low velocity gas must be in front (Hindman 1964;Mathewson & Ford 1984), given that optical absorption is typically associated with the lowvelocity Hi peak (e.g., Danforth et al. 2002;Welty et al. 2012).But the interpretation of these results has not always necessitated the presence of two distinct components.For example, a recent study of gaseous filaments in the SMC establishes that the filaments in the low and high velocity components show a preference for being aligned with each other, which suggests that they must somehow be physically related (Ma et al. 2023).Another alternative interpretation is that the system is a series of overlapping supershells (Hindman 1967;Staveley-Smith et al. 1997).
Despite this complexity, the line-of-sight integrated Hi velocity field of the SMC is broadly consistent with a rotating disk (Stanimirović et al. 2004;Di Teodoro et al. 2019).This disk model has underpinned our concept of the SMC system and has long-defined its total dynamical mass -a key ingredient to numerical simulations.However, the velocity gradient observed in Hi emission is perpendicular to that observed in the radial velocities of young stars (Evans & Howarth 2008;Dobbie et al. 2014).In addition, the distance gradient observed in Cepheids from the North East to South West is oriented ∼ 90 degrees away from the minor axis of the gas disk (Scowcroft et al. 2016).Finally, the 3D kinematics of young, massive stars throughout the SMC, which are young enough to trace the kinematics of their birth clouds in the ISM, are inconsistent with the disk rotation model (Murray et al. 2019).
These observational inconsistencies have traditionally been explained away by the fact that the SMC has been severely disrupted by its recent interactions with the nearby LMC.Although precise observations of stellar proper motions in the LMC and SMC (Kallivayalil et al. 2006b,a;Besla et al. 2007;Piatek et al. 2008;Zivick et al. 2018) indicate that the LMC and SMC are likely on their first infall into the Milky Way (MW) halo (although a second infall scenario is also possible; Vasiliev 2024), spatially-resolved star formation histories show that the two galaxies have been interacting with each   The average Hi brightness temperature spectra (T b (vr)) profiles from spatial bins marked at left (black).Each panel includes shaded grey envelopes denoting the 16 th through 84 th percentile of the T b (vr) within each spatial bin.We observe that the radial velocity structure of Hi emission features multiple, distinct velocity peaks (and typically two, dominant components).
As the more massive of the two galaxies by a factor of ∼ 10 ( Peñarrubia et al. 2016;Erkal et al. 2019), the LMC has survived this eventful interaction history as a nearly face-on rotating disk (van der Marel & Kallivayalil 2014;Gaia Collaboration et al. 2021).As a result, theoretical efforts to reconstruct the formation history of the broader Magellanic System (MS) largely focus on the interaction between the LMC and MW (e.g., Besla et al. 2012;Pardy et al. 2018;Lucchini et al. 2021).However, constraining the precise orbital histories (and futures) of the MCs and their satellites relies on the influence of the SMC (Patel et al. 2020).It is thus critical to develop a model that reconciles these diverse observational constraints.
In this work, we present a new model for gas in the SMC, and show that it is composed of not one, but two structures located at distinct distances along the line of sight, distributed across the inner ±4 • of the SMC.In Section 2, we discuss the data products used in our analysis.In Section 3, we present our analysis procedure, in which we determine the order of ISM components along the line of sight.In Section 4, we present a toy model to describe our observational results.In Section 5, we discuss this new model in the context of the literature and the history of the MS.Finally, we conclude in Section 6.

DATA
In this section, we describe the data products used in our analysis.Specifically, we will combine information from the following sources: • High-resolution observations of the SMC's neutral gas emission at 21 cm from the GASKAP-Hi survey (Pingel et al. 2022), which trace the radial velocity of the bulk gas population in the system.
We also use supplementary information from the APOGEE spectral fits, as well as large-area molecular gas observations of the system to assess its chemical state.

Neutral gas radial velocities
The ISM of the SMC is heavily dominated by atomic gas (Leroy et al. 2007;Jameson et al. 2016), and therefore to trace the bulk distribution of gas, we use the recent survey of 21 cm emission in the SMC from the Galactic Australian Square Kilometer Array Pathfinder (GASKAP-Hi) collaboration (Pingel et al. 2022).This dataset was constructed from 20.9 hours of integration time by ASKAP (Hotan et al. 2021) on the SMC (RA = 00h58m43.280s,Dec = −72d31m49.03s)in December 2019.ASKAP is composed of 36 antennas, each equipped with a phased array feed receiver which combine to give the telescope an instantaneous field of view of ∼ 25 square degrees (Hotan et al. 2021).The resulting data products were combined with single-dish observations from the GASS survey (McClure-Griffiths et al. 2009;Kalberla et al. 2010;Kalberla & Haud 2015) to produce an Hi emission data cube with an angular resolution of 30 ′′ (∼ 10 pc at a distance of 62 kpc), sensitivity to Hi brightness temperature of 1.1 K and a per-channel velocity resolution of 0.98 km s −1 -by far the highestresolution, fully-resolved Hi survey of the SMC to date (Pingel et al. 2022).
In Figure 1, we display the intensity-weighted mean radial velocity of the GASKAP-Hi cube (a.k.a., the first moment) for visualization of the velocity structure of gas in the SMC.There is a clear gradient in radial velocity across the main body of the SMC, which seems to resemble the signature of a rotating disk (e.g., Di Teodoro et al. 2019;Gerrard et al. 2023).In the righthand panel of Figure 1, we display example brightness temperature spectra (T B (v)) drawn from a regular grid of positions overlaid on the first moment map.We observe that these example spectra display complex, multipeaked structure.In particular, across the main body of the SMC, where we see strong emission, there are often two distinct spectral components, one at high velocity (∼ 170 km s −1 ) and one at low velocity (∼ 130 km s −1 ).

Line-of-sight extinction and radial velocity of embedded stars
Although 21 cm emission provides a high-resolution probe of the position-position-velocity structure of the gas, it does not directly probe the line of sight distance axis.To sample this third spatial dimension, we use measurements of extinction by dust toward individual stars.As extinction provides a cumulative measurement of the effect of dust along the line of sight, we can compare extinction values towards stars at different radial velocities in order to determine the relative spatial order of these stars along the line of sight (i.e., "front" and "behind").
As our primary goal is to connect the line of sight structure of the stars with the radial velocity structure of the gas, we make the assumption that the gas, the dust, and the young stars coexist spatially and share bulk kinematics on the same spatial scales probed by each tracer.With these assumptions, we can compare total extinction measured to samples of stars with different radial velocities in order to determine the relative order along the line of sight.Although the assumption that dust and gas are well-mixed is widely accepted, the assumption that the stars trace the same 3D motions as the gas relies on targeting only the youngest stellar populations, whose lives are short enough that their motions are coupled to the motions of their parent gas clouds (e.g., Großschedl et al. 2021).The precise definition of "young" to satisfy this requirement is not well-defined, especially in the low-metallicity, dynamically complex environment of the SMC.In general, stellar populations will lose memory of their birthplaces on a dynamical timescale (i.e., a crossing timescale of about 100 Myr for the SMC).As we will describe below, we endeavor to select stars young enough (ages ∼ 10s of Myr) such that our assumption is valid when averaged across large spatial areas, rather than in individual cases.
For this study, we assemble SMC stars observed by the APOGEE survey and the Gaia mission.Both products include precise radial velocity measurements for thousands of stars, and their selection is described below.

Gaia
To construct an initial sample of SMC stars from Gaia DR3 (Gaia Collaboration et al. 2016, 2023;Babusiaux et al. 2023), we follow the strategy described by Gaia Collaboration et al. (2018b).First, to construct a proper motion filter for SMC stars, we select all stars from DR3 within 3 • of (α 2000 , δ 2000 ) = (13.15833• , −72.8002 • ) with low parallax over error (ϖ/σ ϖ < 10) to minimize foreground contamination.Following Equation 2of Gaia Collaboration et al. (2018b), we compute the orthographic projections of the source coordinates, and further select all stars with mean G-band magnitude G < 19 mag.Next, we compute the median proper motions of these stars in x and y (µ x and µ y ), and exclude all stars with values beyond 4× the robust scatter estimate (Gaia Collaboration et al. 2018b) from the median in both axes.Then, we compute the covariance matrix of µ x , µ y (σ) and ultimately define a filter on proper motions defined as µ T σ −1 µ < 9.21.We acknowledge that genuine SMC stars may be excluded by this cut, as stars outside the inner 3 • may have significantly different proper motions.However, we estimate this to be a small effect on our analysis as our focus is on young stars with kinematics consistent with the bulk SMC population.
To select a complete sample of SMC stars, we select all stars within 8 • of the SMC center (although this is larger than the GASKAP-Hi footprint, here we simply follow the previous analysis of Gaia Collaboration et al. (2018b)) with ϖ/σ ϖ < 10 and G < 20, and apply the proper motion filter defined on the central regions of the SMC above to produce a catalog of 2,006,326 SMC stars.As a last step, we restrict this catalog to the footprint of the GASKAP-Hi SMC cube, yielding 1,687,195 stars.
Since our primary interest is to use stars to probe the velocity structure of gas in the system, we further select only stars with radial velocity measurements from Gaia.This removes the majority of stars in our initial SMC/GASKAP-Hi sample, as radial velocities are only measurable at the SMC distance for the brightest stars.In total, we have 3,783 stars with radial velocities.From these, we select only those within the SMC velocity range (80 < v r < 250 km s −1 ; 3,707 stars).Finally, we convert the Gaia radial velocities from their native barycentric reference frame to the Local Standard of Rest frame to be consistent with the Hi.

APOGEE
In addition to the Gaia stars, we select stars from the second iteration of the Apache Point Observatory Galactic Evolution Experiment (APOGEE-2, Majewski et al. 2017), which was a part of the fourth iteration of the Sloan Digital Sky Survey (SDSS-IV, Blanton et al. 2017).APOGEE-2 consisted of two H-band spectrographs (Wilson et al. 2019), one on the SDSS 2.5m telescope at Apache Point Observatory (Gunn et al. 2006) and one on the 2.5m du Pont Telescope at Las Campanas Observatory (Bowen & Vaughan 1973), the latter of which allows for observations of the SMC due to its position in the Southern hemisphere.SMC stars were targeted through a combination of main survey fields and special fields from "Contributed Programs", as described and detailed in Zasowski et al. (2017); Nidever et al. (2020); Santana et al. (2021).
In this work, we make use of APOGEE radial velocities and metallicities from the 17th Data Release (DR17) of SDSS-IV (Abdurro'uf et al. 2022).APOGEE measures radial velocities by cross-correlating the observed spectra with synthetic spectra, using methodology described in detail in Nidever et al. (2015) and updates described in Abdurro'uf et al. (2022).Jönsson et al. (2020) found that the APOGEE DR16 radial velocities agreed with Gaia DR2 to within 1 − 2 km s −1 depending on the apparent magnitude of the star.We find similar agreement.For the 548 stars in common between APOGEE and Gaia, the RVs agree to within 0.03 +2.2  −1.6 km s −1 (mean and ±1σ errors), which is comparable to the average uncertainty in radial velocity in the Gaia sample of 1.5 km s −1 .APOGEE metallicities are measured using the APOGEE Stellar Parameters and Chemical Abundance Pipeline (ASPCAP, García Pérez et al. 2016), which uses the FERRE code (Allende Prieto et al. 2006) to match the normalized, RV-corrected, and sky-subtracted spectra to a library of synthetic stellar spectra, the generation of which is described in Zamora et al. (2015).Updates to this spectral grid are described in Abdurro'uf et al. (2022), which include spectral grids generated using the SYNSPEC code (Hubeny et al. 2021) with non-local thermodynamic equilibrium corrections from Osorio et al. (2020).We select all APOGEE stars within a radius of 9 • from the SMC center (5,938 stars).This search radius is significantly larger than the GASAKP-Hi footprint, and therefore is aribitrarily selected.We then compute the proper motions in x and y in the same way as for the Gaia sample and exclude all stars with values beyond 4× the robust scatter estimate (Gaia Collaboration et al. 2018b) from the median of the Gaia sample in both axes.We then convert the Heliocentric velocities from APOGEE to the LSR frame and select stars with SMC velocities (80 < v r < 250 km s −1 ) residing within the GASKAP-Hi footprint (2,407 stars).

Young star selection
In Figure 2, we display color-magnitude diagrams (CMDs) for the two samples.We use the Gaia observed magnitudes, BP -RP (a blue-red color index) and G.
In the case of the Gaia sample, the stars with measured radial velocities are all bright (G ≲ 16).In both samples, the stars span a wide range of stellar types, including pulsating variables, red supergiants, and asymptotic giant branch stars (Gaia Collaboration et al. 2021).The stars display a radial velocity gradient across the main body of the SMC with similar magnitude as the gas (Figure 1).
The APOGEE stars span a wider range in magnitude, as this survey included both shallow and deep observations.As a result there is a gap between the faintest and brightest RSGs around G ∼ 15.We will discuss the effect of the magnitude distribution of this sample on our results in Section 3.
From these stars, we need to sub-select according to two conditions: (1) stars that are young enough to support our assumption that their motions trace the bulk motions of their parent gas clouds; (2) stars of a spectral type that is well-understood enough for the line of sight extinction towards them to be accurately inferred with available data.To satisfy these criteria, we select all stars within the region highlighted in Figure 2, as this region of the CMD is dominated by RSGs.As a test of our assumption, we generate a set of synthetic stellar isochrone tables for a range of SMC metallicities ([M/H] = −1.0 to −0.65) from PARSEC (Bressan et al. 2012;Tang et al. 2014;Chen et al. 2014Chen et al. , 2015;;Marigo et al. 2017;Pastorelli et al. 2019Pastorelli et al. , 2020)).In Figure 3, we plot histograms of the ages of stars in our RSG selection (blue), as well as the ages of a range of stellar populations identified in the SMC by Gaia Collaboration et al. ( 2021), including red giant branch stars (RGB), asymptotic giant branch stars (AGB), and blue loop stars (BL).The RSG are of intermediate age (∼ 10−100 Myr), but are among the youngest stars available in the two surveys.In addition, RSGs in the LMC have been shown to closely mimic the kinematics of the Hi (Olsen et al. 2015).
Following the RSG selection, we combine the Gaia and APOGEE samples by cross-matching the catalogs and removing overlapping stars from the Gaia catalog, yielding a total of 1,947 stars (782 APOGEE, 1,165 Gaia).In Figure 2c, we plot the spatial distribution of these samples.Because the APOGEE survey of the SMC was limited to individual pointings, these stars are not uniformly spread across the galaxy as the Gaia stars.We will discuss the effect of the difference in spatial sampling in Section 3.

Estimating Extinction
To estimate the extinction towards each source in the sample in the same manner, we use the "Rayleigh Jeans Color Excess" (RJCE) method (Majewski et al. 2011).RJCE was developed to quantify the amount of reddening towards individual stars based on a single observed color, H − 4.5 µm.These wavelengths sample the Rayleigh-Jeans tail of the stellar SED where red clump and red giant stars have the same intrinsic color.As a result, variations in their observed colors can be attributed to the effect of dust extinction (Majewski et al. 2011).
We note that the RJCE method has not been explicitly validated on the stellar types selected here (RSG).To test the validity of the method, we use the same set of synthetic isochrones from Figure 3 and plot histograms of the H − 4.5µm colors for the RSG sample, compared with the RGB sample (where the color-magnitude selection for SMC RGB is defined in Gaia Collaboration et al. 2021) in Figure 4. We observe here that our RSG selected stars span an even narrower H−4.5µm color range than the RGB stars, granting us confidence that these stars are likely good candidates for the RJCE method.
To compute extinction in K s band, A K , we crossmatch our sample with the catalog of stars from the Spitzer SAGE-SMC survey (Gordon et al. 2011), specifically the "more reliable" catalog of IRAC Epoch 0, 1, and 2 catalog.This catalog includes IRAC 4.5µm photometry, as well as H-band photometry from 2MASS (Skrutskie et al. 2006).We find a match within a search radius of 3 ′′ for 1,934 of our stars.We then compute the extinction in the K s band, A K , using RJCE via Majewski et al. ( 2011) Equation 1, reproduced here, (1) We remove stars with NaN-values (117 stars) or extraneously red 3.6−4.5µmcolors (defined as [3.6−4.5]> 0.15; 10 stars).Our final catalog has 1,801 stars.
Figure 6.AK versus velocity for stars within the same regularly-spaced spatial bins as in Figure 1.The stars are colored by the radial velocity according to the same colorbar as in Figure 1.The average Hi T b (vr) profiles are overlaid in black, and include shaded grey envelopes denoting the 16 th through 84 th percentiles of the T b (vr) within each spatial bin.In general, the stellar radial velocities span the same range as the Hi velocities across the SMC.
To further validate our measurements, we compare A K with the results of Zhang et al. (2023), who estimated stellar parameters and extinction towards the 220 million stars with XP spectra from Gaia DR3 using a forward-modeling approach trained on the LAMOST spectral library.We find that the values of A K are positively correlated with the extinction parameter, E (i.e., a scalar proportional to total extinction), from Zhang et al. (2023) with high significance (Pearson R = 0.52) and the null hypothesis that the two samples are uncorrelated can be rejected with p ≪ 0.001. 1  We emphasize that, as will be elaborated on in the ensuing analysis, our key metric for this work is a differential measurement in extinction, and therefore the precise zeropoint calibration to the RJCE method are not needed for our analysis.In addition, the correction for Milky Way foreground extinction is not needed, as we assume that the foreground extinction does not vary 1 We note that we do not use the Zhang et al. (2023) values in our analysis directly given the uncertainty in validating stellar models trained for Milky Way conditions on the SMC.However, as a test we repeated our full analysis procedure using E instead of A K and found statistically indistinguishable results.This lends confidence that we are not significantly biased by our choice of extinction measurement.
significantly between spatial bins considered here (see Section 3).
In Figure 5, we plot the spatially-binned average radial velocities of the stars in our sample.The stars are distributed according to the regions of highest Hi column density, and display an East-West gradient in v r with the same magnitude as seen in T b (v r ) shown in Figure 1.In Figure 6, we plot A K vs. v r for stars in the same example bins shown, and observe that the stars follow the same radial velocity range of the Hi profiles.Of the 1,801 total stars in the sample, 995 (55%) are at low velocity (i.e., v r < M 1 ) and 806 are at high velocity.

ANALYSIS
The GASKAP-Hi data cube and our selected sample of stars provide two high-resolution probes of the 3D structure of the SMC system.In the following analysis, we match the stars with peaks in T b (v r ) based on their v r (Figure 1) and compare the average extinction towards the Hi velocity structures observed across the SMC.This differential measurement -the average difference in extinction towards the "low" and "high" velocity peaks observed in Figures 1 -determines the order of these components along the line of sight.-therefore, ∆AK < 0 and the "low velocity" gas component is in front of the high velocity component in this direction.Bottom: The average extinction towards stars with "low velocity" (vr < M1) is greater than the average extinction towards stars with "high velocity" (vr > M1)-therefore, ∆AK > 0 and the "low velocity" gas component is behind the high velocity component in this direction.

Matching Stars with Gas
The process of matching the stars with Hi peaks requires a series of selection parameters, which we describe below.

Number of H i components
To identify the number of Hi components and their properties, we follow Murray et al. (2019), who used the first numerical derivative of the brightness temperature spectrum in each pixel of the Hi data cube from McClure-Griffiths et al. (2018) to identify peaks in radial velocity.In Figure 2 from Murray et al. (2019), they conclude that most of the SMC main body features at least two components.However, in detail, this process requires choosing two free parameters: the minimum T B threshold for a component to be considered "significant" (T B,min ), and a Gaussian smoothing kernel width (σ G ) for suppressing the effect of spectral noise on the numerical derivative, which is necessary to determine the number of components along the line of sight.We find that varying these selections (1 < T b, min < 10 K, 2 < σ G < 10) does not have a significant effect on our results, as the majority of the stars are located along high-T b , complex LOS.For the ensuing analysis, we select T b, min = 5 K and σ G = 4 channels.

Spatial binning of stars
Given that the stellar sample is unevenly distributed across the SMC, with most stars lying in the central regions and far fewer in the outskirts (Figure 2c), we use Voronoi binning (Cappellari & Copin 2003), an adaptive binning scheme that allows us to ensure consistent source counts across spatial bins.This process also requires several binning parameters.First is the degree of regular-spaced spatial binning before applying the Voronoi binning algorithm (as in practice, using all ∼ 4 million individual Hi spaxels is prohibitively slow), parameterized by the number of pixels per side per bin (N pix,V ).The second free parameter is the signal-tonoise ratio per bin ((SNR) V ).We select N pix,V = 30, 40, 50, 60 and (SNR) V = 3, 4, 5, 6 and repeat our analysis (below) for each unique combination of these parameters.This allows us to account for the effect of binning choice, while ensuring that there are a uniform number of stars per spatial bin.

Identifying structures via differential extinction
For each combination of selection parameters, in each bin we identify stars at "low" velocity (v r < M 1 ) and at "high" velocity (v r > M 1 ) where M 1 is the intensity-  weighted average radial velocity of the gas along the line of sight.If there is at least one source in the "low" and "high" velocity samples, we compute the average extinc-tion towards the "low" and "high" samples (⟨A K ⟩ low , ⟨A K ⟩ high ).Finally, we compute the difference between these values, defined as "low" minus "high": and assign this value of ∆A K to all stars within the bin.
In Figure 7, we illustrate how this process works with simple cartoons for two cases.When ∆A K < 0, there is less extinction towards the "low" velocity component than the "high" velocity component, which implies that the low velocity component is in front of the high velocity component.Conversely, when ∆A K > 0, there is more extinction towards the low velocity component, which implies that it is behind the high velocity component.
In Figure 8, we plot the average of ∆A K over all 16 unique combinations of binning parameters for each star (⟨∆A K ⟩).In each trial, we find that there is a roughly equal number of stars at low and high velocity within each spatial bin, with a total number per bin that ranges between ∼ 6 and ∼ 40 stars.In addition, the value of ⟨∆A K ⟩ for each star is derived from a non-uniform (Voronoi) binning scheme which ensures that a different set of stars is selected in each trial.This choice allows us to minimize the impact of single outliers on the results.
We observe clear trends in ⟨∆A K ⟩ across the face of the SMC.The points in Figure 8 have sizes that are scaled by the signal to noise in the measurement of ⟨∆A K ⟩.Given the fact that the dust content of the SMC is diffuse (outside of the major star-forming regions), the magnitude of the extinction difference is small (≲ 0.1 mag).However, the fact that the sign of ⟨∆A K ⟩ is constant across large areas (100s of pc, much larger than the spatial bin sizes) indicates that the signal detected at each location is significant.

Selecting "front" and "behind" samples
Building from the spatial trends observed in Figure 8, we select all stars in "front" (i.e., negative ⟨∆A K ⟩ at low velocity and positive ⟨∆A K ⟩ at high velocity) and "behind (positive ⟨∆A K ⟩ at low velocity and negative ⟨∆A K ⟩ at high velocity).
In Figure 9, we plot the spatial distributions of the "front" and "behind" samples.We find 1,014 stars in the front sample and 787 in the behind sample.Of these, 35% and 42% are from APOGEE in the front and behind samples respectively, highlighting the fact that both samples have similar representation from the Gaia-only and APOGEE coverage.

Kinematics
To analyze the kinematics of the front and behind Hi samples, we first correct for the systemic motion of the SMC.We estimate the mean proper motion of the both the front and behind samples to be (µ α,sys , µ δ,sys ) = (0.77, −1.25) mas yr −1 , which is in good agreement with previous estimates of the proper motion of the SMC (e.g., Jacyszyn-Dobrzeniecka et al. 2017;Zivick et al. 2018).We estimate the systemic radial velocity to be the median of the GASKAP-Hi first moment map, equal to v r,sys = 151 km s −1 .We subtract v r,sys , µ α,sys , µ δ,sys from v r , µ α and µ δ in the sample.
For convenience, we next convert the observed coordinates (α, δ) coordinates to a Cartesian frame in the plane of the sky (X, Y ) following van de Ven et al. ( 2006) Equation 1, reproduced here, where ∆α = α − α 0 , ∆δ = δ − δ 0 and (α 0 , δ 0 ) is the assumed center of the SMC (here, from Zivick et al. 2021).In this frame, positive "X" pointing West and positive "Y " pointing North.We convert µ α and µ δ to the Cartesian frame (µ W and µ N ) using the orthographic projection from Equation 2 of Gaia Collaboration et al. ( 2018b), with positive µ W pointing West and positive µ N pointing North.Another important observational effect is "perspective rotation", or the apparent kinematic signature imposed by the systemic motion of the objects with large angular extent van de Ven et al. (e.g., 2006).Correcting for this effect requires an assumed distance to the star along the line of sight, which in this case is not known.Therefore, we incorporate this correction in when we model our results in Section 4.
In Figure 10, we plot the three residual (i.e., systemiccorrected) velocities (v r , µ W and µ N ) for the stars in each component.In the bottom row, we compare the average values as a function of East-West coordinate (X), binned so that there is an equal number of stars per bin.
To compare the velocities between the front and behind samples, we conducted an Anderson-Darling test to compare the distributions of v r , µ W and µ N .The null hypothesis that the two samples are drawn from the same parent population can be rejected with a p-value of < 0.001 for v r , and 0.02 for µ W and µ N .This is confirmed visually in Figure 10, it is clear that the slope of the two distributions of v r are significantly different.
In addition, both front and behind display a strong gradient in µ W .We fit a line to this gradient (Figure 10) and find the slope to be 0.185 (0.187) mas/(yr deg), with an intercept of 0.068 (0.075) km s −1 , in the front and behind components respectively.This gradient is consistent with being caused by the tidal influence of the LMC on the SMC, which we will discuss further in Section 4. To further explore the residual proper motions in the system, we subtract the gradient in µ W and look for signatures of rotation in the residuals.To quantify this, we convert µ W,resid and µ N to radial and tangential proper motions (µ r,resid and µ θ,resid ).In Figure 11, we plot µ θ,resid as a function of radius in both components.Coherent rotation should appear as a nonzero µ θ,resid as a function of radius, with positive or negative values depending on the direction of rotation (e.g., Zivick et al. 2021).We bin the stars so that there is an equal number of stars in each of 10 radial bins (∼ 100 stars), and compute the mean and standard error.
In both components, the residual tangential proper motions are consistent with zero in the innermost radii (r ≲ 1.5 • ).At larger radii in the behind component, µ θ,resid switches from positive to negative, whereas in the front component µ θ,resid remains negative.We interpret these results as there being tentative evidence for rotation in the front component, however the uncertainties are large given that this analysis relies on an assumed galaxy center and orientation of rotation.

Where is the molecular gas?
Resolving the molecular gas content of the SMC is important for gauging the amount of fuel available for star-formation in the system.In addition, the ratio of mass in the atomic and molecular gas phases is sensitive to the influence of physical mechanisms in the system (e.g., radiation fields, shocks, turbulence), and traces the chemical state of the ISM.
The best available tracer for the cold, dense molecular gas phase in the SMC, which represents the conditions necessary for star formation, is carbon monoxide (CO) emission at millimeter wavelengths.Existing CO surveys of the SMC have established that there is very little CO relative to the rich atomic gas reservoir (Leroy et al. 2007;Jameson et al. 2016).In addition, CO in the SMC displays a much simpler radial velocity distribution than the Hi-appearing as a single component with a velocity consistent with one of the "low" or "high" velocity components seen in Hi emission but not the other (Rubio et al. 1991;Israel et al. 1993;Rubio et al. 1996;Mizuno et al. 2001;Muller et al. 2010;Tokuda et al. 2021;Saldaño et al. 2022).In addition, high-resolution observations with the ALMA observatory have revealed a wealth of clouds with weak CO emission and extremely small angular size (< 1 ′′ ; e.g., Tokuda et al. 2021).
To explore the molecular gas content of the front and behind components, we compare our results with observations from Mizuno et al. (2001), who mapped the 12 CO(1-0) transition across the SMC with the NAN-TEN telescope.Although the resolution is coarse (2.6 ′ , ∼ 45 pc at the distance of the SMC), it is the survey with the widest and most contiguous spatial coverage.Due to the large angular size of the SMC, most studies are limited to individual star-forming regions.We compute a 0 th moment map along the declination axis of the NANTEN SMC data cube, and compare the structure in right ascension (α 1950 ) vs. radial velocity (LSR; v r ) with the stars in the front and behind components in Figure 12.
Although there is overlap between the radial velocities in the front and behind components across the middle of the SMC, we observe in Figure 12 that the front component matches the RA vs. radial velocity structure of the 12 CO emission better than the behind However, the interpretation of these results becomes complicated if the dependence of CO emission strength on gas-phase metallicity is considered.The ratio between CO and H 2 is known to decrease steeply as a function of metallicity (Schruba et al. 2012;Accurso et al. 2017).So, if the behind component has significantly lower gas-phase metallicity, CO emission will be much more difficult to detect despite the possible presence of molecular gas.In addition, CO structures at larger distances would also be more difficult to detect due to the effects of beam dilution.

Stellar Metallicity
To explore differences in metallicity, we analyze the APOGEE stars in front and behind components.In addition to radial velocities, APOGEE provides stellar parameter measurements, including stellar metallicity, [Fe/H].The APOGEE survey provides accurate and precise metallicity measurements for RGBs across the MW (see e.g., Jönsson et al. 2020), but there have been no assessments of the accuracy of metallicities derived from RSG spectra.While a detailed comparison of APOGEE RSG metallicities with optical studies is beyond the scope of this work, we use isochrones to determine if the relative metallicities of red supergiants are informative.Figure 13 shows a Kiel diagram colored by metallicity for the APOGEE RSGs (a) and the RSG isochrones described and used in Section 2.2.3 (b).APOGEE finds that the more metal-poor stars are generally cooler and that the more metal-rich stars are generally hotter at a fixed luminosity.Overall, the APOGEE stars show the same qualitative trend of metal-rich stars having higher temperatures at fixed log(g), suggesting that the APOGEE metallicities are at least useful in a relative sense.
We note that the subset of RSGs with low log(g) sitting below the distribution in Figure 13a are all bright (G ≲ 13 mag) and "dusty" (A K > 0.2).These stars are present in both the front and behind components, and excluding them does not affect the analysis presented here (i.e., all observed signals persist).
Histograms of the metallicities derived by APOGEE are shown in Figure 14.Although there is significant overlap, there is evidence that stars in the front component have higher metallicity.A Kolmogorov-Smirnoff test between the two samples rejects the null hypothesis that they were drawn from the same parent population with p = 10 −8 .To check whether the observed difference in metallicity is due to a bias in observed magnitudes (e.g., is the behind component dominated by the brightest, intrinsically lowest-metallicity sources?), we compare the CMDs of the two components.We find that both components have stars spanning the full range of observed magnitudes in the RSG selection (Figure 2).The completeness in our sample at these bright magnitudes is expected to be close to 1, and therefore the difference in metallicity is not likely to be a result of observational selection.
In Figure 15 we plot the spatial distribution of RSGs in the front and behind components, colored by metallicity.The majority of the highest metallicity sources are located in the SW Bar, where the majority of the dust and molecular gas content resides (Jameson et al. 2016).Interestingly, we observe that the lowest metallicity sources are in the behind component and are distributed in the Northern Bar.As a possible explanation, there are known outflows of cold Hi from this region (e.g., McClure-Griffiths et al. 2018), and thus stars which form within outflowing gas may be more metal poor than those forming with the Bar, which is a possibility for both the front and behind components.
In the following section, we construct a simple toy model to constrain the distances to the front and behind components based on their observed kinematics.Although the history of the LMC-SMC system is extremely complex to model, even in the latest hydrodynamical N-body simulations (Besla et al. 2012), our goal is to determine whether a simple, physically-meaningful model can explain our results in broad strokes.
We treat the two components (front and behind) as independent structures experiencing tidal effects from the LMC (evident by the gradient in µ W ; Figure 10) as well as rotation (evident from residual µ θ ; Figure 11).We do not consider the bulk motion of the LMC and SMC in the model, we model all stars in the reference frame of the LMC-SMC system, and calculate the tidal force at a fixed distance to the LMC over time.
For each component, we initialize a population of stars with positions distributed normally with σ = 1 kpc in each of the x, y and z dimensions.This value was chosen to roughly model the observed distribution of stars in X and Y , and we assume for simplicity that the systems are symmetric in Z.The stars are given velocities with random Gaussian amplitude and σ = 20 km s −1 in the x, y, and z directions.This value was chosen based on the average standard deviation in the values of µ N for both samples, converted from mas/yr to km/s assuming a distance of 62 kpc.In the plane of the sky, the LMC is located at a distance of d LMC, sky = 20 kpc, and has a mass of M LMC = 1.4 × 10 11 M ⊙ (Erkal et al. 2019).
Next, we construct a simple prescription for the differential tidal force on the stars due to the LMC.Considering the large uncertainties in the SMC geometry, we only consider a simplistic tidal prescription (i.e., omitting the influence of the Milky Way and angular momentum effects).By integrating the tidal acceleration on the stars over a timescale ∆t, and assuming the distance from the center of the system to the the LMC (r LMC ) is much larger than the distance from each star to the center of mass (r), we can compute the velocity in the LMC direction of each star as a result of tides via: In addition to the tidal velocity, we also include a simple prescription for a rotational component (v rot ) as a function of the observed coordinates, following Di Teodoro et al. ( 2019) for asymptotic velocity v f and scale radius X f (we assume X f = 1 kpc).Of course, this prescription assumes a particular orientation of the rotation in each system, as we do not take inclina-tion or higher-order kinematic effects into account.We emphasize that this is out of scope for this work.Given this simple model, we can project the resulting velocities into the three observable components, v r , µ W and µ N based on the values of the three free parameters, d SMC , ∆t and v f and compare them with the results of Figure 10.Because we declined to correct the observations for perspective rotation, we add this effect to the model (i.e., for which we know the distance to each source), using Equation 6 of van de Ven et al. (2006), reproduced below for clarity: where d SMC is the line of sight distance to the center of the system in our observed reference frame, and is computed from r LMC via, To perform the fit, we use Markov Chain Monte Carlo sampling implemented in the emcee Python package (Foreman- Mackey et al. 2013).Following Bayes' rule, the posterior probability that our observations x (where x ≡ v r , µ W ) are consistent with our model, m(θ) (where θ ≡ d SMC , ∆t, v f ) is proportional to the product of the likelihood P(x|θ) and the prior P(θ), P(θ|x) ∝ P(x|θ) P(θ). (8) Here we assume the likelihood is roughly Gaussian and independent for each of the observed quantities so that, where the subscript x runs between the values of x (i.e., v r and µ W ). We assume a flat prior on each of the three free parameters, and allow them to have values in the range of 50 < d SMC < 80 kpc (used to compute r LMC ), ∆t = 0−1000 Myr and −200 < v f < 200 km s −1 .We initialize the model at the same position for both samples (d smc , ∆t, v f ) = (62 kpc, 600 Myr, 30 km s −1 ), and sample the posterior with 1,000 walkers and 10,000 steps per parameter, excluding a burn-in phase of 1,000 steps.
In Figure 16 we display a "corner" plot of the 2D distributions of the posterior, as well as marginalized 1D PDFs.We compute the 50 th percentile of each 1D PDF and compute the uncertainties based on the 16 th and 84 th percentiles.These results are summarized in Table 1.
In Figure 17, we overlay the contours from one realization of the toy model.We observe that there is good agreement between the observations and the model contours in all panels, which is encouraging given the simplicity of the model (i.e., we did not take into account effects of the components on each other, the Milky Way, or any external effects beyond tidal acceleration from a simplistic point-mass LMC).The shape of the µ W (X) distribution (middle column; Figure 17) is constrained by the d SMC and ∆t parameters, and the shape of the v r (X) distributions (left column; Figure 17) is also constrained by v f .We emphasize again that the structures were modeled independently, without prior information about their order along the line of sight (i.e., flat prior on all model parameters).
It is highly encouraging that we recover the fact that the front component is closer along the line of sight (61 +2 −1 kpc versus ∼ 66 +2 −2 kpc).In addition, our results agree with previous constraints for the distance structure to the young, star-forming components of the SMC.
In particular, we agree with the analysis of Yanchulova Merica-Jones et al. (2017), who inferred the properties of dust along the line of sight by modeling the detailed structure of the red giant branch and red clump traced by UV-IR photometry of the Southwest Bar from the HST SMIDGE survey.They concluded that the dust in that region is localized to a thin layer at ∼ 62 kpc within a broad stellar distribution spread ∼ 10 kpc along the line of sight.The dust in this region is spatially coincident with molecular gas traced by CO emission, which we have shown is consistent with being located in the front component.We find consistent results for the distance to the dust in the front component of ∼ 61 kpc.
The other two free parameters in the model, ∆t and v f are also consistent with expectations.The two components have statistically indistinguishable ∆t, indicating that they have undergone tidal influence of the LMC for the same amount of time, which is reasonable.In the front component, v f is larger than in the behind component.This agrees in broad strokes with observations of µ θ , which is more consistently negative in the front component than the behind.We emphasize that v f is extremely uncertain, as it relies on an assumed rotation axis and structure (i.e., higher-order effects such as inclination are not taken into account).Properly modeling rotation in each system requires a much more sophisticated approach, which we consider to be beyond the scope of this work.We emphasize again that this model was designed to constrain d SMC in particular, which influences all kinematic components (rather than just v r in the case of v f ).

Estimating total H i mass
Based on the results of this simple model, we estimate the total Hi masses of the two structures.The first step is to split the GASKAP-Hi cube into the front and behind structures, which requires an estimate for ⟨∆A K ⟩ across the full survey footprint.To estimate this, we compute the binned average of ⟨∆A K ⟩ in square bins 10 pixels on a side, and for pixels without a measure-  ment, we use the average of the nearest 10 pixels with a measurement.
We then model all spectra in the SMC Hi cube with Gaussian functions using GaussPy+ (Lindner et al. 2015;Riener et al. 2019), an extremely efficient, reproducible method for spectral line decomposition.We use the simplest form of GaussPy+, which leverages the Autonomous Gaussian Decomposition algorithm to identify the positions, widths, and amplitudes of spectral lines based on their numerical derivatives, controlled by only three free parameters: smoothing scales for detecting narrow and broad components (α 1 and α 2 ) and the desired signal to noise ratio (S/N).The optimal values for α 1 and α 2 are chosen by running the algorithm on a training set of known decomposition, and in this case we use the helper functions from GaussPy+ to construct a training set selected from the SMC cube itself, and decompose them using an independent iterative method (for more details see Riener et al. 2019).From this process, we find that the optimal decomposition parameters are: α 1 = 1.08 and α 2 = 4.95 and S/N= 3.
Next, for each pixel, if ⟨∆A K ⟩ < 0, we assign the emission for all fitted Gaussian components with mean velocities less than the first moment to the front component and the rest to the behind component (and vice-versa for ⟨∆A K ⟩ > 0).We assume the average distances to the front and behind structures from the simple model, 61 kpc and 66 kpc respectively, and report the Hi mass in Table 1.
Overall, the structures appear to have roughly the same Hi mass.This is not surprising, given the fact that the amplitudes of the Hi features at low and high velocities are similar throughout the system.Our estimates generally agree with previous estimates for the Hi mass of the SMC, which range from ∼ 3 − 6 × 10 8 M ⊙ (Hindman 1967;Bajaja & Loiseau 1982;Stanimirović et al. 1999;Di Teodoro et al. 2019).
In Figure 18 we display Hi column density maps, computed in the optically-thin limit, for the front and behind components.

DISCUSSION
Taken together, our results paint a picture in which the SMC is composed of two, distinct, star-forming systems which overlap along the line of sight.
At face value, this result is consistent with previous conclusions that multi-peaked Hi velocity profiles of the SMC arise due to the presence of separate sub-systems at different line of sight distances (Kerr et al. 1954;Johnson 1961;Mathewson & Ford 1984).The key difference is that we have shown that the two systems cannot be simply classified as having low or high velocity, as in fact the order of the velocity components along the line of sight changes as you move across the face of the SMC.
Another important difference between this work and previous analysis of the line of sight structure of the ISM is that we probe a wide range of SMC environments.For example, optical absorption lines detected towards SMC stars are typically associated with Hi at low velocities (Danforth et al. 2002;Welty et al. 2012), prompting the conclusion that gas at low velocity must be in front along the line of sight across the system.However, these studies are limited to handfuls of detections, biased towards the regions which feature the highest densities and highest optical depths.Our method allows us to probe the low-density diffuse medium of the SMC, and the smoothness of ⟨∆A K ⟩ (Figure 8) provides supporting evidence that we are sensitive to real changes in the column density structure in different regions.
In addition to our agreement with the only constraint for the dust geometry of the SMC in the SW Bar (Yanchulova Merica-Jones et al. 2017), we have general agreement with the picture painted by older stellar populations in the system.Cepheid variables, stars for which precise distances can be estimated, are highly-elongated (up to ∼ 20 kpc) along the line of sight (Scowcroft et al. 2016).Although they do not show the same spatial distribution as inferred for the star-forming ISM in this work (as they are much more elongated and are generally uniformly distributed), there is clear evidence for the presence of distinct sub-structures with similar average positions along the line of sight, i.e., ∼ 60 and ∼ 66 kpc (Ripepi et al. 2017;Jacyszyn-Dobrzeniecka et al. 2017).These substructures are angled across the field of view such that they overlap along the line of sight when viewed in the plane of sight, making it difficult to conduct detailed comparisons with the morphologies of the front and behind components identified here (which also overlap along the line of sight).In addition, the bimodality in distance observed in red clump stars in the East (Nidever et al. 2013;Subramanian et al. 2017;Tatton et al. 2021;Omkumar et al. 2021;Almeida et al. 2023) underscores the potential for multiple components to the SMC system.The separation between the red clump structures is larger than what we infer here (∼ 10 kpc), and we consider it beyond the scope of this work to determine if they are directly related (as this will require more detailed numerical modeling of the system than performed here).However, overall, we consider it to be encouraging that there is a wealth of ancillary evidence for significant substructure along the line of sight to the system.
In addition to the fact that the reservoir of cold, dense star-forming material, also displays evidence for significant substructure along the line of sight, we find that both systems identified here are actively forming stars.Taking a conservative estimate for the lifetimes of these stars to be 100 Myr (most are much younger; Figure 3), and a rough estimate for the relative speed of the stars with respect to the gas of ∼ 10 − 100 km s −1 , the stars will travel ∼ 0.1 − 1 kpc in their lifetimes.This is only ≲ 20% of the inferred separation of the two systems along the line of sight (∼ 5 kpc), which suggests the RSGs likely formed within the system they are located in today.
One possible explanation for our results is that the two systems are remnants of distinct galaxies.The difference in observed stellar metallicities, as well as the distinct difference between the molecular gas features in the two structures, suggests that they may have evolved in different ways.Although the chemical enrichment history appears similar across the SMC (Almeida et al. 2023), there is evidence for "non-uniform" enrichment in the outskirts of the system (Mucciarelli et al. 2023).
Alternatively, the front structure is the main galaxy (the "SMC Remnant"; Mathewson & Ford 1984) and the behind structure is tidal debris resulting from the interaction with the LMC.The complex nature of the debris field surrounding the SMC, including a bifurcated bridge towards the LMC (Muller et al. 2004), and the bifurcated, trailing Stream (Putman et al. 1998;Nidever et al. 2010), support this picture.Furthermore, although in this work we have made the choice to focus on the properties of the two dominant Hi velocity features, it is clear that in some regions of the SMC there are more than two components.Although these additional components tend to have significantly smaller column density than the two main peaks, their presence highlights the possibility that additional subsystems may be identified.
One particular hypothesis of note is that the behind component is a "counter-bridge" of material pulled from the SMC as a result of its interactions with the LMC (Toomre & Toomre 1972;Besla et al. 2012;Diaz & Bekki 2012).According to numerical simulations of the MS, this structure is expected to sit directly behind the main body of the SMC, and if pulled from the outskirts of the system, may explain its lower observed metallicity.The caveat to this interpretation is that the two systems have roughly equal Hi mass, which is not expected in the counter-bridge model.In addition, the spatial distribution of Hi in the system indicates that the "loop" in the North East, posited to be part of the counter bridge (Muller & Bekki 2007) is actually at the near SMC distance.However, the 3D positions and kinematics of clusters in the SMC periphery have consistent distances with the front and behind components found here, as well as being consistent with the numerical predictions of the counter bridge (Dias et al. 2022).
Ultimately, reconciling these results within the broad landscape of observational constraints for the structure of the SMC requires further work.This includes direct distance constraints to the ISM components, as well as numerical studies of the structure of the SMC in the context of the history of the full MS.However, the implications for our understanding of star formation in low-metallicity dwarf galaxies are already significant.The SMC is used as a template for the evolution of the dust-to-gas ratio with metallicity (e.g., see Clark et al. 2023).Accurately associating the diffuse, atomic gas to the metal-rich ISM is crucial for inferring this parameter, which underlies our understanding of the chemical enrichment of galaxies.

SUMMARY
In this study we investigate the 3D structure of the star-forming ISM in the SMC.We compile a sample of young, massive stars from Gaia DR3 and APOGEE whose kinematics are consistent with distinct Hi velocity components observed by the high-resolution GASKAP-Hi survey.By comparing the line of sight extinction to these stars across the SMC, we find that the order of the low and high Hi velocity components changes between SMC environments.We separate the stellar sample into sub-systems, one in front along the line of sight and one behind.We find that these systems have distinct metallicity distributions, and there is evidence that the molecular gas in the SMC resides mostly in the front component.We construct a simple toy model to explain the kinematics of the two components undergoing the tidal influence of the LMC as well as rotation.We find that the two components lie 5 kpc from each other along the line of sight (average distances of 61 and 66 kpc respectively).These results are broadly consistent with previous estimates of the line of sight structure of the SMC from diverse tracers.

Figure 1 .
Figure 1.Top: Intensity-weighted mean velocity map of the SMC from Pingel et al. (2022).The single contour indicates an Hi column density of 15 × 10 20 cm −2 .The center of the SMC used in this work (defined based on the stellar populations observed by Gaia Zivick et al. 2021) is marked with a magenta cross.Spatial bins are overlaid and identified with numbers.Bottom:The average Hi brightness temperature spectra (T b (vr)) profiles from spatial bins marked at left (black).Each panel includes shaded grey envelopes denoting the 16 th through 84 th percentile of the T b (vr) within each spatial bin.We observe that the radial velocity structure of Hi emission features multiple, distinct velocity peaks (and typically two, dominant components).

Figure 2 .
Figure 2. CMDs of the Gaia (a) and APOGEE (b) samples (before RSG selection).Black contours in each panel denote the 1, 2 and 3σ limits of the total Gaia SMC sample.Red supergiants (RSGs) are indicated by the blue highlighted region in each panel.(c) Spatial distribution of selected RSG stars for the Gaia and APOGEE samples.The Hi distribution is highlighted by a single gray contour at 15 × 10 20 cm −2 (see Figure 1).

Figure 3 .
Figure3.Ages of synthetic stellar populations with a range of SMC-like metallicities.We compare synthetic stars passing our RSG selection criteria (Figure2) with synthetic stars in the RGB, AGB and "blue loop" (BL) (the selection of these populations in color-magnitude space is defined in Section 2.3.2 of GaiaCollaboration et al. 2021).We note that the blue loop sample (BL) overlaps with the RSG selection in Figure1.

Figure 4 .
Figure 4. Synthetic IR color histograms for three subsets of the simulated stellar population: all (black), those passing our Gaia color selection for RSGs (blue), and those passing the Gaia color selection for RGBs defined in (Gaia Collaboration et al. 2021) (red).We observe that the simulated RSGs span an even narrower H − 4.5µm distribution than the RGBs.

Figure 5 .
Figure5.Average radial velocity of the selected stellar catalog from Gaia and APOGEE (vr, LSR) (binned spatially according to the pixel sizes shown).The stars display a radial velocity gradient across the main body of the SMC with similar magnitude as the gas (Figure1).

Figure 7 .
Figure7.Cartoon schematic of the differential extinction estimation process.Top: The average extinction towards stars with "low velocity" (vr < M1) is less than the average extinction towards stars with "high velocity" (vr > M1) -therefore, ∆AK < 0 and the "low velocity" gas component is in front of the high velocity component in this direction.Bottom: The average extinction towards stars with "low velocity" (vr < M1) is greater than the average extinction towards stars with "high velocity" (vr > M1)-therefore, ∆AK > 0 and the "low velocity" gas component is behind the high velocity component in this direction.

Figure 8 .
Figure 8.The difference in extinction towards low velocity gas and high velocity gas (⟨∆AK ⟩), averaged across binning parameters at the location of each source.The symbol sizes are scaled by the signal to noise in ⟨∆AK ⟩ (larger symbols = more confidence).The assumed center of the SMC is marked with an "x" (α0, δ0 = 13.15833,−72.80028Zivick et al. 2021) and the direction towards the LMC is indicated with an arrow.

Figure 9 .
Figure9.Spatial distribution of stars in the "front" structure (⟨∆AK ⟩ < 0 at low velocity and ⟨∆AK ⟩ > 0 at high velocity) and the "behind" structure (⟨∆AK ⟩ > 0 at low velocity and ⟨∆AK ⟩ < 0 at high velocity).As in Figure8, the points are scaled by the significance in the ⟨∆AK ⟩ measurement.

Figure 10 .
Figure10.Comparing the 3D kinematics (vr, µW and µN ; left, middle and right, respectively) as a function of Cartesian East-West (X) position for stars in the front and behind structures, after subtracting the systemic motions of the SMC.The bottom row of panels compares the binned average values for each velocity component in 10 bins of equal size for the front and behind structures, highlighting in particular the difference in vr as a function of X.

Figure 11 .
Figure11.The residual tangential proper motion (µ θ,resid ) as a function of radius in the front and behind components, after subtracting the gradient in µw (Figure10).Here we plot the average (line) and standard error (shading) in bins of equal numbers of stars (6 bins per sample).

Figure 12 .
Figure 12.Comparing right ascension (α1950) vs. radial velocity (vr) in bins of equal number of stars (∼ 100 per bin) in the front and behind components (mean and standard error) with observed structure in molecular gas (NANTEN; Mizuno et al. 2001).The colorbar indicates the integrated 12 CO(1-0) flux density.

Figure 13 .
Figure 13.(a): APOGEE log(g) vs. Teff for the RSGs analyzed in this work, colored by metallicity ([Fe/H]).(b): Binned average log(g) vs. Teff for synthetic stellar isochrones with the same RSG selection.We note that the clump of stars at T eff ∼ 4000 K in panel (a) are all bright, dusty RSGs, whose effect we discuss in the text.component.The correspondence between CO velocities and the front component is particularly true on the Eastern side (α 1950 ∼ 14.5 • ) and in the central region (α 1950 ∼ 12.5 • ).However, the interpretation of these results becomes complicated if the dependence of CO emission strength on gas-phase metallicity is considered.The ratio between CO and H 2 is known to decrease steeply as a function of metallicity(Schruba et al. 2012;Accurso et al. 2017).So, if the behind component has significantly lower gas-phase metallicity, CO emission will be much more difficult to detect despite the possible presence of molecular gas.In addition, CO structures at larger distances would also be more difficult to detect due to the effects of beam dilution.

Figure 14 .
Figure 14.Histograms of [Fe/H] for APOGEE stars in the front and behind components.

Figure 15 .
Figure 15.Spatial distribution of stars in the front (left) and behind (right) systems, colored by metallicity, [Fe/H] (for those with estimates from APOGEE).The single contour denotes an Hi column density of 15 × 10 20 cm −2 .

Figure 16 .
Figure16.Results of the MCMC simple model fit procedure to the front (purple) and behind (green) structures.2D histograms display the 1, 2 and 3σ limits of each posterior distribution.Histograms display 1D marginalized posterior PDFs of each parameter.

Figure 17 .
Figure17.The best-fitting toy model results (black contours, denoting 1, 2 and 3σ limits of the model distributions) overlaid on the observed velocity components of the front and behind structures (Figure10).

Figure 18 .
Figure 18.Hi column density maps of the front (left) and behind (right) components.
the anonymous referee, whose comments improved the clarity of the manuscript.C.E.M. acknowledges insightful conversations with Christopher Clark, Karl Gordon, Peter Scicluna, Lara Cullinane and David French during the development of this project.The authors acknowledge Interstellar Institute's program "II6" and the Paris-Saclay University's Institut Pascal for hosting discussions that nourished the development of the ideas behind this work.J.E.M.C acknowledges STFC studentship.E.D.T. was supported by the European Research Council (ERC) under grant

Table 1 .
Simple Model Fit ResultsNote-Parameters derived from the MCMC fit of our simple toy model to the observed velocities of RSGs in the front and behind systems.