key: cord-0320487-n6nq6yiu
authors: Quirk, A. C. N.; Guhathakurta, P.; Gilbert, K.; Chemin, L.; Dalcanton, J.; Williams, B.; Seth, A.; Patel, E.; Fung, J.; Tangirala, P.; Yusufali, I.
title: The Triangulum Extended (TREX) Survey: The Stellar Disk Dynamics of M33 as a Function of Stellar Age
date: 2022-02-09
journal: nan
DOI: 10.3847/1538-3881/ac5324
sha: a2635cce73450ee9fb28f38513e45045e4602166
doc_id: 320487
cord_uid: n6nq6yiu

Triangulum, M33, is a low mass, relatively undisturbed spiral galaxy that offers a new regime in which to test models of dynamical heating. In spite of its proximity, the dynamical heating history of M33 has not yet been well constrained. In this work, we present the TREX Survey, the largest stellar spectroscopic survey across the disk of M33. We present the stellar disk kinematics as a function of age to study the past and ongoing dynamical heating of M33. We measure line of sight velocities for ~4,500 disk stars. Using a subset, we divide the stars into broad age bins using Hubble Space Telescope and Canada-France-Hawaii-Telescope photometric catalogs: massive main sequence stars and helium burning stars (~80 Myr), intermediate mass asymptotic branch stars (~1 Gyr), and low mass red giant branch stars (~4 Gyr). We compare the stellar disk dynamics to that of the gas using existing HI, CO, and Halpha kinematics. We find that the disk of M33 has relatively low velocity dispersion (~16 km/s), and unlike in the Milky Way and Andromeda galaxies, there is no strong trend in velocity dispersion as a function of stellar age. The youngest disk stars are as dynamically hot as the oldest disk stars and are dynamically hotter than predicted by most M33 like low mass simulated analogs in Illustris. The velocity dispersion of the young stars is highly structured, with the large velocity dispersion fairly localized. The cause of this high velocity dispersion is not evident from the observations and simulated analogs presented here.

The Triangulum Galaxy (M33) is one of two dwarf spirals in the Local Group (Gaia Collaboration et al. 2021) . With a stellar mass of roughly 4.8×10 9 M (Cor-arXiv:2202.04659v1 [astro-ph.GA] 9 Feb 2022 belli et al. 2014) , it is ten times less massive than the Milky Way (MW) and Andromeda (M31) and twice as massive as the Large Magellanic Cloud (LMC). It has a high star formation rate (SFR), ∼ 0.7M yr −1 (Blitz & Rosolowsky 2006) , characteristic of later type spirals. Furthermore, it is a relatively isolated galaxy, making its disk more pristine. As high redshift galaxy observations push down to lower stellar masses, M33 becomes an important local comparison point.

While M33 is sufficiently close to be studied in detail, much remains unknown. For example, there is debate as to its interaction history with M31. Some of M33's properties could be explained by an interaction; such as its highly warped H i disk that extends far beyond the stellar disk, out to 22 kpc (Braun & Thilker 2004; Putman et al. 2009; Semczuk et al. 2018) , the apparent disturbances to the outskirts of M33's stellar disk (Mc-Connachie et al. 2009; McConnachie et al. 2010) , and the mutual burst of star formation 2 − 4 Gyr ago for both M31 and M33 (e.g., Putman et al. 2009; McConnachie et al. 2009; Lewis et al. 2015; Ferguson & Mackey 2016; Semczuk et al. 2018; Williams et al. 2017) . However, new studies using the proper motion of both galaxies from Gaia DR2, the Hubble Space Telescope (HST) and the Very Long Baseline Array (VLBA), and cosmological simulations support the notion that M33 is on its first infall (Patel et al. 2017b,a; van der Marel et al. 2018) and has had no significant interaction with M31.

Furthermore, M33 is more complex then expected for a low luminosity galaxy, which are thought to be dominated by a single component, although the highest mass dwarfs can have a thick disk in addition to a stellar halo (Roychowdhury et al. 2013; van der Wel et al. 2014; Patra Nath 2020; Kado-Fong et al. 2020 . However, M33, is proving to be more complex, with a newly confirmed centrally-concentrated halo (Gilbert et al. in press) and bar (Williams et al. 2021; Smercina & et al. in preparation; Lazzarini & et al. in preparation) and the mysterious warps described above (Braun & Thilker 2004; Putman et al. 2009; Semczuk et al. 2018) .

Resolved stellar spectroscopy is a valuable tool to unlock the history of M33, especially when information from spectroscopic observations is examined as a function of stellar age. This is only possible in the Local Group, in which the distances to galaxies allow us to view their entire disk while resolving individual stars. For example, stellar line-of-sight velocity dispersion as a function of stellar age can point to gradual heating and heating via one event, as seen in M31, where the high velocity dispersion increases monotonically with stellar age, suggesting a constant and violent merger history (Dorman et al. 2015) . Comparing stellar kinematics to gas kinematics can also give insight into a galaxy's dynamical heating history, as was demonstrated in M31 (Quirk et al. 2019) , as the asymmetric drift of stellar populations is correlated with their velocity dispersion in the plane of the disk. Furthermore, comparing this observed asymmetric drift to that in simulations suggests that M31 had a relatively recent 4:1 merger (Quirk & Patel 2020 ). M33's distance of 859 kpc (de Grijs & Bono 2014 ) is comparable to that of M31 (McConnachie et al. 2009 ), so we can measure spectroscopy of individually resolved stars in M33. These techniques may also be capable of constraining the merger history of M33.

We expect M33 to be a particularly interesting target for these studies. Unlike M31's obvious remnants from its violent history, M33 is morphologically less disturbed. It has however a much higher SFR (Blitz & Rosolowsky 2006) and lower disk mass surface density (Corbelli et al. 2014) , placing its disk in a very different regime than M31's. M33 is also the prime environment to observe the effects of internal heating (i.e., bursts from star formation, perturbations from giant molecular clouds and the bar, and density waves from spiral arms) because of its low mass and low inclination (Warner et al. 1973 , ∼ 54 • ). Stellar feedback can cause powerful inflows and outflows of gas and bursts of star formation, which can result in drastic radial migration in low mass galaxies and can even lead to the creation of a stellar halo (e.g., Stinson et al. 2009; Maxwell et al. 2012; El-Badry et al. 2016 ). This internal feedback could have drastically changed the stellar kinematics of M33 since its birth. Stellar disks are fragile (Toth & Ostriker 1992) , and while more massive disks are likely able to survive major merger events (e.g., D'Souza & Bell 2018; Hammer et al. 2018 ) and low mass disks are believed to have to survive many minor events (Helmi et al. 2012) , low mass disks, like that of M33, are unlikely to remain intact after major merging events. This notion paired with the fact that M33 is relatively isolated (∼ 230 kpc from M31), and about half of its stellar mass comes from stars that are ∼ 10 Gyr (Williams et al. 2009 ), means that the disk of M33 is fairly pristine and therefore can give us insight into the evolution of isolated high redshift galaxies.

In this paper, we present the TRiangulum EXtended (TREX) Survey of ∼ 7000 targets, making it the largest stellar spectroscopic survey of M33. The survey spans across the entire disk of M33, out to a deprojected radius of ∼ 11 kpc. It is the first dataset that consists of individually resolved stars that extends across the entire inner and into the outer disk. With this dataset, we examine the kinematics of stars in the disk of M33 as a function of stellar age to measure the dynamical heating of the evolving disk. This analysis, which uses a subset of the total sample, is the first study of disk kinematics as a function of stellar age in M33 using only individual stars and is overall the third of its kind (after the MW and M31). The robust dataset presented here has already been used to confirm the existence of a dynamically hot component in M33 (Gilbert et al. in press) .

This paper is organized as follows. In Section 2 we present our new spectroscopic dataset and archival datasets used in this study. Section 3 describes the separation of stars into broad age bins and the removal of possible halo stars, and Section 4 shows the calculation of local velocity dispersion, rotation velocity, and asymmetric drift. Section 5 highlights a comparison of observed kinematics to the kinematics seen in M33-like simulated galaxies, and Section 6 compares them to the kinematics of M31 and the MW Solar Neighborhood. We summarize the main points of this work in Section 7, and in Section A, we show the kinematics of rare spectral types.

Our study made use of large catalogs of stellar photometry and spectroscopy, as well as imaging data of the M33 gas content. Below we describe the photometric, stellar spectroscopic, and gas imaging catalogs.

We started with large libraries of resolved stellar data, including space-based and ground-based photometry. These, in turn, allowed us to obtain our stellar spectroscopic dataset. We selected targets from a wide range of masses and evolutionary stages, including massive main sequence (MS), massive helium burning (HeB), intermediate mass asymptotic giant branch (AGB), and low mass red giant branch (RGB) stars. The use of these different stellar types is described later in Section 3. These different stages were not prioritized equally over the four observing epochs, see description below and Table 1 . The broad evolutionary stage of a star comes from color magnitude diagrams (CMD) of the photometry catalogs described in Section 2.1.1. We describe this process in detail below.

Our strategy for target selection relies on selecting bright isolated stars from resolved stellar photometry. The high precision of this photometry allows us to target stars in the crowded central regions of M33 with confidence that we were observing isolated stars instead of blended light from multiple sources.

Over the four years of spectroscopic observations, our stellar selection varied in response to the available photometry and evolving scientific opportunities, as stated in Table 1 . We relied on a mix of photometry from the Hubble Space Telescope (HST) and the Canada-France-Hawaii Telescope (CFHT). Targets observed in 2016 were selected using archival HST data with broad bands F475W + F814W, or F606W + F814W or where there was a gap in HST coverage, using data from MegaCam on CFHT with i− and g− bands. The HST fields were observed with the Advanced Camera for Surveys (ACS) and the reduction is described in Williams et al. (2009 Williams et al. ( , 2014 . Each of the 2016 masks overlapped with multiple of these ACS fields. The CFHT/MegaCam data was reduced using the MegaPipe reduction pipeline (Gwyn 2008) . The primary targets for these masks were RGB stars, but also included some AGB and red HeB stars and a small number of MS and BHeB stars.

The 2018 and 2019 slitmasks had targets selected from HST photometry from the Panchromatic Hubble Andromeda Treasury: Triangulum Extended Region (Williams et al. 2021, PHATTER) and CFHT/MegaCam (same as described above). The PHATTER survey observed stars in the Andromeda and Triangulum galaxies with six filter photometry using ACS and WFC3: F275W, F336W, F475W, F814W, F110W, and F160W (Dalcanton et al. 2012; Williams et al. 2021) . The photometric catalogs are described in Williams et al. (2014) ; Williams et al. (2021) .

The availability of six filter photometry allows us to more precisely divide stars into broad age bins, as described later in Section 3. With the six filter photometry, we were able to target a range of stellar evolutionary stages for the 2018 masks: MS, HeB, AGB, and RGB stars. To sample a broad range of stellar ages, we preferentially targeted rarer stars, including a large number of bright HeB stars. These stars were identified from PHATTER CMDs. CFHT/MegaCam data was used to fill in the outer parts of the DEIMOS slitmasks that extended beyond the HST PHATTER footprint, into the low density outer disk where HST resolution is less needed. For the 2019 masks, RGB stars were the primary targets. These stars were identified from PHATTER CMDs if in the PHATTER range. If they came from the CFHT/MegaCam MegaPipe reduction, we prioritized stars with g − i > 0.5 and 21 < i < 22.

The 2020 slitmasks were positioned to probe the outer disk, and they are beyond the PHATTER footprint and any other continuous HST coverage. For these outer slitmasks, we used the catalog from the Pan-Andromeda Archaeological Survey (PAndAS) (McConnachie et al. 2018) . PAndAS used CFHT/Megacam to observe >400 square degrees of sky centered on Andromeda and Triangulum with i− and g− bands. The observations and data reduction are described in McConnachie et al. (2018) . Only objects that were flagged by the PAndAS team to most likely be stars were included in our target list. We prioritized placing RGB stars on these masks (g − i > 0.5 and 20.5 < i < 22.5).

To avoid blends, especially in the crowded central regions, we applied an isolation criterion I neighbor < I target − ( d 0.8 ) 2 + 3 to exclude stars with neighbors that are too close and/or too bright and therefore might contaminate the target's light during spectroscopic observation with DEIMOS (Dorman et al. 2012) . We applied this criterion to all of the photometry catalogs used, although it was most critical for the crowded regions targeted using PHATTER photometry. If a target candidate has a single neighbor that fulfills this criterion, it is excluded from the slitmask.

Even with this criterion, it is possible to have multiple objects in a given slit. The majority of these serendipitous observations have good quality spectra and wellmeasured velocities that did not interfere with the main target's spectrum. However, since we do not have easily paired photometry for these objects, we do not include them in this particular study, although they will eventually be incorporated. The total number of targets is 7684, with 2381 from 2016, 2457 from 2018, 906 from 2019, and 1940 from 2020. Adding in the serendipitous targets would increase the sample by ∼ 27%.

The spectroscopic data comes from four epochs of observing. All observations were taken with the DEIMOS Spectrograph (Faber et al. 2003 ) on the Keck II 10 meter telescope. The program uses thirty-six DEIMOS slitmasks and two different grating setups-one to target a wide range of stellar evolutionary phases and one to target older, redder stars. The first ten slitmasks were observed in 2016 using the 600 line mm −1 grating, which has a resolution of R∼ 2000, and a central wavelength and wavelength range of 7200Åand λ ∼ 4600-9800Å, respectively. This setting allows us to target a wide range of spectral types. In 2018, we observed eleven slitmasks across the central disk using the same configuration. In 2019, we obtained four additional slitmasks of spectroscopic data. The first used the 600 line mm −1 grating configuration, and the remaining three were observed using the 1200 line mm −1 grating (R∼ 6000) to account for additional moon light and to target older stars. With this grating, we used a central wavelength of 7800Å and a wavelength range of λ ∼ 6300-9800Å, focusing on the redder part of the spectrum where RGB stars emit significant flux. In the Fall of 2020, we observed the last eleven slitmasks in the outer disk with the 1200 line mm −1 configuration.

The layout of the thirty six slitmasks can be seen in Figure 1 . Because of the high density of stars in the inner regions, we were able to have slitmasks targeting different stars at the same slitmask location and orientation. Some of the targets on different slitmasks were repeated to get higher signal for faint targets and to help calibrate velocity measurement errors. The positions of the 2019 slitmasks are the same as the two most northern 2018 slitmasks. Each slitmask was observed for approximately 2 hours. Table 1 lists the names, positions, orientations, exposure times, numbers of stars observed, years observed, gratings used, photometry sources for target selection, and the main targets for each DEIMOS slitmask.

The DEIMOS spectra were reduced with the spec2d and spec1d programs (Cooper et al. 2012 . This software has been adapted from its original use for the Sloan Digital Sky Survey to be used on DEIMOS spectroscopy. The resulting onedimensional spectra were flat-fielded and sky subtracted and then cross correlated against stellar template spectra to measure the line-of-sight velocity of the target star (Simon & Geha 2007) . The velocity measurements were confirmed visually using the zspec software. At this step, each measurement is given a quality rating, and rare stars and MW foreground stars are identified (more details below). We then shifted the velocities to the heliocentric frame. We account for possible miscentering of the star in the slit width direction, which causes a systematic wavelength shift. To do so, we calculated the wavelength shift of the atmospheric absorption line at 7600Å. We call this shift the A-band correction and applied it to the measured velocity of the star. We found that the A-band correction varies across the mask depending on the slit's position along the length of the mask, possibly because of a slight positional and/or rotational misalignment of the mask in the sky. To account for the spatial variation, we fitted a polynomial to the A-band velocity as a function of mask position for the stars with the best quality spectra. The polynomial was then applied to all stars to calculate the A-band correction based on the stars' positions along the mask. A typical A-band correction is ∼ −1.3 km s −1 and varies by ∼ 7 km s −1 across a mask.

The systematic uncertainties for the old stars observed with the 600 line mm −1 and 1200 line mm −1 gratings were calculated as in Simon & Geha (2007) ; Collins et al. (2011) , giving 5.6 km s −1 for the 600 line mm −1 and 2.2 km s −1 for the 1200 line mm −1 grating. We also estimate random uncertainties, derived from duplicate velocity cross correlation measurements of RGB stars (1.65 km s −1 for the 600 line mm −1 and 1.85 km s −1 for the 1200 line mm −1 grating). The final error is the result of adding the estimated random uncertainties to the systematic uncertainties in quadrature. We do not yet have enough duplicate observations of young and intermediate age stars to calculate an estimate of the systematic uncertainty. Initial analysis of the duplicate young stars suggest that a typical velocity measurement error for these stars is ∼ 12 km s −1 . For now, we take the velocity errors from the spec1d pipeline (Cooper et al. 2012; Newman et al. 2013 ) for the young and intermediate age stars. A typical velocity error measurement is assumed to be ∼ 6 km s −1 , but this is likely an underestimate of the true velocity uncertainty. We will rectify this in the future after obtaining a larger sample of repeat observations.

MW foreground stars are not identified during target selection. Instead, they are removed from our sample if there is Na i absorption present during visual inspection of their spectra, as this indicates the star is likely a dwarf star (Gilbert et al. 2006) . Once these visually classified foreground stars are removed, we compared the line-of-sight velocities of the remaining stars to a Besancon model (Robin et al. 2003 (Robin et al. , 2014 Robin et al. 2017; Amôres et al. 2017) at the location of the TREX Survey with a color and magnitude cut similar to the range of the targets, shown in Figure 2 . The MW foreground stars have radial velocities that peak at −39 km s −1 , with 57% of the distribution at > −50 km s −1 . Only 18 stars, or 0.70% of our final sample (described in Section 3) has line-of-sight velocities > −50 km s −1 , so our study is largely free of MW contamination.

Targets are also removed from our sample if their spectra suggests the target is an extended source (i.e., background galaxy) or if the quality of the spectrum is too poor to return a well measured velocity. We also eliminate stars with |v LOS | > 500kms −1 or that are miscentered in the slit enough to cause a needed correction on the order of 80kms −1 (Sohn et al. 2007; Simon & Geha 2007 ) based on the polynomial fit estimate. Stars with extreme velocities compared to M33's systematic velocity are unlikely to be members of M33 or to have properly measured velocities. With foreground stars, poor quality targets, and duplicate measurements removed, our spectroscopic dataset consists of 4118 stars in M33. In this specific study, we further narrow our sample based on velocity measurement error (Section 4), CMD placement (Section 3), and, for the older stars, the probability of belonging to the stellar disk component (Section 3.1). With the final sample of resolved spectroscopy, we can examine the line-of-sight velocity for individual stars across the disk in M33.

Unlike stars, which retain their kinematic memories of past dynamical heating, dense gas's comparatively rapid cooling dissipates energy, leaving the gas as a low velocity dispersion tracer of the disk's gravitational potential. In this study, we use velocity measurements of H i, CO, and Hα to make comparisons between the dynamics of the gas and stars in the disk of M33.

The H i measurements from the Very Large Array (VLA) are described in Gratier et al. (2010) and have 5 -25 resolution (FWHM), depending on the specific pointing (see their Table 2 with ages derived from spectroscopic analysis instead of CMD divisions; see Appendix A for more details. We list the adopted average age for each bin in the legend. All magnitudes have been extinction corrected.

The CO(2-1) data were observed using the Institute for Radio Astronomy in the Millimeter Range (IRAM) 30 meter antenna by Gratier et al. (2010) ; Druard et al. (2014) . The observations began in 2008 and build off of those from Gardan et al. (2007) . The angular resolution of the data is 12 with a spectral resolution of 2.6 km s −1 (RMS).

The Hα measurements were observed by Kam et al. (2015) at Observatoire du Mont Megantic with a 1.6-m telescope and the Fabry-Perot interferometer in September 2012, producing an angular resolution of ≤ 3 and a typical velocity measurement uncertainty is ∼ 10 km s −1 (FWHM).

Further details on the observations and data reduction are described in the references listed above. All of the gas velocity measurements have been shifted to the heliocentric frame. Each star corresponds to a specific a single pixel of the gas maps, which allows us to make local and direct comparisons of the gas and stellar kinematics. In Section 4.1, we discuss how we locally average the stellar kinematics to better match the resolution of the gas imaging so that one star does not represent the stellar kinematics of an entire pixel in a gas map.

The extent of the H i dataset is vast, so we are able to pair almost all stars with the nearest (in projection) H i velocity measurement. The percentages of stars paired with H i measurements are 91%, 79%, and 71% for the young, intermediate age, and old stars, respectively. (We discuss the division of stars into age bins in Section 3). The CO and Hα velocity measurement maps have a smaller extent so only stars with a projected radius of ≤ 4 kpc are able to be paired to a CO and Hα measure-ment. For Hα, this includes 24%, 12%, and 4% for the young, intermediate age, and old stars. For the CO, it covers 34%, 25%, and 12% for the young, intermediate age, and old stars.

We divide stars loosely into three age groups based on average stellar age at present day. First we use color magnitude diagrams (CMD). Different regions on a CMD are dominated by stars of different masses and ages. For example, the MS-dominated region we target consists almost entirely of massive stars with short stellar lifetimes. Regions dominated by evolving AGB stars are populated by intermediate mass stars with present day ages older than the main sequence, but not as old as the targeted RGB stars, which occupy a region dominated by older low mass stars with present day ages > 2.5 Gyr.

We can use stellar population models to estimate average ages for each CMD stellar region (Williams et al. 2014; Williams et al. 2021 ). In the rest of the paper, we will refer to stars in the AGB region as "intermediate age," and stars in the RGB region as "old." We combine MS stars, blue HeB, and red HeB, and into a single broad age group that we will refer to as "young." Even though the RHeB stars are far to the red, they are put into the young group because we targeted high mass ones with short lifetimes. See Figure 2 for the approximate location of each stellar lifetime division for our sample.

After the CMD division, we re-categorize weak CN and carbon stars based their spectroscopic information, regardless of their CMD location. Both the intermediate age carbon and young weak CN stars are identified using a combination of visual inspection and machine classification of stellar spectra; they are discussed in greater detail in Appendix A. We assign weak CN stars to the young age group and carbon stars to the intermediate age group because the average age of these stars are consistent with the young and intermediate age group, respectively. We have marked these stars in Figure 2 with black outlines to distinguish them from the CMD divisions.

We assign each broad bin an average age using sim- (2009) for specific star formation histories of regions in M33.) Additionally, these age bins have some overlap and contamination due to the approximate nature of CMD divisions. However, the average ages for each bin are distinct enough to broadly study stellar kinematics as a function of age, which is the goal of this work. We compare the dynamics of these three broad groups and look for trends with stellar age. Gilbert et al. (in press) provide evidence for the existence of a dynamically hot component in M33, using a subset of the spectroscopic dataset used in this paper. They do not find evidence for this component in their young stellar sample, made up of weak CN stars, which are best described by a single dynamically cold component. The stars of the hot component make up ∼ 22% of the total old sample of the TREX Survey, which we correct for using the model described in Gilbert et al. (in press) to remove likely halo contaminants from the old disk population. Gilbert et al. (in press ) model the disk and halo assuming a tight kinematic connection to the H i disk. They compare the line-of-sight velocities of stars to the lineof-sight velocity of the H i at the same radius using the titled ring model in Kam et al. (2017) , rather than individual H i measurements. They model the disk and halo as Gaussians in a transformed line-of-sight velocity space defined as the difference between a star's velocity and the calculated line-of-sight velocity of the disk or halo component at the star's disk location, assuming the fractional rotation speed for that component. This allows each component to rotate at a fraction of the speed of the H i disk model. The best fit model from Gilbert et al. (in press) then returns a probability that a given star's kinematics belong to a dynamically cold component. Although the Gilbert et al. (in press) analysis focuses only on the old stars, we see preliminary evidence that the intermediate age population may host a similar component, and thus we apply the same model to remove candidate kinematically hot AGB stars. The model was run separately on the AGB stars utilizing the same model formalism and procedure used by Gilbert et al. (in press) for RGB stars. We use the probability from the model to keep all intermediate age and old stars with velocities that are at least 80% likely to belong to the dynamically cold component, eliminating velocity outliers and producing a more pure disk-like population. We assume all young stars are disk stars. Figure 3 shows a map of the intermediate age and old stars color coded by probability their kinematics are consistent with a cold component.

Removing stars with disk probabilities which are < 80% eliminates ∼ 14% of the initial intermediate age bin and ∼ 23% of the initial old age bin. For the old stars, the percentage of stars eliminated from the disk sample is consistent with the expected fraction of RGB halo stars from Gilbert et al. (in press ). Future work, utilizing an increased AGB sample, will characterize the kinematics of the AGB population as a whole and explore the nature of the AGB stars which have velocities well removed from the M33 disk velocity. In Section 5 and 6, we explore the implications of not removing possible halo stars.

With the quality cuts and the elimination of old halo stars and intermediate age halo star candidates, our study consists of 952 young stars, 521 intermediate age stars, and 1088 old stars for a total of 2561 stars.

We get line-of-sight velocities for the stars in our sample from the stellar spectroscopic observations. In this section, we describe how we use the line-of-sight velocities to calculate local line-of-sight velocity dispersion, construct rotation curves, and calculate asymmetric drift for each stellar age bin.

We calculate the local velocity dispersion by examining the velocities of neighbors around each star. We start with a 50 radius aperture for selecting neighbors and then grow in radius by 5 until there are at least fifteen stars of the same broad age bin within the circle. The sizes of the circles used are illustrated in the last column of Figure 4 . For the young group, the median radius was 85 , for the intermediate age group the median radius was 120 , and for the old age group the median radius was 100 , which is 0.35 kpc, 0.5 kpc, and 0.42 kpc, respectively, at the distance of M33.

If the circle gets to 300 and still does not contain fifteen members, we do not calculate a line-of-sight velocity average or velocity dispersion at that location. The velocity of the skipped star can still contribute to the analysis of its neighbors however. This smoothing technique is similar to that used in Dorman et al. (2015) , Quirk et al. (2019) , and (Quirk & Patel 2020) . After smoothing, we have local velocity dispersion measurements centered on 879 young stars, 462 intermediate age stars, and 1053 old stars.

The resulting velocity maps are shown in Figure 4 . The second column shows the locally averaged line-ofsight velocities for each age bin. The third column shows the local velocity dispersion, the calculation of which we describe below. The fourth column shows the size of the smoothing circle that was used for each center, along with the size of the smallest and largest smoothing circle used for that age bin.

We calculate the weighted mean of the line-of-sight velocities and the inverse variance weighted root mean square dispersion (RMSE) of the members (Equation 1 and 2). In these two equations, the weights are normalized and derived from the velocity measurement errors (σ err,i in the following equations) such that, n i=1 w i = 1 and w i = σ −2 err,i . The velocity measurement error model is discussed in Section 2.1.2.

We use the inverse variance weighted RMSE as the dispersion. We only consider stars of the same age bin in the averaging and dispersion calculation.

Along with local velocity dispersion, we compute the median velocity dispersion for each of the three age bins, which is reported in Table 2 . We can compare these values to the global models fit from Gilbert et al. (in press ), who find a global velocity dispersion of ∼ 16 km s −1 for the young stars and ∼ 21 km s −1 for the old stars, which is similar to the median local velocity dispersions reported here.

We also show details of the velocity dispersion distributions in Figure 5 . The left panel shows the median value, interquartile range, and outliers for the three age bins across the full extent of the TREX Survey. Overall, we find that velocity dispersion does not vary strongly with stellar age, as the median values of each population do not vary significantly. Furthermore, the median values of each age bin are relatively low and are roughly twice the average dispersion of the H i (Chemin et al. 2020, 8 km s −1 ). The low magnitude of velocity dispersion and the lack of an increase with stellar age are not expected when compared to expectations from simulations of slightly more massive disk galaxies (Martig et al. 2014) or observations of star clusters in M33 (Beasley et al. 2015) . However, these studies do not remove dynamically heated populations that are likely to belong to a halo, so it is not an exact comparison. Martig et al. (2014) find that the shape of the agevelocity dispersion relation for young and intermediate age stars is dependent on recent merger or other heating activity, whereas the velocity dispersion for the oldest stars is more dependent on birth kinematics. They also find that uncertainties in ages can obscure possible trends in the age-velocity dispersion relation. Since the age bins in this work are broad, we could be missing a more subtle trend, but this would not explain the low magnitudes.

While the medians of the distributions are similar, the distributions themselves are broad enough that they may be widened by trends with radius. In the right panel of Figure 5 , we show the distributions of velocity dispersion for each age group broken into an inner (r < 5 kpc) and outer (r < 5 kpc) subgroup. For all age bins, the distribution of velocity dispersion shifts lower in the outer region compared to in the inner region, although the median velocity dispersion for the young stars is higher in the outer region. The number of outliers is higher in the inner region for the young and intermediate age stars.

We look more closely at velocity dispersion as a function of radius in Figure 6 , which shows the distributions of velocity dispersion and radius. There is a clear downward trend in the intermediate age and old stars populations as one goes out to greater radii. This suggests the outer disk is dynamically cooler than the inner disk, which is consistent with the findings of Gilbert et al. (in press) and disk galaxies in general (e.g. Bottema 1993) . The outliers are stars with a velocity dispersion that are at least 1.5 σ from the first or third quartile value, which is the distance marked by the whiskers. Velocity dispersion does not vary significantly with stellar age and is on average higher in the inner region. Koch et al. (2018) find that the velocity dispersion of the H i decreases by ∼ 2 km s −1 from the center of M33 to ∼ r = 8 kpc. For the young stars, any trend is less clear because of the extended area of high dispersion, which we examine below, however there is a slight trend of increasing velocity dispersion with radius.

Although there is not much trend between velocity dispersion and age, there is an extended area of extreme velocity dispersion for the youngest group of stars. Figure 4 shows this extended area of high velocity dispersion in the young star population. We examined this anomalous region in detail to test if the substructure is real or potential velocity measurement artifact. First, we looked at young stars with extreme velocities that contribute to the high velocity dispersion in this extended region. Their velocities are well measured and pass our quality cuts. Secondly, this region is also not dominated by the largest size smoothing circle. Because of this, we think the areas of high velocity dispersion in the young stars are real effects and will investigate it further using duplicate velocity measurements as future work. A smaller portion of this same area is also coincident with high velocity dispersion for the old stars that do have well characterized velocity measurement errors.

We also compare the location of the high velocity dispersion region to a UV map of M33 in Figure 7 . The yellow ellipse in the right panel shows the approximate location of the extended area of high velocity dispersion with respect to the UV emission. Because there is ongoing star formation across the disk, there does not appear to be anything notable about the specific high dispersion area and no correlation between stellar velocity dispersion and whether a location is UV dark or UV light. The left panel shows velocity dispersion as a function of position compared to H ii regions from Hodge et al. (1999) . There also does not seem to be anything coincident about the H ii regions around the high stellar velocity dispersion. Koch et al. (2018) examine the atomic ISM across the disk of M33 and find a filament of H i that extends to 8 kpc across the southern disk (see their Figure 18 ). While the filament overlaps the area of high velocity dispersion reported here, the young velocity dispersion is coincident with a void in the filament. It is further unlikely that the filament is related to the high stellar velocity dispersion because the line-of-sight velocities of the young stars would need to be blueshfited by > 30 km s −1 , which they are not. There are some localized bright areas of H i with broad line widths that are close to the area of high velocity dispersion shared by the young and old stars, but it is unclear if they are related to the high stellar dispersion. It is possible that other phenomenon are causing the high velocity dispersion in the young stars and small number of old stars. For example, there could be unmarked massive gas clouds or a significant number of stellar binaries. Additionally this area could contain substructure from a relatively recent minor merger that lies in front of the disk. If M31 and M33 did have an interaction in the past, perhaps the interaction could have increased the dispersion of the H i and young stellar disk.

We use the weighted average line-of-sight velocities to calculate the rotation velocities of the stars, which we compare to the rotation velocities of gas tracers such as H i, CO, and Hα. To calculate the rotation velocity, we convert the line-of-sight velocity to a circular velocity using the tilted ring model described in Kam et al. (2017) . Kam et al. (2017) divide the H i disk of M33 into forty-nine annuli from r = 0 to r = 96 and measure the position angle (PA) and inclination (i) of each annulus, as tabulated in Table 4 of their paper. The rings in Kam et al. (2017) have a width of 0.2 . We interpolate this table to create thinner rings. We then calculate a star's deprojected distance from M33's center using the PA and i of M33. With that distance, we match the star/gas measurement to a tilted ring from the interpolated table and assign it the corresponding ring's PA and i. We recalculate the deprojected distance, and reassign the star/gas measurement to another ring if needed. We repeat this process twice before adopting the final PA and i for the star/gas measurement.

The above equation projects a star onto a circular orbit. We use v sys = −180kms −1 (Kam et al. 2017) . Theta is the azimuthal angle in the plane of M33, and i is the inclination that comes directly from matching a star/gas measurement to a tilted ring. Theta is calculated with θ = β × [α cos(i T R )] −1 where α = η cos(PA TR ) + ξ sin(PA TR ) and β = ξ cos(PA TR ) − η sin(PA TR ). We use the assigned PA of the ring that a star/gas measurement lies in and the deprojected coordinates of the star/gas measurement centered on M33. Since each star is paired with a gas measurement, the star and the gas measurement share the same deprojected geometric parameters but can have different values for rotation velocities, as they have different line-ofsight velocities.

With the above equations, we construct a rotation curve for each of our broad stellar age bins. The three rotation curves are shown in Figure 8 . These rotation curves show the rotation velocity as a function of deprojected radius for the three stellar age bins and for the H i along the same line-of-sight of each star for the full extent of the TREX Survey. We also plot the inner rotation curves compared to the CO and Hα datasets, which do not extend beyond 4 kpc so are not shown in the full rotation curve. These inner rotation curves are plotted in Figure 9 for the inner 4 kpc for the three age bins and H i, CO, and Hα.

In both Figures 8 and 9 it is clear that the rotation curves of both the stars and the gas show significant scatter. Quirk et al. (2019) demonstrate the deprojection factor of the tilted ring model, or the denominator, approaches zero along the minor axis, regardless of inclination because of the cos(θ) factor, which can explain some but not all of the scatter in the rotation curves. Other sources of scatter in the stellar rotation curves could be from poorly measured velocities, particularly in the young stars, or disturbances from M33's bar Lazzarini & et al. in preparation) and spiral arms. Scatter in the gas rotation curves could reflect the impact of star formation or turbulence in the ISM. The amount of scatter in the stellar rotation curves is largest for the youngest group and decreases with stellar age, which suggests M33's high star formation rate could be causing turbulence not only in the gas but also in the birth kinematics of stars born from that gas. This could also for the stars (gas). The deprojection effects around the minor axis cannot explain the full amount of scatter, especially for the gas. The bar is believed to extend to 0.5 kpc (Williams et al. 2021; Smercina & et al. in preparation; Lazzarini & et al. in preparation) .

be causing some of the extreme velocities in the young stars.

To quantify the scatter, we have added median lines to each stellar and gas rotation curve and have made the marker size proportional to the cos(θ) factor. The scatter cannot entirely be attributed to deprojection effects along the minor axis, especially for the gas, which appears to have the most scatter away from the minor axis. We used binning of 0.5 kpc to remove likely outliers while preserving local discrepancies for the median lines. Any velocity difference between the stellar and gas rotation curves is a visual representation of asymmetric drift and will be explored further in the next section.

We use the stellar and gas rotation curves in Figures 8  and 9 to calculate the asymmetric drift (AD or v a ) of the three broad age bins. Often, AD is defined as the difference between the circular velocity derived from the potential of a galaxy and the rotation velocity of the stars (Strömberg 1946) . For the purpose of this study, we define AD to be the difference between the rotation velocity of the gas and that of the stars (Quirk et al. 2019; Quirk & Patel 2020) . This choice allows us to make local and empirical AD measurements without relying on models of the potential of M33.

AD measurements can be used to measure dynamical heating (Quirk et al. 2019) , as gas, which is collisional, can fairly easily dissipate energy and maintain a low energy orbit (Sellwood & Moore 1998) , while stars retain a non-circular orbit if they have been perturbed onto an eccentric orbit (Leaman et al. 2017; Sellwood & Binney 2002) . In M31, observed AD compared to that in sim- ulated IllustrisTNG-100 M31-like analogs has provided evidence that M31 experienced a relatively recent major merger (Quirk & Patel 2020) . In our own galaxy, AD is routinely used to correct the local standard of rest and to predict rotation curves outside the solar neighborhood before Gaia (Golubov 2014; Huang et al. 2016) . For every star and gas measurement pair, we calculate an AD value, v a , from their respective rotation velocities, using v a = v rot, gas − v rot, . We make AD measurements with respect to H i, CO, and Hα. (This equation is in rotation velocity space, whereas the offset used in Gilbert et al. (in press) for halo star identification is in line-of-sight space.) The CO and Hα measurements only exist for the inner ∼ 4 kpc so we calculated AD in the inner regions with respect to all three kinds of gas but only with respect to H i for the full survey extent (or out to ∼ 11 kpc). The distribution of these values are plotted in Figure 10 and 11, and the median values, width of the distributions, and percentage of outliers are listed in Table 3 stars, which means that on average, the Hα is lagging behind the stars and therefore the stars could be more settled on a circular orbit than the Hα. However, both of these AD values are consistent with 0 within 1σ. We expect AD to increase with stellar age (e.g., Sellwood & Binney 2002; Westfall et al. 2007; Gaia Collaboration et al. 2021; Quirk et al. 2019) , so seeing a high magnitude of AD in the young stars suggests there is something perturbing the young stars on short timescales, which resulted in their kinematics diverging from the gas.

Overall, the old stars tend to have the greatest magnitude of AD, but the distributions of AD values are narrower. In the inner 4 kpc, the young and intermediate age stars have the largest widths, that are at times double that of the old stars for each of the three types of gas used in calculating AD, which is the opposite of what is expected for a canonical model of steady dynamical heating of stars. The old stars are thus more systematically offset from the nearby gas in terms of the kinematics, whereas the young and intermediate stars are more randomly offset from the kinematics of the nearby gas. Over the full radial range, even though the young stars still have the widest AD distribution, the three groups of stars have more similar widths. This suggests the stars in our sample are most similar between radii 4 an 11 kpc.

The percentage of outliers (stars with AD values that are 1.5σ or greater away from the 16 th or 84 th percentiles) in each distribution shows no clear pattern with median AD, width of the distribution, or stellar age. The intermediate age stars have the highest percentage of outliers in the distribution of AD values calculated with respect to H i in the inner 4 kpc.

It is important to note that the three gas datasets used come from different telescopes and were reduced using different methods. Subtle differences in the deriva- Table 4 . Stats for the distribution of AD for the three age bins with respect to H i, Hα, and CO for the inner 4 kpc. Median values of AD, width (sigma) of the distribution, and the percentage of outliers in the distribution are shown. The errors on the median value represent the difference between the 16 th and 84 th percentiles divided by √ N, where N is the number of stars. The outliers are stars with a velocity dispersion that are at least 1.5 σ from the first or third quartile value. tion of the gas velocity measurement could be causing variations in the AD median values, not just physical phenomena. However, subtle differences would systematically affect the three age bins equally, which means these potential differences cannot be obscuring a lack of trend between age and asymmetric drift. In summary, M33 does not show the expected increase in AD with stellar age because of the low AD of the intermediate age stars and the high AD value for the young stars. The distributions of AD get narrower with stellar age, showing that the younger stars are more randomly offset from the gas than the older stars.

We analyze M33-like analogs from the IllustrisTNG50-1 cosmological simulation to comment on the uniqueness of the results presented above. The IllustrisTNG Project is a suite of N-body and hydrodynamic simulations from redshift z = 127 to z = 0. It uses the AREPO movingmesh code (Marinacci et al. 2018; Naiman et al. 2018; Springel et al. 2018; Pillepich et al. 2018; Nelson et al. 2018; Springel 2010) . The simulations have a cosmological volume of (51.7 Mpc) 3 and have the following ini- We use data from the IllustrisTNG50-1 simulation (hereafter IllustrisTNG), the smallest but highest resolution simulation in the project that includes both dark matter and baryons. IllustrisTNG follows the evolution of 2160 3 dark matter particles and 2160 3 hydrodynamical cells, resulting in a a baryonic mass resolution of m bary = 8.1 × 10 4 M and a dark matter particle mass resolution of m DM = 4.4 × 10 5 M . Halos and subhalos are identified using SUBFIND (Springel et al. 2001; Dolag et al. 2009 ).

We identify general M33-like analogs as halos that fit the following criteria: the subhalo is the central/primary subhalo at z = 0 in Friend of Friend (FoF) group with a virial mass of M vir = 1 − 2.5 × 10 11 M (Patel et al. 2018 , and references within); the subhalo has a stellar mass of M = 2.8 − 5.6 × 10 9 M (Guo et al. 2010) , and the subhalo's maximum dark matter particle circular velocity is < 70 km s −1 . These cuts produce a sample of 224 analog galaxies. We eliminated eight of these due to lack of baryon particles across the inner 10 kpc, leaving 216 M33-like galaxy analogs.

For each analog, we rotate the particles/gas cells to be face-on so that the "line-of-sight" direction is the z component for all analogs. We exclude particles with |z| > 10 kpc and with a galactic radius of > 20 kpc to target star particles that likely belong to the disk. We then locally average the velocities to mimic the resolution from observations, following Section 3 of Quirk & Patel (2020) . We divide the stellar particles into four groups based on exact stellar age: < 1 Gyr, 1-5 Gyr, 5-10 Gyr, and > 10 Gyr. Unlike the rough age division used for the observational part of the analysis, these age bins use a star particle's exact age and spans the full cosmological time. We use the gas cells to calculate AD for each age bin as well.

For each analog, we use the smoothed kinematics to construct an azimuthally averaged rotation curve using the equations below.

We limit the analysis to the 2D xy plane to mimic the observed line-of-sight velocities. We choose to azimuthally average the rotation curve for each. To azimuthally average the rotation curves, the particles/cells are placed into 0.1 kpc bins based on their distance from the center of the analog. We then take the median of the rotation velocities in each bin for the final rotation curve. This process creates a rotation curve similar to the median lines shown in Figures 8 and 9 .

To calculate a single set of rotation curves and median AD values as a function of stellar age for the entire analog sample, we first calculate both the rotation curve and AD for each individual analog, as described above. Since the rotation curves are radially binned, the cumulative rotation curve for the entire sample of analogs is made by using the median rotation velocity for each stellar particle age group at each radial bin. We calculate AD values for the entire sample using this cumulative rotation curve. The set of cumulative rotation curves is shown in the left panel of Figure 12 . Also shown in this figure are the median rotation curves for the observed M33 young stars, intermediate age stars, and old stars in the right panel.

The set of observed M33 rotation curves is consistent with the rotation curves for star particles with ages 1−5 Gyr from the simulated analogs. However, these rotation curves are not a perfect comparison, as the observational one comes from a tilted ring model, and the simulated ones come from calculating the actual tangential velocity of the stellar particles.

The cumulative AD values for the simulated M33 analogs as a function of stellar particle age are shown in Figure 13 along with the observational median M33 AD measurements. For the simulated analogs, AD increases with stellar age. A gradient in AD with stellar age is also seen in M31-like simulated analogs in the same simulation (Quirk & Patel 2020) . There is no such gradient in the observed AD in M33.

For the young and intermediate age stars, the observed AD is higher than what is seen in simulated M33-like analogs. The AD for the oldest stars is consistent within the observed age errors. While some of the differences in the observed M33 AD and the simulated analogs AD is most likely from the differences in the way AD is calculated for each, it is likely that this would affect all of the age bins equally. Thus, the difference in AD between the observed intermediate and the simulated analogs could be a nonphysical systematic offset. However, the observed young stars have an AD value that is more offset from the AD of simulated analogs than the offset for the other two age groups. Thus even if there is a systematic offset, the young observed stars still break the gradient in age and AD that is seen in the simulated M33-like analogs.

The high AD observed in the young stars in M33 suggests that some phenomenon is dynamically heating the young stars. Of the M33-like analogs, seven have young stars with an AD greater than 10 km s −1 . Six have young star AD values between 12 and 19 km s −1 , and one has an AD value of 33 km s −1 for the young stars. Of these seven, four have AD distributions that are similar for all stars 0-10 Gyr (the age bins that best represent the observed bins). For three analogs, the intermediate age stars (1) (2) (3) (4) (5) have the lowest AD. Upon visual inspection of the star and gas particles used in these analogs, only one analog shows a slight disturbance in the disk, and only one has had any mergers (minor or major) in the past 4 Gyr. Thus, there is no immediate connection between these analogs that could point towards a cause of the high AD in the young stars. Figure  13 shows the AD values for the analog that is closest to the observed AD values. This analog has not had any mergers within the past 4 Gyr and does not show a disturbance or asymmetry in its disk. Closer inspection of the stellar assembly history of this analog is future work. Compared to the whole sample of analogs, the AD from observations is significantly higher for the youngest stars, than that is seen in the simulated analogs.

those for the simulated star particles were not based on a kinematic model. Including all halo stars roughly doubles the AD value for the intermediate age and old bins (see Section 6 and Table 6 ). It is possible that the gap between the AD in the observations and the simulated M33-like analogs is larger for the intermediate age and old stars, but even if the simulated AD halved for the older bins, the largest gap would still be for the young stars.

Resolved stellar spectroscopy of individual stars allows us to study a galaxy in detail. For example, one can examine dynamics as a function of age, spatial position, and metallicity, which puts narrow constraints on origin and evolution. Studies of this kind have already proved valuable for the comparative evolution of the LG's two massive spirals, M31 and the MW. In spite of having similar masses (M vir ≈ 10 12 M ), the disk kinematics of the MW and M31 differ significantly. The MW has thin disk stars, while the majority of disk stars in M31 belong to a thick disk or kicked up disk (Dorman et al. 2013; Dalcanton et al. 2015) , with all the analysis pointing towards M31 having a more violent merger history than the MW (e.g. Tanaka et al. 2009; Mackey et al. 2019; Ferguson & Mackey 2016) . Leaman et al. (2017) use data from the literature to compare the age-velocity dispersion relation of eight LG members, including the MW, M31, and M33 to model the evolving interstellar medium (ISM) turbulence as a galaxy experiences heating from various sources. In this section, we focus on the observed velocity dispersion of the three most massive LG members and add measurements for M33 that include a robust sample of individually resolved stars across the inner and outer disk of M33. We also present the first comparison of AD as a function of stellar age of LG members. It is important to note that possible halo stars were not removed from the M31 velocity dispersion or AD analysis. Dorman et al. (2012) find that ∼ 15% of stars in the M31 sample are part of a dynamically hot component. Dorman et al. (2013) characterize the bulge, disk, and halo of M31 using the luminosity function of stars from the PHAT survey, line-of-sight velocities of RGB stars from the SPLASH survey, and surface brightness profiles. Using these criteria, Dorman et al. (2013) find that the hot component in the M31 sample is best described by a disk-like surface brightness profile and lumi-nosity function but spheroid-like kinematics, suggesting that ∼ 30% of the dynamically hot stars belong to a kicked up disk and ∼ 10% of the M31 sample is halo contamination. Since Dorman et al. (2013) did not fit for a both a disk and a kicked up disk along with a halo component, it is possible that the halo contamination is lower.

The MW data presented here consists of nearby F and G dwarf stars in the Solar Neighborhood that were targeted using Hipparcos parallaxes as part of the Geneva-Copenhagen Survey (Nordström, B. et al. 2004) . Nordström, B. et al. (2004) use kinematics to identify thick disk stars and distances to target nearby stars but do not formally remove possible halo stars. The MW halo fraction in the Solar Neighborhood is very small (0.15% Du et al. 2018) , so it is possible the sample contains few halo stars but unlikely that these bias the velocity dispersion measurements significantly.

It was important for this work to remove possible M33 halo stars because M33's halo is centrally concentrated and dense, so there is more overlap between the halo and disk stars' sightlines, giving a larger halo fraction in the RGB stars in M33. For the sake of comparison, though, we show the median velocity dispersion and AD values for M33 if halo stars are not removed. We also acknowledge there is not necessarily a clear distinction between a galaxy's halo and kicked up disk, even if they were formed from different mechanisms and that differentiating between a disk, kicked up disk, and a halo can be arbitrary. Figure 14 shows velocity dispersion as a function of stellar age for the MW solar neighborhood (Nordström, B. et al. 2004; Holmberg, J. et al. 2009 ), the northern disk of M31 (Dorman et al. 2015) , star clusters in M33 (Beasley et al. 2015) , and stars in the disk of M33 (this work). M31 has both a greater velocity dispersion (Nordström, B. et al. 2004; Holmberg, J. et al. 2009; Sysoliatina et al. 2018; Dorman et al. 2015; Budanova et al. 2017 ) and a steeper gradient in the age-velocity dispersion than the MW and M33 (Dorman et al. 2015; Bhattacharya et al. 2019) . The median velocity dispersion values are significantly lower in M33, which is expected, as M33 is a less massive and therefore more fragile disk.

The velocity dispersion reported in this work is not consistent with that measured using star clusters in M33 (Beasley et al. 2015) . The star clusters show a higher magnitude that is consistent with the MW Solar Neighborhood and a gradient in velocity dispersion and stellar age. One potential explanation for the difference is that we remove stars with line-of-sight velocities that are consistent with a dynamically hot component. If we do not remove these stars, we do recover a sloped age velocity dispersion relation and higher magnitudes, although not high enough to match the dispersion of the star clusters. It is possible that a larger fraction of the star clusters than individual stars belong to the hot component. It is unexpected that this work does not see a trend between velocity dispersion and age, as it is predicted and seen in the MW and M31 (Leaman et al. 2017; Beasley et al. 2015; Nordström, B. et al. 2004; Dorman et al. 2015) . This lack of trend is what Leaman et al. (2017) calculate for dwarf spheroids in the LG with masses ∼ 1000 times less than M33, not for a more massive galaxy like M33. Their models, which do not distinguish between galactic components, predict M33 should have a velocity dispersion of ∼ 40 km s −1 for stars with ages of several Gyr, which is also consistent with the observations of M33 star clusters (Beasley et al. 2015 ). If we do not remove likely halo stars, the oldest M33 stars in this work have a velocity dispersion that approaches 40 km s −1 . This demonstrates the effect that contaminate halo stars can have in a measured age-velocity dispersion relation. The velocity dispersion of the individual stars and the clusters were also measured in two different ways in two different samples, which results in different systematic uncertainties in the intrinsic dispersion calculation and could contribute to the discrepancy.

The existence of a substantial kinematically hot halo component in M33's inner disk (Gilbert et al. in press) may suggest internal heating mechanisms, like gas flows from stellar feedback (El-Badry et al. 2016) , could have caused enough disk heating to substantially contribute to a stellar halo, which could also cause stellar clusters and individual stars that were born in the disk to be heated or kicked up into the halo component.

Additionally, because of their respective inclinations, the observed line-of-sight velocity dispersion of M33 includes a larger z component than the observed line-ofsight velocity dispersion of M31. Martig et al. (2014) find that in more massive disks, the vertical component of velocity dispersion is largely dependent on merger history but these trends can be obscured when there are large uncertainties on age estimates. It is possible our wide age bins combined with the large influence of the vertical component in M33's line-of-sight velocity dispersion are obscuring a subtle trend in the age-velocity dispersion relation, but it is unlikely to fully explain the low magnitude of velocity dispersion for the intermediate age and old stars.

We also make the first comparisons of the AD as a function of stellar age observed in M31 and M33. There are no current measurements of AD as a function of stellar age in the MW or other galaxies in the LG. Table  5 compares the AD (with respect to H i) in M31 and M33. Like velocity dispersion, M31 has greater values of AD and a steeper rise in AD as a function of age than observed in M33, suggesting that it has experienced several significant and ongoing heating events throughout its relatively recent history. If possible halo stars are not removed from the M33 dataset, there is a trend between velocity dispersion/AD and stellar age. However, this trend disappears if the likely halo stars are removed. Furthermore, the AD values roughly double when possible halo stars are not removed for the intermediate age and old stars.

There are other measurements of AD from Integralfield-units (IFUs): Martinsson et al. (2013) use IFUs to measure the AD with respect to ionized gas in face-on spiral galaxies and find stars lag behind the ionized gas by ∼ 11 ± 8%, which is consistent with the young stars lagging on average by 20%, the intermediate age stars by 9%, and the old stars by 18%. This amount of lag is similar to studies of AD in other inclined local galaxies (Ciardullo et al. 2004; Herrmann & Ciardullo 2009; Westfall et al. 2007 Westfall et al. , 2011 and in the MW (Ratnatunga & Upgren 1997; Olling, Rob P.; Dehnen 2003) .

In this work, we present the largest stellar spectroscopic survey of M33, the TREX Survey, and present initial analysis of the stellar disk kinematics as a function of stellar age using only individually resolved stars, which is the first of its kind in M33. Below we summarize our main findings and conclusions of the complex, yet low mass, galaxy.

• 1. The TREX survey consists of ∼ 4500 stars with good quality spectra across the entire disk of M33, ranging from several evolutionary stages: massive MS and HeB stars, intermediate mass AGB stars, and low mass RGB stars. This work uses a subset (2561 spectra) of the full survey.

• 2. We find that M33's stellar disk has an average velocity dispersion of ∼ 16 km s −1 , which is significantly lower than what is observed in the disk of MW and M31 (Holmberg, J. et al. 2009; Dorman et al. 2015) and lower than what is measured using star clusters (Beasley et al. 2015) . The average magnitude of AD is on the order of ∼ 16, which is also lower than what is observed in M31.

• 3. Velocity dispersion and AD do not increase with stellar age in the disk, which is unexpected. We highlight the importance of removing potential halo stars when measuring the age-velocity dispersion and age-AD relation.

• 4. This analysis suggests that M33 has experienced a significantly different dynamical heating history than M31 and the MW, which may have been dominated by internal heating mechanisms rather than external ones.

• 5. The young stars are as dynamically hot as the older stars in the stellar disk component. These young stars also have a wider distribution of AD values than the old disk stars. They are more dynamically hot than simulated M33-like analogs predict. Possible mechanisms for the heating of the young stars could include turbulence from ongoing star formation or scattering from giant molecular clouds. It is also possible that these stars are remnants of a relatively recent minor merger or other type of galaxy interaction and lie in front of the disk.

A. SPECIAL STARS: WEAK CN AND CARBON STARS As mentioned in Section 3 we identified two classes of rare M33 stars with unusual spectroscopic features, carbon stars and "weak CN" stars in our TREX survey dataset, using the same combination of visual inspection and machine classification of spectra that our team has used to identify such stars in M31 (Hamren et al. 2015 (Hamren et al. , 2016 Guhathakurta et al. 2017; Kamath et al. 2017; Masegian et al. 2019; Kotha et al. 2020; Girardi et al. 2020; Chauhan et al. 2021; Rodriguez et al. 2021) . The carbon stars are identified by their strong CN, CH, and C 2 molecular absorption bands, the strongest of which is the "W"-shaped CN band at 8000Å. By contrast, the "weak CN" spectrum displays a much weaker 8000Å CN absorption feature while the rest of the spectrum resembles that of a normal (i.e., O-rich) star of comparable color in that it displays the near-infrared Ca triplet and, for redder stars, a set of TiO bands.

We have compared the positions of these rare M33 stars in CMDs to theoretical stellar tracks from the PARSEC library (Bressan et al. 1993 (Bressan et al. , 2012 , for a variety of combinations of the PHATTER filters (all except F275W). As expected, this comparison indicates that carbon stars are intermediate mass/intermediate age AGB stars with masses of M ∼ 2-4 M and ages of a few Gyr. In contrast, the "weak CN" stars are massive (M ∼ 4-M ) and short lived (t < 100 Myr) in the RGB/AGB evolutionary phase.

We identify 97 carbon stars and 363 weak CN stars that pass our quality and disk star cuts. For the bulk of this analysis, we have assigned the weak CN stars to the young star group and the carbon stars to the intermediate age group. In this section, we examine the kinematics of these carbon and "weak CN" stars separately to see if they are unique beyond their spectral features. Figure 15 shows the location and line-of-sight velocity of the weak CN and carbon stars. As there is a lack of young stars in the outer disk in our sample due to different target selection epochs, the weak CN stars are centrally concentrated while the carbon stars are present at extended radii. We have calculated the velocity dispersion for these two groups ( Table 6 ). The weak CN stars have a lower velocity dispersion (11.9 km s −1 ) than the young stars, intermediate stars, and old stars (15.9 km s −1 , 15.2 km s −1 , and 16.5 km s −1 , respectively; Table 2 ), and the carbon stars have a slightly higher velocity dispersion (20.3 km s −1 ) than all the three age groups.

We also construct rotation curves and calculate AD for the weak CN and carbon stars. Because most of the carbon stars are located at large radii, we only calculate AD with respect to the H i. The median AD values, width of the distributions, and percentage of outliers in the distribution are shown in Table 6 . The AD of the weak CN stars (8.0 km s −1 ) most closely resembles that of the intermediate age stars (8.0 km s −1 ). The carbon stars' median AD value (22.4 km s −1 ) is more consistent with the old (24.5 km s −1 ) than the intermediate age stars, even though they 58.8 0.0 Table 6 . Median values of velocity dispersion and AD for the weak CN and carbon stars with respect to H i. Also shown is the width (sigma) of the AD distributions and the percentage of outliers in the distributions. The errors on the median value represent the difference between the 16 th and 84 th percentiles divided by √ N, where N is the number of stars.

are believed to be intermediate age. The width of the distribution of AD for the weak CN and carbon stars are both closest to the width of the intermediate age stars (67.7 km s −1 ), which is wider than the old stars. There does not appear to be anything kinematically unique about these rare spectral types.

New Aspects of Magellanic Cloud Research

American Astronomical Society meeting 237

Instrument Design and Performance for Optical/Infrared Ground-based Telescopes

Substructure and Tidal Streams in the Andromeda Galaxy and its Satellites

American Astronomical Society meeting 235

The authors would like to thank everyone who contributed to the general public's safety, health, and well being during the ongoing COVID-19 pandemic. A very large thank you to all of the essential workers and med-ical professionals. Thank you also to the folks at Keck who worked to make PJ observing happen.The authors recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.Thank you to Bruno Villasenor, without whom, our code would have been too slow. ACNQ and PG are grateful to the Sulphurous AA Overspray project for their generous support.We would like to thank the anonymous referee for improving the clarity of this paper. Support for this work was provided by NSF grants AST-1909066 (KMG), AST-1909759 (PG), DGE-1842400 (ACNQ), GO-14610 (BW), and ANID/FONDECYT Regular Project 1210992 (LC). The analysis pipeline used to reduce the DEIMOS data was developed at UC Berkeley with support from NSF grant AST-0071048.

Keck Data Availability Statement: The data used in this article is from the PAndAS catalogue, which is available to download at: http://www.cadc-ccda.hia-iha.nrccnrc.gc.ca/en/community/pandas/query.html.The other data in this work come from sources available on the Keck and HST archive or are not publicly available.