Rotation and Lithium Confirmation of a 500 Parsec Halo for the Open Cluster NGC 2516

Recent analyses of the Gaia data have identified diffuse stellar populations surrounding nearby open clusters. It is important to verify that these"halos","tails", and"strings"are of similar ages and compositions as stars in the denser part of the cluster. We present an analysis of NGC 2516 ($\approx$150 Myr), which has a classical tidal radius of 10 pc and an apparent halo of stars spanning 500 pc ($20^\circ$ on-sky). Combining photometry from Gaia, rotation periods from TESS, and lithium measurements from Gaia-ESO and GALAH, we find that the halo of NGC 2516 is the same age as the cluster's core. Two thirds of kinematically selected halo members out to 250 pc from the cluster center have rotation periods consistent with a gyrochronological age of 150 Myr. A comparison sample of field stars shows no such trend. The lithium abundances of stars in the halo are higher than in the field, and are correlated with the stellar rotation rate and binarity fraction, as has been noted in other young open clusters. Broadly speaking, this work supports a new paradigm wherein the halos of open clusters are often more populous than their cores. We highlight implications for spectroscopic survey targeting, open cluster dispersal, and planet searches around young stars.


INTRODUCTION
Star clusters form when dense filaments in hierarchically structured molecular clouds collide and collapse (Shu et al. 1987). Over the first 10 Myr, feedback effects including protostellar outflows, photoionization, radiation pressure, and supernova shocks disperse the gas (Krumholz et al. 2019).
Since the gas represents about 99% of the mass of the original cloud, gas dispersal enables the majority (∼90%) of stars in the cluster to escape the cluster's gravitational well (Lada & Lada 2003).
From 10 to 1000 Myr, the cluster remnants that survive gas dispersal suffer an onslaught of supernovae, AGB winds, and close stellar encounters that almost always leads to dissolution (Lamers et al. 2010). Extrinsic to the cluster, collisions with giant molecular clouds (Spitzer 1958), and perturbations from the galactic tide in both the radial and vertical dimen-Corresponding author: L. G. Bouma luke@astro.princeton.edu * This is paper 3 of the Cluster Difference Imaging Photometric Survey. sions further promote stellar escape (e.g., Fukushige & Heggie 2000;Bergond et al. 2001).
Finding the stars that have dispersed from their clusters into the galactic field is a pressing topic for a few reasons. One is to understand the spatial extent and duration of star formation events (e.g., Wright & Mamajek 2018). Another is to understand the Sun's birth environment (Adams 2010). The plausibility of finding the stars that formed with the Sun, whether through chemical tagging (Freeman & Bland-Hawthorn 2002;Hogg et al. 2016;Ness et al. 2018), kinematic analyses, or both, can be tested by first using the remnants of open clusters that have only recently dissipated into the field.
A separate project that benefits from the new discoveries of dispersed stellar populations is that of finding young transiting exoplanets (e.g., Mann et al. 2016;Ciardi et al. 2018;David et al. 2019;Livingston et al. 2019;Bouma et al. 2020;Rizzuto et al. 2020;Plavchan et al. 2020;Newton et al. 2021;Tofflemire et al. 2021;Zhou et al. 2021). Young transiting planets are hard to find because young stars are rare, and mostly reside in the crowded galactic plane (Kharchenko et al. 2013;Piskunov et al. 2018). If the dispersed halos of nearby star clusters could be reliably identified, this could expand the census of nearby young stars by up to a factor of 10, based on the expected fraction of stars thought to be lost during gas dispersal.
Although it has been possible to detect dispersed stellar associations for a long time, Gaia is now enabling major advances (e.g., de Zeeuw et al. 1999;Bergond et al. 2001;Zuckerman & Song 2004;Oh et al. 2017;Cantat-Gaudin et al. 2018;Kounkel et al. 2018;Zari et al. 2018;Kounkel & Covey 2019;Fürnkranz et al. 2019). The populations found by the most recent studies can be summarized as follows. On one end are groups with no discernable cores (e.g., Psc-Eri and µ-Tau, Curtis et al. 2019b;). On the other are hierarchically structured associations with many over and under-densities (e.g. the Sco-Cen and Vela associations, Pecaut & Mamajek 2016;Cantat-Gaudin et al. 2019). Here, we focus on an intermediate regime: lowdensity halos associated with a single overdensity, typically an open cluster that was known before Gaia (see Kounkel & Covey 2019;Kounkel et al. 2020;Meingast et al. 2021). In some cases, these halos may actually be tidal tails, as have been reported in the Hyades Röser et al. 2019), the Ursa Major moving group , and Coma Berenices (Tang et al. 2019;Fürnkranz et al. 2019).
One point of difficulty however is that different algorithms for identifying clusters yield differing levels of sensitivity and precision (Hunt & Reffert 2020). For example, using a Gaussian Mixture Model tautologically yields open clusters that are elliptical (e.g., Wallace 2018). Some unsupervised methods yield dispersed and asymmetric structures with number densities down to 100 times lower than the field around the same regions (e.g., Covey 2019 andMeingast et al. 2021). Since different techniques yield quantitatively and qualitatively different outcomes, the purity of kinematically selected samples is important to verify through independent lines of evidence, including rotation periods and spectroscopy.
As part of a Cluster Difference Imaging Photometric Survey (CDIPS, Bouma et al. 2019), we have been making TESS light curves of stars reported to be members of coeval populations and searching them for planets (Bouma et al. 2020). Our analysis of TESS Sectors 1-13 yielded light curves of 483,407 candidate cluster members in the Southern Ecliptic hemisphere, which are available on MAST. 1 As part of this project, we focus here on a single southern open cluster, and ask: is the cluster halo truly coeval with the core? We chose NGC 2516 for this analysis because it is young (100-200 Myr) and close enough (d ≈ 400 pc) to facilitate measurements of stellar rotation periods using TESS. Healy & McCullough (2020) in fact already used TESS rotation periods to study the stellar inclination distribution of stars in the core of NGC 2516. The cluster's halo however was reported by Kounkel & Covey (2019) to span 1 https://archive.stsci.edu/hlsp/cdips ≈ 20 × 10 × 350 pc, which is a larger volume than has been considered in previous rotation analyses (Irwin et al. 2007;Fritzewski et al. 2020;Healy & McCullough 2020). We want to know: is this halo real or is it just a collection of interlopers, foreground and background stars? To what extent can we use Gaia alone to reliably identify age-dated needles in the haystack of field stars? What are the implications, observationally and theoretically, if we can identify the dispersed halos of open clusters?
To assess whether the stars in the cluster's halo are the same ages as stars in the cluster's core, we apply three different age-dating techniques: isochrones, gyrochrones, and lithium depletion (see Soderblom 2010;Soderblom et al. 2014, for reviews). Gyrochronology has now been empirically calibrated for FGK stars with ages spanning ≈10 Myr to ≈4 Gyr (e.g., Rebull et al. 2018;Curtis et al. 2019aCurtis et al. , 2020Barnes et al. 2016). Lithium equivalent widths can be used for relative age-dating of dwarf stars over similar age ranges, though for a narrower range of spectral types (Sestito & Randich 2005). One complication when deriving relative lithium ages of the K dwarfs comes from the correlation between lithium abundance and stellar rotation rate, which has recently been reviewed by Bouvier (2020), and is a subject that we explore in the latter portions of our analysis.
A brief note on terminology. Low-density stellar associations connected to a dense population (a "core") have been described as being in "halos", "strings", "coronae", "snakes", "outskirts", and "tidal tails" (e.g., Chumak & Rastorguev 2006a;Davenport & Sandquist 2010;Kounkel & Covey 2019;Röser et al. 2019;Tian 2020;Meingast et al. 2021). The latter term implies a particular model for the formation of the dispersed group. "Halo" is model agnostic, but is not an ideal term because it can connote spherical symmetry, which is rarely the case. The halo of NGC 2516 is most accurately described as consisting of a leading and trailing tail. We call it a halo for simplicity.
Section 2 summarizes the astrometric analyses of the Gaia data that led to our interest in NGC 2516. Section 3 measures the age of the cluster's halo and core using Gaia photometry (Section 3.1), TESS gyrochronology (Section 3.2), and lithium equivalent widths (Section 3.3). In Section 4 we discuss the implications of our analysis for open cluster dispersal, stellar spin-down, and the lithium-rotation correlation. Section 5 gives our conclusions.

A DISPERSED HALO AND A CORE?
We selected candidate NGC 2516 members based on those reported in the literature. While some pre-Gaia analyses were available (Jeffries et al. 2001;Kharchenko et al. 2013), the purity and accuracy of the Gaia-derived results are the current state of the art. We therefore adopted what we viewed as the most important Gaia-based samples to compare: those of Cantat-Gaudin et al. (2018), Kounkel & Covey (2019) and Meingast et al. (2021). We refer to these studies in the fol-  The core (black) was analyzed by Cantat-Gaudin et al. (2018) using Gaia DR2, and is coincident with the visual overdensity of stars canonically accepted as "the cluster". The halo (blue) is a concatenation of studies by Kounkel & Covey (2019) and Meingast et al. (2021), which used less restrictive membership assignment algorithms described in Appendix A. The field sample is randomly drawn from an (α, δ, π) cube centered on the cluster. The low volume density of the halo stars makes it difficult to visually distinguish the field and the halo samples. lowing as CG18, KC19, and M21 respectively. A useful visualization of these samples is available online. 2 The Gaia clustering studies each used different selection functions, which yielded different results. CG18 considered stars brighter than G = 18 mag within a few degrees of the center of NGC 2516, and reported 1106 candidate cluster members. KC19 and M21 considered stars up to ≈1 mag fainter, and reported 3003 and 1860 members respectively. The unsupervised clustering techniques that each of these studies applied to the second Gaia data release are discussed in Appendix A, as is the overlap between their resulting membership catalogs. Figure 1 shows the cluster members reported by each study in the space of observed positions, proper motions, and radial velocity when available. In Figure 1, and for the remainder of the study, we describe the CG18 members as the "core" of the cluster, and any non-overlapping KC19 and M21 members as the "halo". This defintion implies that there are 1,106 core members, and 2,192 halo members. The distinction is tautological, in the sense that CG18 did not extend their search for members out to tens of degrees from the cluster center. Nonetheless, the CG18 membership catalog is consistent with that of many earlier investigators (e.g., Jeffries et al. 2001;Kharchenko et al. 2013), and is consistent with the general visual overdensity one sees when viewing NGC 2516 in the sky: it appears to span ≈2 • , not ≈ 20 • . Outside of the core and halo, we also define a set of nearby field stars in the "neighborhood" of NGC 2516. Based on the observed distribution of halo members, we drew these stars randomly from the following intervals of right ascension, declination, and parallax: (2) π [mas] ∈ [1.5, 4.0]. (3) We imposed a magnitude limit of G = 19 mag, and ran the queries using the astroquery.gaia module (Ginsburg et al. 2018). We allowed the number of stars in the comparison sample to exceed that in the cluster sample by a factor of ≈5, to ensure broad sampling of stellar masses and evolutionary states. We also required the comparison sample to not overlap with the cluster sample. The style of visualization given in Figure 1 does not make it clear, in our eyes, why the clustering algorithms have decided to associate certain stars with the halo. Canonically, open clusters are spherical, and span ≈10 pc (e.g., Kharchenko et al. 2013). Inverting the parallaxes in Figure 1 shows that members have been reported from line-of-sight distances ranging from ≈300 to ≈600 pc. Is this structure real? What explains its existence? An initial step in visualizing the structure and kinematics of the candidate cluster members is to consider their Cartesian coordinates (Figure 2). We computed these positions by transforming from (α, δ, π) to galactic (X,Y, Z) assuming the astropy v4.0 coordinate standard (Astropy Collaboration et al. 2018). The Galaxy rotates in the +Ŷ direction. The cluster orbit (gray line) was evaluated by taking the median parameters for core members for which CG18 reported membership probabilities exceeding 70%. We then integrated the orbit using gala and the MWPotential2014 potential (Bovy 2015;. The general shape of the cluster itself seems to include a central overdensity and two tails (the halo). One tail is leading the cluster's orbit, and is angled toward the center of the galaxy when viewed top-down. The other tail is trailing the cluster's orbit, and is pointed away from the center of the galaxy. The elongation of the cluster in both the (X,Y ) and (Z,Y ) planes is correlated with the direction of the LSRcorrected median cluster velocity. Possible explanations for this overall morphology are discussed in Section 4.1.
3. AGE-DATING THE HALO OF NGC 2516 3.1. HR Diagram from Gaia The first check on whether the candidate cluster members share the same age is to analyze the Gaia Hertzsprung-Russell (HR) diagrams. Comparable analyses have previously been performed by CG18, KC19, and M21. Figure 3 shows the HR diagram in the space of absolute Gaia G magnitude as a function of observed G BP − G RP color. Magnitudes are from Gaia EDR3; we performed the Gaia DR2 to EDR3 cross-match using the pre-computed gaiaedr3.dr2_neighbourhood table available at the Gaia archive, and required the closest (proper motion and epoch-corrected) source to be the single match. Spectral types are interpolated from the Pecaut & Mamajek (2013)  and metallicity, and varying mass. The halo stars show a similar sequence, though with greater scatter. One possible explanation for the additional scatter is that the halo is more contaminated by field stars. For instance, ≈ 5 red giants in the halo must be field interlopers, because their isochronally implied ages would be inconsistent with that of the cluster.
Another possible explanation for scatter in the halo's HR diagram is differential reddening across different sightlines. The reported halo spans 20 • on-sky, and varies in position from about b = −12 • to b = −20 • , with the stars closest to the galactic plane also being further from the Sun by up to 300 pc ( Figure 1). An HR diagram binned by galactic latitude did show some minimal evidence for differential reddening, and so we expect that both field star contamination and differential reddening could play a role in the larger scatter of the halo relative to the core. A third possibility that we have not explored is whether the binary fraction could also differ between the different regions of the cluster.
To determine the reddening correction across NGC 2516, we queried the 2MASS DUST service 5 and retrieved the total line-of-sight extinction parameters at the positions of each NGC 2516 member. The mean and standard deviation of the E(B − V) values for the Schlegel et al. (1998) and Schlafly & Finkbeiner (2011) maps were consistent, so we adopted the average from Schlegel et al. (1998): E(B − V) S = 0.206 ± 0.039. Bonifacio et al. (2000) found however that the Schlegel et al. (1998) maps overestimate the reddening values when the color excess exceeds about 0.10 mag. We therefore applied the correction proposed by Bonifacio et al. (2000), and additionally corrected for the mean cluster distance by a factor exp(−|d sin(b)|/h) where d is the average cluster distance, b is the average galactic latitude, and h is the scale height of the galactic disk, assumed to be 125 pc. This yielded our adopted E(B−V) = 0.102±0.019, where the uncertainty reflects the standard deviation of the individual Schlegel et al. (1998) values. Finally, to convert to Gaia colors we used the calibration of Stassun et al. (2019), namely E(G BP − G RP ) = 1.31E(B − V) and A G = 2.72E(B − V). This yielded E(G BP − G RP ) = 0.134 ± 0.025, which is used in the plots. To redden the isochrones, we assumed R V = 3.1, and applied the O'Donnell (1994) SED-dependent extinction law star-by-star through the PARSEC interface.
For the metallicity, we considered a range of super and sub-solar metallicities to fit as much of the locus from 0.5 < G BP − G RP < 1.5 as possible, and settled by eye on the slightly super-solar [M/H] = 0.06 (Cummings 2011). Subsolar metallicities led to model predictions too blue along the main sequence by ≈0.1 mag. Our adopted metallicity agrees with the spectroscopic metallicity from Cummings (2011, Sec 4.4.4), and with the iron abundance determined by the Gaia-ESO team (Baratella et al. 2020). The latter result represented an up-revision from an earlier sub-solar metallicity (Randich et al. 2018). We caution though that other subsolar metallicities have been reported (Bailey et al. 2018), and note that the cluster metallicity and reddening are degenerate; if the reddening is lowered, the implied metallicity will increase. While a detailed redetermination of these parameters from the Gaia photometry is beyond our scope, our adopted values for both the reddening and metallicity are within the range of previously reported values in the literature.
Overall, the data and models agree for masses above ≈0.45 M . Below this mass, the data and models diverge at colors redder than G BP − G RP ≈ 2.2, in the sense that the model isochrones have lower luminosities and bluer colors than the observations. 6 The MIST isochrones showed a comparable disagreement (Choi et al. 2016). Proposed explanations for the discrepancy between the models and observations include starspots and incomplete molecular line lists (e.g., Stauffer et al. 2003;Feiden & Chaboyer 2013;Rajpurohit et al. 2013;Mann et al. 2013;Choi et al. 2016). Ultimately, we adopted the PARSEC isochrones because Chen et al. (2014) implemented an empirical temperature-opacity calibration, which leads the PARSEC isochrones to diverge at slightly lower mass than the MIST models. Regardless, for purposes of age-dating the cluster we focus on the mainsequence turn-off (MSTO), since this is where the models have maximum sensitivity to age. Cummings & Kalirai (2018) presented techniques for mitigating some of the complexities of MSTO age-dating (e.g., sparse turnoffs, stellar rotation, high binarity fractions, and the presence of blue stragglers). Combining a UBV color-color analysis with Gaia DR2 cluster memberships, they found MSTO ages for NGC 2516 ranging from 165 to 195 Myr, depending on the choice of model isochrone (Y 2 , PARSEC, MIST, or SYCLIST).
Our goal is to determine whether the ages of the core and halo are consistent. The left panel of Figure 3 suggests that isochronally, they are consistent: stars past the turnoff in both the halo and core are on the same locus. The right panel of Figure 3 demonstrates the precision with which the claim can be made. The MSTO stars are consistent with the 178 and 316 Myr models, but are bluer than the 100 Myr model. The red-giant-branch (RGB) stars are consistent only with the 178 Myr (log 10 t/yr = 8.25) model. Based on the assumption that the width of the MSTO can be attributed to binary stars, we are most interested in the blue edge (blue stragglers are, however, a concern). The blue edge appears most compatible with the 178 Myr model. These results are consistent with the MSTO age range of 165 to 195 Myr found by Cummings & Kalirai (2018), and we refer to their work for the more precise model-dependent comparison. Based on Table 3 of Cummings & Kalirai (2018), we therefore adopt a MSTO age for NGC 2516 of 175 ± 35 Myr, where our quoted uncertainty is a quadrature sum of the statistical and systematic components from Cummings & Kalirai (2018).

Methodology
We began our rotation analysis by considering all 3,298 candidate members of NGC 2516. For each source, we retrieved all available CDIPS light curves, on a per-sector basis. This yielded 2,205 stars with at least one sector from a CDIPS light curve. Each of these stars is brighter than the CDIPS magnitude limit of G RP = 16 mag. 2,270 of the stars from the NGC 2516 source list have G RP < 16. The difference (65 stars) is due to 35 stars from Meingast et al. (2021) which were not available at the time of the CDIPS reductions, and 30 stars falling on the TESS chip gaps. At the distance of NGC 2516, the G RP = 16 mag cutoff imposed during the CDIPS processing corresponds roughly to (G BP − G RP ) 0 of 2.2, or a spectral type of ≈M2V.
We removed systematic trends common to the light curves on each CCD in each individual sector, and stitched together the resulting light curves before searching for rotationinduced periodicity. Details regarding our detrending approach are presented in Appendix B. After detrending, we proceeded with a few cleaning steps: we masked 0.7 days at the beginning and end of each spacecraft orbit, and ran a sliding standard-deviation rejection window over the light curve, which removed any outlying points within ±3×MAD of the median in each window.
We then measured the rotation period from the resulting light curve using the Lomb-Scargle periodogram implemen-    tation in astrobase (Lomb 1976;Scargle 1982;Bhatti et al. 2018). We used the CDIPS aperture radius that, based on theoretical expectations, was expected to give the optimal balance between light from the target star and the sky background (Sullivan et al. 2015). This typically resulted in a circular aperture radius of either 1 or 1.5 pixels. We recorded the top five periodogram peaks, and their corresponding powers. Finally, as a check on crowding, we recorded the number of stars within the aperture with greater brightness than the target star, and with brightness within 1.25 and 2.5 TESS magnitudes of the target star. The resulting rotation periods and periodogram powers are reported, regardless of detection significance, in Table 1. To clean these measurements, we designed automated cleaning criteria, which yielded sets we denote A and B. For Set A, we considered light curves for which the peak Lomb-Scargle periodogram period was below 15 days, the normalized Lomb-Scargle power exceeded 0.08, and for which no companions with brightness exceeding one-tenth of the target star were in the aperture (i.e., excluding visual binaries with ∆T < 1.25, according to the Gaia DR2 source catalog). These limits were chosen after inspecting the light curves visually, to eliminate low-quality rotation periods while retaining high completeness. For Set A, this yielded 987 stars, out of 1,641 stars that met the initial crowding requirements. For Set B, we additionally required that the Lomb-Scargle period fell below the boundary drawn in Figure 4, which was constructed using the Pleiades polynomial from Curtis et al. (2020). This yielded 892 light curves. Boolean flags for each set are included in Table 1. A set of field stars was also analyzed for comparison. The results are discussed in Appendix C. The rationale behind the additional cut in Set B is that NGC 2516 is already known to have few if any stars above the slow sequence (Fritzewski et al. 2020). Stars above the slow sequence are therefore most likely either not cluster members, or they are cluster members, but with unresolved binary companions contaminating the rotation period measurement (e.g., Stauffer et al. 2016, Section 5.1). For the purpose of assessing which stars are rotationally consistent with being in the cluster, this is a viable filter; for the purpose of exploring whether any rotational outliers might be bonafide cluster members, Set A would be the better cut. Figure 4 shows the resulting rotation periods, with Set A shown on top and Set B on the bottom. While Set A is more complete than Set B, this completeness comes at the expense of greater contamination. To facilitate a comparison against the Pleiades and Praesepe, we used the rotation periods and dereddened colors from Table 5 of Curtis et al. (2020), which drew on data from Rebull et al. (2016) and Douglas et al. (2019) respectively. The first order conclusion is that the Pleiades and NGC 2516 appear gyrochronologically coeval for colors spanning 0.5 < (G BP − G RP ) 0 < 1.2. The implications of this comparison for the age of NGC 2516 are discussed in Section 4.4. Figure 5 shows the same data, with the points colored by indicators of binarity. To assess astrometric binarity, we used the renormalized unit weight error (RUWE) from Gaia EDR3. For plotting purposes, we labeled anything with RUWE exceeding 1.2 as an astrometric binary (11% of the overall cluster sample; see e.g., Belokurov et al. 2020). To assess photometric binarity, we fitted a spline to Figure 3, fixing the nodes by hand. We then labeled any points brighter than 0.3 mag above the cluster sequence as photometric binaries. This included 20% of the overall cluster sample. Many of the resulting candidate binaries, though not all of them, appear below the slow rotation sequence. In Set B, over the color range where the slow and fast sequence are present (0.5 < (G BP − G RP ) 0 < 1.3), there are 523 stars, of which 174 (33%) show signs of binarity. Dividing these 523 stars into "fast" and "slow" sequences by eye, the fraction of stars showing signs of binarity in the slow and fast subsamples are 27% (106/289) and 51% (68/134), respectively. We speculate on possible explanations in Section 4.2.

Comparing Kinematics and Rotation
The TESS rotation periods provide an independent check on the Gaia-derived kinematic cluster memberships. We explored this by cross-matching the stars with detected TESS rotation periods against our original target list of 3,298 Gaia members.
To visualize the results, we again opted for Cartesian coordinates, and supplemented them by calculating the tangential stellar velocities relative to the cluster center. Appendix D   Figure 2) and 2D-tangential velocities for halo members. Stars with detected rotation periods in Set B are shown in blue; halo members for which rotation periods should have been detectable, but were not detected, are shown in orange. Stars in the core are shown as black points. In the top right, the ratio of detected to expected stars with rotation periods in Set B is shown versus 3D separation from the core and 2D tangential velocity difference. Bin widths are chosen to enforce the same number of stars in the denominator (N = 54) per bin. The bottom panels show the innermost and outermost 100 stars with detected rotation periods in position and velocity space (from Set A). discusses the projection effect correction required for calculating the tangential velocities: for a star at position (α, δ), we compare the observed proper motion (µ α , µ δ ) with what the proper motion at the star's position would have been if the star were comoving with the core of NGC 2516. This yields a quantity we denote ∆v * , per M21. We convert these proper motion differences to physical units by multiplying by the measured parallax.
The results are shown in the top two rows of Figure 6. Alternative visualizations in the canoncial space of observed positions, parallaxes, and proper motions are presented in Appendix E. In all these figures, to ensure a fair comparison, we only show stars with 0.5 < (G BP − G RP ) 0 < 1.2 for which our TESS pipeline succeeded in making light curves. In other words, the stars in the base sample are those for which we would have expected to detect a rotation period. This color range is preferred because the the slow sequence becomes less defined for spectral types later than ≈K4V. The bins in the histograms (top-right) were chosen to ensure that the same number of stars were in each bin. Rotation periods consistent with a gyrochronological age of 150 Myr are detected for 68% (471/692) of stars reported to be in the cluster. The average detection rate inside of 25 pc, 72% (290/404), is slightly higher than outside of 25 pc, where the detection rate is 63% (181/288). Within ≈3.5 pc, the period detection rate of 59% (26/44) is anomalously low, due to crowding and saturation in the TESS images.
The lower panels of Figure 6 give a second view of how field star contamination affects the outskirts of the cluster. Of the outermost 100 stars in position space with detected rotation periods in Set A, ≈8 appear as outliers above the slow sequence, compared to just a few for the innermost 100 stars. For the outermost 100 stars in velocity space, ≈13 appear as outliers. After accounting for the stars that were expected to show rotation periods and did not, this suggests field contaminant rates of ≈50% for the outermost cluster members in our initial target list.
Despite this level of contamination, the rotation period measurements show that the halo extends to separations of ≈250 pc in physical space from the cluster core. The most widely separated F2V-K2V stars with rotation periods (blue points in Figure 6) are separated by ≈430 pc. The total length of the structure is therefore 400-500 pc, depending on which members of the halo are chosen as the "tips" on either end. This agrees with the overall structure of the halo reported by KC19, with some minor caveats (orange circles) visible in the top-left panels of Figure 6.
In projection-corrected tangential velocity space, the fraction of stars with rotation period detections remains high out to roughly 5 km s −1 . Meingast et al. (2021) by comparison required a physically motivated cut in tangential velocity space of 1.5 km s −1 . Our results show that at the expense of higher field star contamination rates, bonafide members can be identified even out at higher velocity separations.

Lithium from Gaia-ESO and GALAH
The final approach we took for assessing the youth of the halo population of NGC 2516 was an analysis of the Li I 6708 Å doublet. For NGC 2516, two spectroscopic datasets seemed important: Gaia-ESO (Gilmore et al. 2012) and GALAH (Silva et al. 2015). At the time of analysis, Gaia-ESO DR4 and GALAH DR3 were the most relevant (Randich et al. 2018;Buder et al. 2020). The target selection and results from each survey were as follows.
The Gaia-ESO collaboration chose candidate NGC 2516 members to observe with GIRAFFE and UVES based on previously reported literature members and publicly available photometry. Since the existence of the NGC 2516 halo was not known at the time of target selection, very few halo stars are in the sample. Based on the data, Randich et al. (2018) determined stellar parameters (including lithium equivalent widths and metallicities) for 796 stars that they considered possible NGC 2516 members. Cross-matching these against our kinematic list of 3,298 candidate members by position and imposing a 0. 5 maximum separation limit yielded 459 kinematic members with available spectra, 417 in the core and 42 in the halo of NGC 2516. The lop-sided ratio is due to the Gaia-ESO selection function. An inspection of the effective temperature vs. lithium equivalent width based on Table 2 of Randich et al. (2018) showed the surprising feature that stars with T eff 4000 K often had non-zero equivalent widths, despite the expectation that their lithium should be fully depleted. Whether these measurements were significant was not clear, and so we opted to download the spectra from the ESO archive to re-measure equivalent widths.
The GALAH DR3 target selection is discussed by Buder et al. (2020). The relevant aspects for our analysis are that the survey targeted 12 < V < 14 stars at δ < +10 • in select fields, provided the stars were at least ten degrees from the galactic plane. We identified the candidate NGC 2516 members for which spectra had been obtained by searching the GALAH_DR3_main_allstar_v1 catalog, after excluding stars with the stellar parameter bit flags 1, 2, 3. This excludes spectra with unreliable broadening, low S/N, and unreliable wavelength solutions (see Table 6 of Buder et al. 2020). Of our 3,298 candidate NGC 2516 members, 107 had spectra in GALAH DR3. 51 were in the core, and 56 were in the halo. We downloaded 7 the GALAH DR3 spectra for all 107 entries.
We measured the lithium equivalent widths using the specutils package (Earl et al. 2020). We did this by first shifting the spectra to the stellar rest frame, and then focused in on a 15Å window centered on the 6707.835Å lithium doublet. We continuum normalized the spectra with a third-degree Chebyshev series, excluding any regions that showed absorption lines (e.g., Fe I is present at 6703.58 Å and 6705.10 Å). We proceeded by fitting a Gaussian to the continuum normalized spectra, considering only a ±1 Å window centered on the Li doublet. The equivalent widths were then evaluated by numerically integrating the fitted model over the same window. Our approach therefore includes the Fe 6707.44Å blend in the reported Li equivalent widthswhich leads to systematic overestimates of 10 to 15 mÅ (e.g., Bouvier et al. 2018). Finally, to derive uncertainties on the EWs, we repeated the procedure twenty times, but added noise to the spectra, drawn from a normal distribution with a scale set by the standard deviation of the continuum. The reported uncertainties are then drawn from the 84 th and 16 th percentiles of the resulting EW distribution, with minima imposed at 5 mÅ for FGK stars, and 20 mÅ for M dwarfs. They are included in Table 1. We verified the overall scale of our results by comparing our measurements with those of Randich et al. (2018). Most stars with T eff > 4500 K followed a 1-1 relation. For cooler stars, our EW measurements are systematically lower, as expected given the depletion timescales for stars at these ages (e.g., Soderblom et al. 2014). We visually verified from the spectra that, for the most part, we prefer our measurements. Figure 7 shows the concatenation of the Gaia-ESO and GALAH results, with the points colored according to whether the stars are in the core or the halo of the cluster. The GALAH spectra span a color range of 0 < (G BP − G RP ) 0 < 1.5, due to the V < 14 brightness cutoff of the survey. The overall increase of Li EW from 0 < (G BP − G RP ) 0 < 1.0 is driven by the temperatures in the stellar chromospheres: hotter stars fully ionize Li out of its ground state (e.g., Figure 4 of Soderblom et al. 1993b). The depletion of stars redder than (G BP − G RP ) 0 1.4, i.e., later than K4.5V (0.71 M ), is visible relative to the earlier K dwarfs: these stars have burned their lithium. For comparison, we have also plotted the Li EWs measured by Berger et al. (2018) for planet-hosting stars in the Kepler field, and Pleiades members from Bouvier et al. (2018) (we interpolated from their effective temperatures using the Pecaut & Mamajek 2013 table). For clarity, we have omitted the few upper limits reported by Bouvier et al. (2018). The field star distribution peaks at late F-type stars. The majority of kinematically selected stars in the core and the halo show lithium equivalent widths substantially in excess of these field stars.
We can go one step further, and match our sample of kinematically selected stars with spectra against the TESS rotators. The result is shown in Figure 8. At fixed stellar mass, the rapid rotators have lithium equivalent widths an order of magnitude larger than the slow rotators. This effect is mostly apparent in the K dwarfs. Similar trends were noted in the Pleiades over three decades ago (Butler et al. 1987), and are thought to be caused by differences in lithium abundance (Soderblom et al. 1993b). Alternative explanations, such as differences in line-formation conditions  (e.g., chromospheric temperatures, microturbulent velocities, or the presence of starspots) are incompatible with the available data. The lithium-rotation correlation is discussed further below (Section 4.3).
4. DISCUSSION 4.1. How did the halo form?
Any theory for the structure of NGC 2516, the Psc-Eri stream, and the strings and halos found by KC19 and M21 needs to explain a few observations. First is their aspect ratios. The halo of NGC 2516, for instance, extends over at least 500 pc, despite being only ≈25 pc wide. Second is their axial tilts in the plane of the Milky Way -why their leading and trailing arms tend to conform to the Galaxy's differential rotation (e.g., Figure 2 or Figure 13 of M21). Third is the correlation between the mean LSR-subtracted cluster velocity and elongation axes (blue arrows in Figure 2). Such a correlation has also been noted in Coma Ber by Tang et al. (2019); the other halos may show similar trends.
The most likely explanation seems to be that the observed halo structures are tidal tails (e.g., Chumak & Rastorguev 2006a;Krumholz et al. 2019). The idea is that stars near the cluster's tidal radius escape slowly due to the galactic tide, and subsequently form leading and lagging arms due to differential rotation in the Galaxy. The contraction rate of the cluster can affect this process. Whether the exact contraction rate of NGC 2516 (Healy & McCullough 2020), and the corresponding evaporation rate of the core are sufficient to explain the formation of the halo has yet to be assessed.
A second possibility is that the clusters form in larger and more dispersed star formation complexes: the stars in the halo need not have formed in the same "clump" as those in the cluster core. In other words, the formation environment might be a giant molecular filament, rather than a giant molecular cloud. A sample of such filaments has been collected and analyzed by Zucker et al. (2018): if the current ≈20:1 aspect ratio of NGC 2516 were primordial, its structure would match the aspect ratios of what they term either "elongated dense core complexes", or "bone candidates". A more immediate example is the Orion A cloud. Großschedl et al. (2018) showed using young stellar objects as tracers that Orion A is 90 pc long, and that it has a dense "head" and a lower density "tail".
Any differences in the stellar mass function between the core and halo could be informative, since the tidal tail explanation predicts that the mean tail star should be less massive than the mean core star (e.g., Chumak & Rastorguev 2006a). Figure 9 shows  . Masses and absolute magnitudes for the innermost and outermost stars in NGC 2516. The 3,298 candidate NGC 2516 members were split into three groups according to their distance from the cluster center. The inner and outermost groups, each containing 1,099 stars, are counted in 1 mag bins. The mean star in the innermost sample has a 10% higher mass than the mean star in the outermost sample, consistent with expectations for tidal tails.
star has (M G ) 0 = 7.61 (≈0.63M ). While the difference does appear significant given the number of stars involved, the theoretical implications are not entirely clear, since the "molecular filament" explanation might be able to accommodate mass differences between the core and tail through density gradients across the initial filament. Observational incompleteness at the low-mass end seems unlikely to alter the result, since the innermost and outermost regions are subject to similar selection effects. One possible approach for distinguishing the two explanations could be to observationally establish an age-size correlation. If the ≈100 Myr halos are produced primarily through evaporation and expansion, then the younger clusters should be significantly smaller (Chumak & Rastorguev 2006a,b). A theoretical approach along complementary lines would be to perform a dynamical analysis of whether filaments like Orion A, or those in the Zucker et al. (2018) sample are at all expected to evolve into structures resembling NGC 2516.

Stellar rotation at 100-200 Myr
From their rotation period analysis of NGC 2516, Fritzewski et al. (2020) demonstrated that between 100 and 200 Myr, young clusters (NGC 2516, Blanco 1, Psc-Eri, M 35, M 50, and the Pleiades) all show similar rotation period distributions. The FGK stars show a "slow sequence", a population of "rapid rotators", and in some cases stars in the "void" between the two populations (e.g.,, Psc-Eri and NGC 2516 show more stars in the void than the Pleiades). The late M dwarfs with detected rotation periods tend to rotate rapidly (< 3 days), though there may also be late M dwarfs with slow (> 20 days) rotation periods . Although these late M dwarfs are beyond the magnitude limit of our analysis, they have been explored in depth by Irwin et al. (2007) and Fritzewski et al. (2020). A comparison of our results is shown in Figure 10. The bottom panel of this plot includes, in order of precedence: rotation periods from our analysis, from "Class 1 and 2" Fritzewski et al. (2020) sources, from Table 1 of Healy & McCullough (2020), and from Irwin et al. (2007). The main contribution of our study beyond these previous analyses is our sensitivity across both the core and halo of NGC 2516, and a corresponding expansion in the number of stars that can be analyzed. In the following section, we therefore focus on the origins of the features in the rotation-color diagram, and whether our expanded sample can provide any new insights.
The existence of a slow sequence is the main basis of gyrochronology (e.g., Barnes 2003). While the spindown relation (P rot ∼ t 1/2 ) proposed by Skumanich (1972) is not an accurate enough description, the general idea is correct, and the first-order effect can be understood through magnetic braking (Weber & Davis 1967).
The other features of the rotation-color diagram are indications that magnetic braking is not always the only important effect. The question of why the rapid rotators exist, for instance, requires additional physics. Models that vary coreenvelope decoupling timescales and pre-MS disk lifetimes can reproduce each type of population (Irwin et al. 2007;Gallet & Bouvier 2013, 2015; the question then shifts to why these timescales should differ for stars that are otherwise identical. Another possibility is that proposed by Matt et al. (2015): the saturation of the magnetic dynamo yields dP/dt ∼ P when P > P sat , dP/dt ∼ P −1 for P < P sat , and a maximal spin-down rate in the intermediate regime. This could naturally produce the period bimodality observed at ages of ≈100 Myr (though it fails to explain the widening of the distribution observed at later times Godoy-Rivera et al. 2021). Other significant processes could be external to the star entirely: Qureshi et al. (2018) for instance have reported that mergers between giant planets and the star could explain a non-neglible fraction of the rapid rotators.
A separate hypothesis that we can test using the existing data is the idea that rapid rotation can be explained through binarity. Correlations between rapid rotation and binarity in both young clusters and the field have previously been noted (Meibom et al. 2007;Stauffer et al. 2016;Simonian et al. 2019;Gillen et al. 2020). Related effects could include tidal synchronization for the shortest period binaries, or alternatively pre-MS magnetic locking between the disk and star (e.g., Koenigl 1991;Long et al. 2005).
The lower panels of Figure 4 show our attempt at identifying unresolved binaries in the NGC 2516 rotation sample. Binaries resolved by Gaia were already excluded in our initial definitions of the samples. In NGC 2516 we see that the fraction of stars showing signs of binarity in the slow and fast rotation subsamples are 26% (106/289) and 51% (68/134); the fast rotators show a preference for binarity. This correlation  is in line with the earlier findings noted above. While many rapid rotators in NGC 2516 (66%) have no evidence for being binaries, the cuts we used to select binaries (RUWE > 1.2; photometric excess >0.3 mag) leave open a significant fraction of parameter space. A 0.3 mag flux excess, for instance, corresponds to mass ratios M 2 /M 1 0.70, assuming the luminosity L scales as M 3.5 .
Comparing the explanations of disk-locking against tidal synchronization, the population statistics seem to rule out tidal synchronization. A quarter of the NGC 2516 members are rapidly rotating. In comparison, in the field half of Sunlike stars are binaries, and only ≈9% of these binaries have periods below 100 days (Raghavan et al. 2010). If we assume that all binaries with sub-100 day periods become tidally synchronized, then we can only explain a rapid rotator occur-rence rate of ≈5%. Elevated binarity fractions in pre-MS stars likely do not change this picture (see Section 4.4 of Duchêne & Kraus 2013).
A separate effect is that on the slow rotation sequence, Figure 4 shows that the binary stars appear to be either preferentially redder, or to have faster rotation periods than the single stars. One likely explanation for this could be that the unresolved binaries have a component contributing additional red light to the system, skewing the color measurement of the primary. Whether any physical effects could be at play remains a question for future work.

The lithium-rotation correlation
The lithium-rotation correlation (Figure 8) carries over beyond just equivalent widths, and to the actual abundance of lithium in the photosphere (Soderblom et al. 1993b). Non-LTE corrections do not change the overall relationships between Li abundance and stellar temperature or rotation (Carlsson et al. 1994;Lind et al. 2009). The correlation has been seen in the Pleiades, the Psc-Eri stream, and M 35, among other clusters (Bouvier et al. 2018;Arancibia et al. 2020;Jeffries et al. 2020;Hawkins et al. 2020). In NGC 2516, Jeffries et al. (1998) tentatively observed it as well, when analyzing 24 stars in the core of NGC 2516. Our concatenation of the Gaia-ESO and GALAH-DR3 spectra represents a significant expanasion in volume and color range from this earlier analysis.
What causes the correlation between lithium abundance and rotation at fixed stellar mass? Hints likely lie in shared correlations between lithium abundance, rotation rate, radius inflation, internal mixing, and magnetic field strength (Chabrier et al. 2007;Somers & Stassun 2017;Jeffries et al. 2020). For instance, one interpretation of the lithiumrotation correlation based on pure hydrodynamics is that the rapid rotation leads to less efficient convective penetration ("overshoot"), which lowers the mixing efficiency at the convective-radiative boundary (Baraffe et al. 2017). The magnetic field itself may even be sufficient to inhibit the convection (Ventura et al. 1998). Alternatively, the star's temperature and density profiles could be at fault: rapid rotators have the strongest magnetic fields and the most spotted surfaces. These spots lower the photospheric flux, and drive the stellar radius to expand while lowering the core temperatures, which in turn slows the rate of lithium-destroying proton capture reactions (Feiden & Chaboyer 2013;Somers & Pinsonneault 2015).
None of the internal processes noted above answer the initial question of why some of the G and K dwarfs rotate rapidly to begin with. One clue could be that the rapid rotators tend to be in photometric or astrometric binaries. The majority of such binaries have wide (a 200 au) or intermediate (0.5 -100 au) separations (Raghavan et al. 2010).
It is expected that the circumstellar disks in these binaries have different properties than those of single stars, due to e.g., gap formation in the circumbinary disk, disk truncation, and faster disk clearing (Artymowicz & Lubow 1994;Moe & Kratter 2020). It would therefore be reasonable to assume that disk lifetimes in the binary systems are shorter than in single star systems. Eggenberger et al. (2012) showed that longer disk lifetimes should lead to an extended phase of disk-locking, which leads to a greater amount of differential rotation generated in the radiative zone, more efficient shear mixing, and consequently a lower photospheric lithium abundance. This provides a plausible explanation for the overall trend that the rapid rotators are often in unresolved binaries, and that they are also lithium rich. A careful observational analysis including an assessment of the completeness to different separations and mass ratios of binaries would be useful to clarify whether or not all of the rapid rotators are in binaries. If not, this scenario would be falsified.

The age of NGC 2516
Age from Gyrochronology -Comparing the slow sequence of the Pleiades and NGC 2516 more closely, the top row of Figure 4 shows that the two sequences overlap from 0.5 < (G BP − G RP ) 0 < 1.2. At redder colors, 1.2 < (G BP − G RP ) 0 < 1.7, the dispersion in rotation periods increases. The maximum rotation periods seen in NGC 2516 extend up to ≈11 days, rather than the ≈8.5 day upper limit seen in the Pleiades. This is consistent with NGC 2516 being slightly older than the Pleiades.
Fitting a model to substantiate the claim that NGC 2516 is gyrochronologically older than the Pleiades (e.g., Mamajek & Hillenbrand 2008;Angus et al. 2019;Spada & Lanzafame 2020) is, unfortunately, an exercise in tautology. Gyrochronological models are empirically calibrated against the Pleiades and Praesepe. At (G BP − G RP ) 0 = 1.35, we observe rotation periods in NGC 2516 ranging from 7 to 10 days. The 7 day rotation periods are consistent with the Pleiades; the 10 day rotation periods are not. Fitting formulae from the previously cited studies give ages for a star with {(G BP − G RP ) 0 = 1.35, (B − V) 0 = 1.10, P = 7 day} of 204 Myr, 316 Myr, and 107 Myr respectively. Given a Pleiades age of 125 ± 20 Myr (an average of MSTO and rotation-corrected lithium-depletion results; see e.g. Stauffer et al. 1998, Soderblom et al. 2014, and Cummings & Kalirai 2018, only the Spada & Lanzafame (2020) model has an absolute scale that seems to be well-calibrated at these ages and colors. For an 8-day rotation period of a star with the same color, this model quotes a 193 Myr age; at 9 days, 388 Myr. Given this degree of model uncertainty, we prefer the relative statement that NGC 2516 appears slightly gyrochronologically older than the Pleiades, and much younger than Praesepe. A gyrochronological age estimate for NGC 2516 between 150 and 200 Myr would therefore appear reasonable. Much older, and the issue of why the F and G dwarfs do not spin down becomes pressing. We caution however that the spindown rate of FGK stars between 100 and 300 Myr has not yet been empirically calibrated, and if "stalling" were to occur, it could affect our age assessment (see Curtis et al. 2020).
Age from Lithium Depletion -The lithium depletion boundary for NGC 2516 has to our knowledge not been measured. Given the distance to the cluster, this is not surprising: at ∼150 Myr, the LDB is expected to occur at a stellar mass of ≈0.088 M (Soderblom et al. 2014). This would correspond to (G BP − G RP ) 0 ≈ 4.7, well beyond the limits of available cluster membership lists and spectroscopy.
An easier relative measurement to make of lithium depletion is from the EW-color diagram (Figure 7). The NGC 2516 and Pleiades sequences in this diagram appear indistinguishable. A relative lithium-based age estimate for NGC 2516 would therefore be "Pleiades-age". The uncertainties on this estimate are large because lithium depletion timescales are long for Sun-like stars past 100 Myr. For example, NGC 2516 FGK stars are on average more depleted than stars in IC 2602 (≈ 40 Myr), but only marginally so (Soderblom et al. 2014). On the older end though, Figure 7 does show that relative to Praesepe (≈650 Myr), the lithium EWs of NGC 2516 are significantly enhanced. Given these comparisons, a plausible range of ages based on lithium for NGC 2516 would therefore be between ≈100 and 200 Myr.
Adopted Age -As noted in Section 3.1, a color-color and color-magnitude analysis of the main sequence turnoff (MSTO) by Cummings & Kalirai (2018) found that NGC 2516 is 40-60 Myr older than the Pleiades. Given a Pleiades age of 125 ± 20 Myr, this yields a MSTO age for NGC 2516 of 175 ± 35 Myr. The gyrochronological age we have argued is likely older than the Pleiades, i.e., in the range of 150 to 200 Myr. The lithium age based on the depletion of the G and K dwarfs is more uncertain, but is consistent with the Pleiades and could be a bit older (≈100 to 200 Myr). Averaging these three different indicators gives an age for NGC 2516 of 167±20 Myr, if we set the uncertainties to match the absolute Pleiades age uncertainty of ±20 Myr (Soderblom et al. 2014). Whether this average age is at all meaningful, we leave the reader to judge.

CONCLUSION
The combination of astrometry from Gaia, photometry from TESS, and spectroscopy from GALAH and Gaia-ESO has clarified a few things about the structure of NGC 2516.
• Over-densities observed in the Gaia 3D positions and 2D velocities revealed a halo of stars spanning ≈500 pc in the plane of the galaxy. The earliest reference to this halo that we can find is by Kounkel & Covey (2019). The halo is more precisely described as a leading and trailing tail (Figure 2), with the leading edge angled toward the galactic center, relative to the cluster's orbit in the Galaxy. This is consistent with the direction and amplitude of the Milky Way's differential rotation (≈-0.2 deg Myr −1 kpc −1 ), and similar halos/tails have been observed around ≈10 other nearby open clusters (Meingast et al. 2021).
• Isochronal, rotational, and lithium dating show that the halo of NGC 2516 is coeval with its core. In short, all the data are consistent with the halo being real. The Gaia EDR3 data show a main-sequence turnoff that suggest an isochronal age of 150 Myr for both the core and halo (Figure 3). The faintest M dwarfs in the core and halo are also brighter than field stars of the same color. Gyrochronally, stars in NGC 2516 spanning 0.5 < (G BP − G RP ) 0 < 1.2 (spectral types F2V-K3V) overlap with the slow sequence of the Pleiades (Figure 4). At redder colors, from 1.2 < (G BP −G RP ) 0 < 1.7 (spectral types K3V-K6V), there is a larger dispersion in rotation period. The upper envelope of the rotation period distribution at these redder colors extends to longer rotation periods than in the Pleiades (11 vs. 8 day rotation periods). This is consistent with NGC 2516 being slightly older than the Pleiades (i.e., ≈150 Myr). The lithium equivalent widths of the kine-matic cluster members (Figure 7) are consistent with this assessment.
• Bonafide NGC 2516 members exist out to ≈250 pc in separation and out to ≈5 km s −1 in tangential velocity separation from the cluster core. Figure 6 shows the overlap between kinematic and rotational cluster members used to make this assessment. The majority of the halo members do however exist within tangential velocity separations of ≈2 km s −1 .
• The field star contamination rate increases at larger physical and velocity separations from the cluster core.
The right and lower panels of Figure 6 quantify the statement: ≈70% of the innermost candidate members have detected rotation periods consistent with a gyrochronological age near that of the Pleiades. At physical separations of ≈200 pc, this fraction decreases to ≈50%. While the latter contamination fraction may seem high, a comparison sample of field stars (Appendix C) has a detection rate three times lower, and a distribution in the period-color plane inconsistent with that of a population with a single age.
• The average star in the outskirts has lower mass (by ≈10%) than the average star in the cluster center. This is shown in Figure 9, and discussed in Section 4.1. While this is consistent with expectations for a "tidal tail" origin of the halo, the possibility of a primordial halo also merits exploration.
• The rapidly rotating K dwarfs show elevated lithium abundances, and elevated binarity fractions. Each trend has been noted in comparable stellar populations (e.g., Soderblom et al. 1993b;Meibom et al. 2007;Jeffries et al. 2020;Gillen et al. 2020). The rough orders of magnitude of the effects are ≈1 dex in lithium abundance for a factor of 10 in rotation period, and a factor of 2 enhancement in binarity fraction for fast vs. slow rotators. Section 4.3 discusses one plausible connection between the two observations: the disk-locking phase may be truncated in binary systems, which could free the stars to spin up as they contract on the pre-MS. This spin up could result in less efficient internal mixing and a longer-lasting presence of photospheric lithium (Eggenberger et al. 2012).
The existence of the halo itself is not particularly surprising. A star with a 1 km s −1 velocity difference from the cluster average will move away from the center at a rate of 1 pc Myr −1 . For a ∼100 Myr cluster, a halo with characteristic size ∼100 pc is therefore expected, given that the typical velocity dispersions of open clusters are of order kilometers per second.
On the other hand, our ability to reliably identify members of these halos is rather surprising. These stars have remained hidden in the background of the Galaxy throughout centuries of modern astronomy. Their discovery represents an important step toward understanding the dispersal of open clusters. Combined with wide-field photometric surveys, the stars themselves enable advances in studies of stellar angular momentum evolution. Combined with wide-field spectroscopic surveys (e.g., Kollmeier et al. 2017), the halos could also improve our understanding of the chemical evolution of low-mass stars. Finally, the halo stars provide an expanded sample of age-dated stars to search for planets (e.g., Newton et al. 2021), which may help in benchmarking our understanding of planetary evolution. ACKNOWLEDGMENTS L.G.B. is grateful to G. Zhou, B. Tofflemire, A. McWilliam, E. Newton, M. Kounkel, A. Kraus, L. Hillenbrand, and K. Hawkins for the discussions on young stars, rotation, and lithium that encouraged this analysis. L.G.B. and J.H. acknowledge support by the TESS GI Program, programs G011103 and G022117, through NASA grants 80NSSC19K0386 and 80NSSC19K1728. L.G.B. was also supported by a Charlotte Elizabeth Procter Fellowship from Princeton University. This study was based in part on observations at Cerro Tololo Inter-American Observatory at NSF's NOIRLab (NOIRLab Prop. ID 2020A-0146; 2020B-0029 PI: Bouma), which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. This paper also includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). Funding for the TESS mission is provided by NASA's Science Mission directorate. We thank the TESS Architects (G. Ricker, R. Vanderspek, D. Latham, S. Seager, J. Jenkins) and the many TESS team members for their efforts to make the mission a continued success.
Software: astrobase (Bhatti et al. 2018), astropy  Figure 11 is a Venn diagram of the three membership catalogs concatenated in this study. 99% of the CG18 sample overlaps with KC19, M21, or both. By comparison, 36% of the KC19 sample, and 15% of the M21 sample do not overlap with any of the other catalogs. The data-Gaia DR2-used by all the studies was the same. What are the differences in methodology that produce these different membership lists?
CG18 applied a procedure that yielded what we colloquially call "the core". Their membership assignment algorithm was to first query a Gaia DR2 cone around the previously reported {α, δ} of the cluster center, and within ±0.5 mas of its previously reported parallax. The outer radius of their cone was the previously reported angular radius of the cluster (0.71 • ; Kharchenko et al. 2013). No proper motion cut was applied. CG18 then applied an unsupervised classification scheme to G < 18 mag stars within this cone (UPMASK; Krone-Martins & Moitinho 2014). The UPMASK algorithm first performs a k-means clustering on the astrometric parameters, {µ α , µ δ , π}. A "veto" step is then applied to assess whether the groups of stars output from the k-means clustering are more concentrated than a uniform distribution, by comparing the total branch lengths of their respective minimum spanning trees. This binary flag is converted to a membership probability by redrawing values for {µ α , µ δ , π} using the reported uncertainties and covariances. The final membership probability for any given star is then the fraction of draws for which the star passes the veto step. In the case of NGC 2516, the resulting radius within which half of the cluster members were found was 0.50 • . We included all CG18 NGC 2516 members with reported membership probability exceeding 10% in our list of candidate cluster members.
KC19 applied a different unsupervised clustering method using the 5D Gaia DR2 positions and proper motions. KC19 began with a sample of ≈2 × 10 7 stars, mostly within ≈ 1 kpc and with G 18 mag. They then applied the hierarchical density-based spatial clustering of applications with noise algorithm (HDBSCAN; Campello et al. 2015, McInnes et al. 2017) on the entire sample, setting the minimum number of stars per cluster to be 40. KC19 then iterated over a few different cutoff parallaxes to improve sensitivity to structures with distances from the Sun ranging between ≈100 pc and ≈1000 pc. They then fitted ages to the resulting clusters using two different methods: isochronal, and a convolutional neural network. Regarding membership proabilities, KC19 reported the binary "member" or "not" from HDBSCAN. We included all KC19 NGC 2516 members in our list of candidate cluster members.
Finally, M21 considered ten known open clusters within the nearest kiloparsec. One of these clusters was NGC 2516. Their initial list of ≈3 × 10 7 stars was similar to that of KC19, mainly requiring well-behaved astrometry and photometry from Gaia. Given the mean positions and velocities of the clusters determined by CG18, M21 then selected all stars with proper-motion correct tangential velocities within 1.5 km s −1 of the cluster mean. They then applied DBSCAN (Ester et al. 1996) on both the observed galactic {X,Y, Z} coordinates, as well as to a separate deconvolved spatial coordinate distribution discussed in their Section 3.3. The results of their procedure, as well as those of CG18 and KC19, are visible at the website of S. Meingast: https://homepage.univie.ac.at/stefan.meingast/coronae.html. Their reported membership proabilities were binary, similar to those of KC19. All of the M21 NGC 2516 members were included in our list of candidate cluster members.  B. CONSIDERATIONS FOR REMOVING SYSTEMATIC VARIABILITY FROM TESS LIGHT CURVES In detrending before searching for stellar rotation signals, our goal was to preserve astrophysical variability, while removing systematic variability. One particular concern for the TESS light curves is systematic variability at the timescale of the 14-day satellite orbit, mostly induced by scattered light from the Earth and Moon.
We opted to use a variant of the principal component analysis (PCA) approach discussed by Bouma et al. (2019). This PCA approach uses a set of "trend stars" selected from across each CCD according using ad-hoc heuristics that on average lead the trend star light curves to be dominated by systematic variability. The resulting principal component vectors, also referred to as the eigenvectors, are rank-ordered by the degree of variance that they predict in the training set of trend stars.
We then posit that any given target star's light curve is described as a linear combination of the eigenvectors. We also considered the inclusion of additional systematic vectors that could affect the light curve, discussed below. These can be treated as additional "features" in the linear model. To determine the coefficients of the linear model after the full set of eigenvectors (plus optional "sytematic" vectors) had been asssembled, we explored two possible methods: ordinary least squares, and ridge regression. Ridge regression is the same as ordinary least squares, except it includes an 2 norm with a regularization coefficient. The regularization coefficient that best applied for any given target light curve was determined using a cross-validation grid search, through sklearn.linear_modelRidgeCV (Pedregosa et al. 2011). Each target light curve was mean-subtracted and normalized by its standard deviation, as were the eigenvectors. The linear problem was then solved numerically, and the light curve was reconstructed by re-adding the original mean, and re-multiplying by the standard deviation to ensure that the variance of the light curve did not change.
We found that the choice of using ordinary least squares versus ridge regression did not significantly affect the resulting light curves. In other words, the inclusion (or lack thereof) of a regularization term did not strongly alter the best-fitting coefficients. In the spirit of keeping it simple, we opted for ordinary least squares. A few other choices seemed to be more important: • To smooth, or to not smooth the eigenvectors. If the systematic trends are smooth in time, the eigenvectors should be as well. They should not contain residuals from e.g., eclipsing binaries in the template set, and they should also not be intrinsically noisier than the target star. If either of these is the case, extra variability can be introduced into the "detrended" light curves. To address this problem, we opted to smooth the eigenvectors using a sliding time-windowed filter (with a "biweight" weight scheme, implemented in wotan by Hippke et al. 2019). One issue with this is that systematic sharp features ("spikes") are no longer captured, and could be present in the detrended light curves. Since they can be filtered out in postprocessing (e.g., using rolling outlier rejection), we prefer this approach to having systematic features being injected by the PCA detrending.
• How many eigenvectors to use. A larger number always leads to greater whitening. In Bouma et al. (2019), we performed a factor analysis cross-validation to determine the number of eigenvectors to use. The typical number adopted based on this analysis was 10 to 15. While this approach should in theory prevent over-fitting, in our experience, for stellar rotation it still often lead to distorts the signals, especially for rotation signals with small amplitudes and periodicities of 3 days. Shorter signals typically are not distorted, since the eigenvectors do not contain the high-frequency content that leads to the distortions. For the present analysis, we therefore set the maximum number of eigenvectors to be 5.
• Which supplementary systematics vectors to use. We considered using the BGV, CCDTEMP, XIC, YIC, and BGV vectors, packaged with the CDIPS light curves. We found that the background value measured in an annulus centered on the aperture, BGV, tended to produce the best independent information from the PCA eigenvectors, and so we adopted it as our only "supplementary" trend vector. We opted to not smooth it, assuming that it would provide direct complement to the smoothed PCA vectors.
After all these considerations, for every target star, we ultimately decorrelated the raw (image-subtracted and backgroundsubtracted) light curve using a linear model with ordinary least squares. The vectors used in the model were the 1 unsmoothed background flux vector, and 5 smoothed PCA eigenvectors. Figure 12 shows fifty light curves randomly drawn from the resulting Set B. A supplementary figure set that enables inspection of each individual star is available at https://lgbouma.com/notes.

C. ROTATION PERIODS FOR COMPARISON FIELD STARS
To verify the significance of the trends seen in Figures 4 and 5, we searched the CDIPS calibration light curves for rotation periods. This light curve database comprises all G RP < 13 stars that fell on TESS silicon. The need for separate field and cluster star magnitude limits is based on our processing capabilities and is discussed by Bouma et al. (2019). Over the southern sky (Sectors 1-13 of TESS), this corresponded to a sample of 9,619,784 stars. Cross-matching these against the 13,843 randomly drawn stars in the neighborhood of NGC 2516 yielded 1,987 unique stars. The magnitude cut of G RP < 13 at the distances of the neighborhood sample corresponds to an extinction-corrected color cutoff of (G BP − G RP ) 0 0.80, or spectral types G1V. This reaches sufficiently far down the slow sequence to enable a comparison against the cluster star sample.
We performed the same light curve stitching and period-search procedure discussed in the text on the field stars. Within 0.4 < (G BP − G RP ) 0 < 0.8, the same requirements for crowding resulted in 656 field stars for which rotation periods could have been detected. Imposing the same Lomb-Scargle power cutoff (LSP > 0.08) and color-period cutoff as that used to define Set B yielded 107 period detections (16.3%). Within the same color range, 164 of 327 kinematically identified candidate cluster members yielded period detections (50.1%). The period-color distributions are also quite different: the candidate cluster members show an overdensity coincident with the expected gyrochronal age, while the field stars show a much broader rotation period distribution ( Figure 13). This confirms our expectation that our period measurement procedure is not biased in any way that could produce the features we observe in NGC 2516.  D. PROJECTION EFFECTS FOR TANGENTIAL VELOCITIES To compare any given star's proper motion against the cluster's mean, one must consider the position of the star on the sky. This is because two stars sharing identical Cartesian velocities will have different velocities when projected on the sky, as was emphasized by Meingast et al. (2021). Figure 14 shows the magnitude of the correction in equatorial tangential velocities. Figure 15 is an alternative visualization of the data shown in Figure 6. The rotation periods in this diagram correspond to "Set B" as described in Section 3.2. One issue with the visualization as displayed however is that the layering hides the non-rotators  Figure 15. Core and halo of NGC 2516 in position and velocity space, cross-matched against the rotators. The plotted cluster members are those with 0.5 < (GBP − GRP)0 < 1.2 that meet the crowding restrictions described in Section 3.2: they should have been sufficiently bright and non-crowded to enable rotation period detection. Stars in Set B are shown in blue; those which should have had rotation periods detected but did not are in orange. We caution that the appearance of fewer non-rotators being present toward the core is due the layering of the plot: quantitatively, stars toward the cluster center do not all show rotation periods in our analysis (see Figure 6). toward the cluster center: we did not detect rotation periods for roughly one in four stars at the cluster center (see upper-right panel in Figure 6).