Revisiting the Membership, Multiplicity, and Age of the Beta Pictoris Moving Group in the Gaia Era

Determining the precise ages of young (tens to a few hundred Myr) kinematic (``moving") groups is important for placing star, protoplanetary disk, and planet observations on an evolutionary timeline. The nearby $\sim$25 Myr-old $\beta$ Pictoris Moving Group (BPMG) is an important benchmark for studying stars and planetary systems at the end of the primordial disk phase. Gaia DR3 astrometry and photometry, combined with ground-based observations and more sophisticated stellar models, permit a systematic re-evaluation of BPMG membership and age. We combined Gaia astrometry with previously published radial velocities to evaluate moving group membership in a Bayesian framework. To minimize the effect of unresolved stellar multiplicity on age estimates, we identified and excluded multi-star systems using Gaia astrometry, ground-based adaptive optics imaging, and multi-epoch radial velocities, as well as literature identifications. We estimated age using isochrone and lithium-depletion-boundary fitting with models that account for the effect of magnetic activity and spots on young, rapidly rotating stars. We find that age estimates are highly model-dependent; Dartmouth magnetic models with ages of 23$\pm$8 Myr and 33$^{+9}_{-11}$ Myr provide best fits to the lithium depletion boundary and Gaia $M_G$ vs. $B_{P}$-$R_{P}$ color-magnitude diagram, respectively, whereas a Dartmouth standard model with an age of 11$^{+4}_{-3}$ Myr provides a best fit to the 2MASS-Gaia $M_{K_S}$ vs. $B_{P}$-$R_{P}$ color-magnitude diagram.


INTRODUCTION
Associations of nearby, kinematically coherent, coeval stars, known as young moving groups (YMGs), are indispensable laboratories for studying early stellar and planetary evolution.Studies of YMG members have informed research on early evolution of stellar interiors, rotation, magnetic activity, and high-energy emission (Kastner & Principe 2022), protoplanetary disk-dissipation timescales (Silverberg et al. 2020;Higashio et al. 2022), and the occurrence of giant planets (Nielsen et al. 2019).Accurately describing these evolutionary sequences relies on precise and accurate ages for these groups, i.e., by fitting model isochrones to stars in a color-magnitude diagram, modeling the depletion of lithium in the lowest-mass members, or back-tracking stars to a common spatial origin based on their ★ Contact e-mail:renaalee@hawaii.edu kinematics.The ages of individual members (and their planets) can then be assigned to the ensemble age of the host YMG1 .In recent years, the advent of all-sky photometric and astrometric surveys, especially by the Gaia mission, has enabled the identification of new members of known YMGs, as well as discovery of new groups (Faherty et al. 2018;Gagné & Faherty 2018;Ujjwal et al. 2020;Lee et al. 2022).
The  Pictoris Moving Group (BPMG), as the youngest, nearest (median distance ∼50 pc) YMG, permits some of the most detailed studies of pre-main sequence stars (pre-MS), including the intrinsically faintest (lowest and substellar-mass) members.Projected angular separations between binaries are larger, facilitating the identification and characterization of stellar and substellar companions to the members.The age of the group, estimated to be ∼25 Myr by lithium depletion boundary (LDB) fitting (Messina et al. 2016), isochrone fitting, and dynamical traceback (Mamajek & Bell 2014), corresponds to an epoch immediately following the dissipation of protoplanetary disks and the initial phase in the evolution of fully-formed stars and planets.The BPMG's eponymous star,  Pictoris, hosts one of the first directly imaged debris disks (Lecavelier Des Etangs et al. 1993), as well as a rare directly imaged multi-planet system (Lagrange et al. 2020).Among other BPMG members are 51 Eridani, hosting a directly imaged Jupiter analog (Macintosh et al. 2015), PSO J318.5-22, a free-floating planet (Liu et al. 2013), and AU Microscopii, hosting an edge-on debris disk and two transiting planets (Kalas et al. 2004;Plavchan et al. 2020;Cale et al. 2021).
The age of the BPMG has been estimated by fitting model isochrones to sequences in color-magnitude diagrams (CMDs), dynamical traceback analyses, and lithium-depletion-boundary (LDB) model fitting.The derived values range between 10 and 50 Myr, with considerable variation in the uncertainties (Fig. 1, Mamajek & Bell 2014).The most recent estimates cluster around 20-25 Myr, placing the BPMG at an intermediate stage between the Upper Scorpius star-forming region (10±3 Myr, Pecaut & Mamajek 2016), and the Carina YMG at 38-56 Myr (Bell et al. 2015;Gaidos et al. 2022).Dynamical traceback studies of BPMG kinematics suggest a single star-forming origin (e.g., Crundall et al. 2019;Miret-Roig et al. 2020), and thus is no evidence for a physical age spread among BPMG members.
Several effects can contribute to systematic error in the estimated ages of moving-groups.Contamination by non-members from the field can increase the scatter in color-magnitude diagrams and cause a group to appear older.Stars in other, older moving groups within the solar neighborhood, which are more likely to have space velocities more similar to BPMG, may also erroneously be identified as members.Non-members contaminating the sample can be identified and removed based on Galactic space motions, but this requires precise parallaxes, proper motions, and radial velocities (RVs).Unresolved binary companions can make members appear more luminous and hence younger than they actually are on CMDs (Sullivan & Kraus 2021), affect the kinematics used for trace-back analysis, or influence the depletion of Li through the effect of a stellar companion on disk dissipation and stellar rotation (Somers & Pinsonneault 2015).Previous studies have highlighted the importance of identifying and excluding binary systems from group age analyses (e.g., Binks & Jeffries 2014;Alonso-Floriano et al. 2015;Messina et al. 2017).Recently, Miret-Roig et al. (2020) used Gaia DR2 astrometry and ground-based RVs to assess the BPMG age, identifying binaries based on the astrometry and RV variability.
An additional source of uncertainty in deriving ages from model comparisons comes from the elevated activity of members of young groups.At young (< 0.5 Gyr) ages, most low-mass (< 0.8 M ⊙ ) stars are rapidly rotating (P < 10 d), and the surface spottedness and magnetic activity of pre-MS stars are observed to be enhanced compared to more quiescent MS stars (Morris 2020).It is not known if any observed internal age gradients within stellar groups are physical, as the mass-dependence of star formation timescales is not well established.Multiple epochs of star formation may contribute to actual scatter in group age estimates (Krumholz et al. 2019), but as mentioned, this case has not been supported by dynamical studies of BPMG (Crundall et al. 2019;Miret-Roig et al. 2020).
An individual star's membership in a YMG is typically as-sessed by comparing its position in Galactic 6-d (3-d position and 3-d space velocity) phase space to that of the established group mean and dispersion.Previous BPMG membership studies have cataloged candidate members (e.g., Gagné & Faherty 2018) using incomplete astrometric (5-D; 3-D position and 2-D motion) information, followed by measurement of the 6th (RV) dimension to confirm membership (e.g., Schneider et al. 2019).Additionally, several studies have used RVs to identify binaries among BPMG members (Alonso-Floriano et al. 2015;Messina et al. 2017;Miret-Roig et al. 2020).However, multi-epoch observations were often not available, meaning that single-lined spectroscopic binaries with variable RVs could not be identified.
The advent of the Gaia mission (Gaia Collaboration et al. 2016) has provided the precise astrometry and, for many sufficiently bright stars, RVs, to calculate space motions and vet members of YMG.Incomplete RV information can be complemented by ground-based surveys such as the Radial Velocity Experiment (RAVE, Steinmetz et al. 2020).Gaia also provides precise and homogeneous optical photometry, i.e. magnitudes in the -band (similar to but much wider than Cousins ) and a  P −  P color constructed from its wavelength-dispersing photometer.These can be used to construct color-magnitude diagrams to compare with model isochrones, and to identify unresolved binaries or contaminated single-star sources via over-luminosity relative to known single stars.Gaia Data Release 3 provides a metric of astrometric error, the Reduced Unit Weight Error (RUWE), relative to a single-star solution, that can be used to identify unresolved binaries (Belokurov et al. 2020;Wood et al. 2021).This can be complemented by high-angular resolution adaptive optics (AO) imaging obtained by searches for substellar companions; these can readily detect stellar companions at separations of ∼10-1000 au.
New, diverse evolutionary models are available to compare with photometry and measurements of Li abundance, in particular those addressing the effects of magnetic activity and starspots on the structure, convection, and luminosity of low-mass pre-MS stars (Feiden 2016;Somers et al. 2020).The contribution of magnetic pressure to the internal gas equation of state and total near-surface pressure leads to inflated predicted radii for a given mass and age relative to non-magnetic models, or, for a given luminosity and effective temperature, to an older model age.
In this work, we combined precision astrometry, photometry, and RVs from Gaia DR3 with ground-based RV measurements and AO imaging to assess and confirm moving group membership, identify and exclude close binaries, and distill a catalog of singlestar members of the BPMG suitable for age analysis by the CMD and LDB methods.Further, we compared our estimate of the multiplicity fraction from this work with those of other young stellar associations.We fit predictions from a suite of evolutionary models designed for spotted, magnetically active stars to precise, homogeneous measurements of these stars to revisit the age of the BPMG, and consider the implications of our new, model-dependent ages.

Catalog Construction
We obtained confirmed and candidate members of the BPMG from the Gagné et al. (2018) 1 in Gagné & Faherty (2018) have been plotted for comparison and context.The sources for ages marked with a "*" do not report errors for their ages.Kiss et al. (2011), Lee & Song (2018), Messina et al. (2017), andShkolnik et al. (2017), yielding an initial catalog of 433 (assumed individual) stars.Binaries are discussed in Sec.2.2.This sample was cross-matched with Gaia DR3 (Gaia Collaboration et al. 2021) and astrometry and photometry were incorporated for 415 stars and RVs for 99 stars.In total, RVs of 238 candidates were gathered, from RAVE DR6 (Guiglion et al. 2020), and previous RV surveys of low-mass stars and YMGs (i.e., Binks & Jeffries 2016;Malo et al. 2014a;Shkolnik et al. 2017;Terrien et al. 2015), in addition to the Gaia RVs.We calculated probabilities of BPMG membership using BANYAN Σ, a multivariate Bayesian classifier (Gagné et al. 2018).BANYAN Σ calculates  space velocities and   galactic positions from coordinates, proper motions, parallaxes, RVs, and their errors, and compares these to the six-dimensional group means and standard deviations of the 27 nearest identified young moving groups.For candidates with multi-epoch RV measurements, we used the median values in the BANYAN inputs.For candidates with an even number of epochs, the mean of the median values were used.Based on the distribution of membership probabilities computed for candidates with complete kinematic data, we chose a threshold probability of 0.9; this selected 148 of 238 candidates with RV measurements.Separately, we identified 63 of the 195 candidates lacking RVs which would have membership probabilities >0.9 for an "optimal" RV that maximizes that probability and is plausible, i.e., is within the distribution of the measured RVs of confirmed members.These are high-priority targets for ongoing RV follow-up but are not further analyzed in this work.
Main sequence stars could be interlopers if their space motions coincide by chance with the BPMG.To identify and filter these, we compared  > 0.9 candidates with the empirical main sequence of Pecaut & Mamajek (2013) in a absolute  magnitude vs.  P −  P color CMD (Fig. 2).We computed the minimum offset in  magnitude of each candidate from the MS.Based on the distribution of offsets, we chose a cutoff of 0.55 mags and removed nine stars within 0.55 mags of the MS locus and  P −  P > 2, corresponding to a mass of < 0.55 ⊙ (Pecaut & Mamajek 2013), as interlopers.We retained stars with  P −  P < 2 since hotter, more massive stars can be close to or on the MS by the age of the BPMG (Krumholz et al. 2019).
We cross-matched BPMG stars with the 2-Micron All-Sky Survey (2MASS) Point Source Catalog (Skrutskie et al. 2006) to obtain   -band photometry in the infrared, and we corrected for extinction and reddening based on the Leike et al. (2020) 3-d ISM dust distribution using Dustmaps (Green 2018).As expected for nearby stars, ISM reddening is small; the mean reddening is ( − ) = 0.0005 mag, and the mean extinction values are   = 0.0003 mag and    = 7.1 × 10 −5 mag.
The 4" resolution of the 2MASS survey means that many binary stars will be unresolved.Unresolved or partially resolved sources can lead to erroneously brighter and somewhat redder photometry that in turn can affect age-dating (Sullivan & Kraus 2021).To minimize this systematic, we corrected the   -band photometry of those binaries that were resolved by Gaia (see Section 2.2) but not resolved by 2MASS (noted in Table 2).We predicted the 2MASS   magnitude of each Gaia-resolved component using the  magnitude of each compnent and a relation between  −   and  P −  P based on the empirical main sequence of Mamajek & Bell (2014) 2 .We then compared the predicted and the measured   mags, and based on the overall distribution of the differences, we flagged systems with a difference >0.25 mag as significantly contaminated by a companion.Based on the contrast of the predicted   magnitudes, we apportioned the measured   -band fluxes between the binary components, and apparent   mags from those fluxes.This salvaged 9 systems for use in the   vs.  P −  P CMD fitting in Sec.3.1.

Binary Identification
We retrieved visual and spectroscopic binary systems from Shkolnik et al. (2017) and spectroscopic and RV binary systems from Miret-Roig et al. (2020).We identified additional binary/multiple systems that appear as Gaia-resolved components in the DR3 astrometric catalog, unresolved systems identified by Gaia astrometric error, systems that are resolved in high-resolution AO imaging, and unresolved systems via multi-epoch RV variability or overluminosity relative to the single-star locus in a CMD.Fig. 3 shows the approximate coverage of these different methods in magnitude contrast-separation space.
Gaia resolved binaries: We identified sources within 30" with similar parallaxes and proper motions; this search radius corresponds to a binary separation of 1500 au at the median distance of BPMG members, well beyond the peak of the M dwarf binary separation distribution at 50 au (Susemiehl & Meyer 2022), or ∼5 au (Winters et al. 2019), and thus encompassing the vast majority of physical companions.For each pair of member and potential companion we computed likelihoods  1 and  2 that the latter is a common proper-motion companion or an interloper from the field, respectively.We then compute a Bayesian probability   /(  +  ), where  ,  are both evaluated as: except that while the sum for   consists only of a contribution from a single putative companion, the sum for   is over all stars from an equivalent search at an uncorrelated (distant) location from the BPMG member.In Eqn. 1, Δ  and Δ  are the difference in parallax and proper motion, respectively,   and   are the standard errors in parallax and proper motion of a stellar pair, added in quadrature, and  orb is the typical proper motion due to orbital motion, included as a "softening" term.This was estimated by considering the peak separation of nearby (< 25 pc) M dwarf binaries, 5 au (Winters et al. 2019), and 2•5 ≈ 30 au orbit.For missions such as Gaia, a five-year baseline would allow detection of ∼6 au yr −1 orbital motion, which at the BPMG's median distance of 50 pc, is on the order of 100 mas yr −1 .Because of orbital inclinations and eccentricities, realistically there would be a decrement in the observed orbital motions on the order of a half.For our computations, we adopted  orb = 10 mas yr −1 , which follows from an orbital motion of ∼0.5 au yr −1 .Wider separations would correspond to smaller orbital motions that are not detectable.
Based on the distribution of  values we adopted a cutoff of  > 0.98.Our Bayesian analysis identified 30 likely multi-star systems; 29 binaries and one triple.We confirmed 14 of these companions as BPMG members based on their astrometry, independent RVs, and BANYAN Σ.The separations and magnitude contrasts of these 14 systems are shown in Fig. 3 as purple squares, labeled "Gaia Resolved."Binary systems resolved by Gaia and 2MASS were included in subsequent analyses; those resolved by Gaia but not by 2MASS were included if the latter were resolved with AO imaging (see below) or corrected for source confusion as described above.
AO imaging: We collected archival AO images of confirmed BPMG members obtained with any of three instruments: NAOS+CONICA (NACO) on the ESO VLT UT1 (Lenzen et al. 2003;Rousset et al. 2003), NIRC2 on Keck-2 (Service et al. 2016), and NIRI on Gemini-North (Hodapp et al. 2003) from the ESO, Keck Observatory, and Gemini Observatory archives, respectively.In each archive, we performed a 10" cone search (the typical field of view) of the coordinates of each member.We retrieved at least one AO image for 117 BPMG members, and we inspected each image for additional sources.In Fig. 3 (left) we plot the 16 AO-resolved binary components as pink triangles.An injection-and-recovery analysis was performed to assess the completeness and false positive rate of our visual inspection, i.e. by adding simulated companions to randomly selected AO images containing a single (real) star.The point spread function of the target star was scaled and used for the simulated (mock) companion.The magnitude contrast of the simulated companion was randomly sampled from a uniform distribution in magnitude, and the separation was sampled from the Winters et al. (2019) orbital distribution of nearby M dwarf systems, assuming a distance of 50 pc.Separations were limited to 1.5" to focus on a region of the parameter space not resolved by Gaia; the recovery and false-positive rates at wider separations is expected to be the same as at 1.5" out to at least ∼5", half the field of view of the NIRC2 images.The results of the data injection recovery are shown in Fig. 3 (right).We performed a kernel density estimation using the Python package sklearn.neighbors.KernelDensity (Pedregosa et al. 2011) on the recovered (green points) and notrecovered (open points) simulated companions to establish a 50% detection efficiency contour in separation-magnitude contrast space (Fig. 3 (right)).From the injection-and-recovery analysis, we found the visual inspection was reliable for detecting companions having up to 4.5 magnitudes of contrast with the primary and separation as small as 40 mas.
Gaia astrometric error: Gaia RUWE is a measure of astrometric error relative to a single-star fit that is analogous to reduced  2 .Stars with RUWE > 1.4 are almost invariably unresolved binaries (A.Kraus, pers. comm., Belokurov et al. 2020).Fitton et al. (2022) found that single stars hosting protoplanetary disks also exhibit elevated RUWE, but BPMG members are known to host only much less substantial debris disks (Lagrange et al. 2010).Twentynine members were identified as Gaia-unresolved binaries based on RUWE.
Radial velocities: Sources with ground-based RVs at more than one epoch were assessed for binarity based on a reduced  2 test.Based on the distribution of values, we chose a criterion of  2  > 20, flagging 6 systems.We show in Fig. 3 an estimated RV detection limit contour based on the separations and magnitude contrast of the resolved BPMG binaries and the average RV amplitude from members with multi-epoch RVs.We computed the expected separation and magnitude contrast of unresolved binary systems assuming a typical primary mass of 0.3 M ⊙ (for M dwarfs) and maximum separation of 1000 au.Additionally, the probability of the null hypothesis (-value) for RV variability are included for many sufficiently bright stars in Gaia DR3 (Chance et al. 2022) and 14 sources for which  < 0.01 were flagged as spectroscopic binaries.
Photometry: Additional candidate binaries were identified by the over-luminosity of stars with respect to the nominal single-star locus for BPMG (described in Sec.2.1).We assumed the stars to be equal-brightness, the most probable configuration among M dwarfs (Duchêne & Kraus 2013) and constructed a binary locus for the BPMG by shifting the single-star locus by -0.75 mag.We quantified the absolute offset in   , Δ  , of individuals from this locus and identified 18 likely photometric binaries based on a cut-off of Δ  < 0.7 mag.
In summary, our membership catalog began with 433 candidate sources compiled from the literature, of which 238 were tested for BPMG membership with full kinematic data.Among the 148 confirmed members, we identified 30 resolved binary systems; 14 from the Gaia proper motion matching and 16 from AO imaging.We identified and excluded 61 unresolved binary systems from subsequent analyses; 29 from Gaia RUWE, 14 from RVs, and 18 from photometry.Since they are not resolved, they are not shown in Fig. 3 (left), though we provide approximate detection limits for these methods.This left 64 single stars and 30 resolved binary systems for subsequent analyses.

Isochrone Fitting
We estimated the age of the BPMG by fitting isochrones from different model suites to the confirmed members that are single or resolved binaries in the two CMDs plotted in Fig. 2. We assessed four combinations of Gaia and 2MASS colors and magnitudes:   vs.  P −  P ,    vs.  P −  P ,   vs. G-K, and    vs. G-K.To quantify the scatter in each locus, we computed the  2 offsets from best-fit 5th order polynomials for each CMD.The final sample for isochrone fitting is plotted in   vs.  P −  P and    vs.  P −  P CMDs (Fig. 2), which are the two best-performing choices of colormagnitude pairs for CMD analysis.We used the pre-MS models of Baraffe et al. (2015), the Dartmouth standard and magnetic models (Feiden 2016), and the SPOTS models of Somers et al. (2020).The Baraffe et al. (2015) models (hereafter, BHAC) are calibrated for pre-MS and MS low-mass stars.The evolutionary calculations are based on standard input physics of stellar interiors (Chabrier & Baraffe 1997), but with updated atmosphere models compared to the predecessor evolutionary models of Baraffe et al. (1998).While the BHAC models are calibrated for pre-MS stars, they do not account for the effects of stronger magnetic fields occurring at the surfaces of young, rapidly rotating stars.We interpolated model values of Gaia and 2MASS magnitudes over the provided masses and ages to produce a finer age model grid (spaced by 1 Myr) than is provided in Baraffe et al. (2015) to better intercompare the precise colors and absolute magnitudes.
Dartmouth stellar evolution models are suitable for modeling low-mass main sequence (Dotter et al. 2008) and pre-MS stars (Feiden 2016), with or without the effects of stellar magnetic fields (Feiden & Chaboyer 2012;Feiden 2016).The magnetic and nonmagnetic models were used to age-date BPMG in Malo et al. (2014b).Their interior structure and evolution calculations adopt standard input physics (e.g., solar-calibrated mixing length, general ideal gas equation-of-state, and MARCS (Gustafsson et al. 2008) model atmospheres) comparable to other modern low-mass stellar models (e.g., BHAC; Dotter et al. 2008), with updates to improve the accuracy of model predictions for very-low-mass stars and the pre-main-sequence phase of stellar evolution (e.g., Malo et al. 2014b;Feiden et al. 2015).Models including magnetic fields account for the impact of magnetism on the gas equation of state and on the efficiency of thermal convection (Feiden & Chaboyer 2012;Feiden 2016).Surface magnetic field strengths are assumed to be in pressure equipartition with the gas near the stellar photosphere (Feiden 2016).The model grid assumes solar metallicity (appropriate for the BPMG) and solar composition (Grevesse et al. 2007) and includes individual mass tracks with 0.1 ≤ / ⊙ ≤ 2.0, and isochrones with ages from 1 Myr -10 Gyr.Synthetic Gaia photometry is calculated using synthetic spectra from MARCS model atmospheres (Gustafsson et al. 2008), 2MASS   passbands (Skrutskie et al. 2006), and revised Gaia DR2 passbands (Evans et al. 2018) assuming the extinction   is zero.
SPOTS evolutionary models (Somers et al. 2020) address the distorting effects of activity and surface spots on the structure of magnetically active stars.The model grid includes evolution tracks and isochrones for 0.1 -1.3 M ⊙ stars with discrete surface spot coverage fractions of f = 0.0, 0.17, 0.34, 0.51, 0.68, and 0.85.The SPOTS isochrones are empirically calibrated from the color transformations of Pecaut & Mamajek (2013).However, Gaia magnitudes and colors are not provided for objects redder than  P −  P ∼2.5 for models and  > 0, where many BPMG members lie.We extend the fits to  P −  P = 3 for  > 0 cases using the model luminosity, gravity, mass, and dual-temperature regime (for hot and cool surface regions) to obtain bolometric corrections (and hence absolute magnitudes) for Gaia G, B  , and R  and 2MASS Ks magnitudes from the YBC database of bolometric corrections (Chen et al. 2019), assuming A  =0 and solar metallicity.The  P −  P colors were generated by taking the difference between the absolute  P and  P magnitudes.
We computed the spot fractions  , which are specified in terms of surface area, to fractions in luminosity  ′ using: (2) The two model temperatures,  hot and  cool from SPOTS were used to separately compute magnitudes and relative fluxes for hot and cool surface regions, and an average flux was computed based on our modified spot fractions.Objects redder than  P −  P =3 are still included in the fits.
To fit isochrones we used the  2 metric (Naylor & Jeffries 2006), a form of 2-D maximum-likelihood fitting. 2 behaves essentially as a  2 statistic, with an added dependence on the model density of stars in color-magnitude space ( in Eqn. 5 of Naylor & Jeffries 2006).We show in Fig. 4 the results of our isochrone fitting analysis for each of the nine models (BHAC, Dartmouth standard, Dartmouth magnetic, and SPOTS models for each of six values of spot fraction), in Gaia  P −  P vs.   and    CMDs.To better visualize the details of model performance, the difference in absolute magnitude vs.  P −  P for each model relative to the 5th degree polynomial fits to the main BPMG locus are plotted in Fig. 4. (See Fig. 2 for an uncluttered look at the sample and polynomial fits).We report the best-fit ages and  2 values from each of the fits in Table 1.

Lithium Depletion Boundary Analysis
We compared observed and model abundance of lithium vs. effective temperature as an independent constraint on the age of the BPMG.Measurements of the equivalent width (EW) as well as upper limits for the strength of the Li I doublet at 6708Å were gathered from the literature (Bowler et al. 2019;Mentuch et al. 2008;Messina et al. 2016;Zuckerman et al. 2001) for single BPMG members.For this analysis, we chose to exclude binary systems due to confusion in deriving accurate abundances for each binary component, as well as the potential effect of a companion on rotational spin-down and Li depletion timescales compared to single stars (Somers & Pinsonneault 2015;Bouvier et al. 2018).EWs were converted to normalized abundances based on the curves of growth of Palla et al. (2007) in the  eff range 3100 -3600 K and Soderblom et al. (1993) in the  eff range 4000 -6500 K.We scaled these nomralized abundances to absolute values assuming A(Li) = 3.30 dex on an A(H)=12 scale, a value based on meteoritic abundances, a proxy for Li in the protosolar value (Asplund et al. 2009).A typical abundance uncertainty was computed from the typical EW error and  eff error added in quadrature.We show the converted Li abundances as a function of Gaia  P −  P , a proxy for  eff , in Fig. 6, along with predicted abundances from each of the nine models.For each of the models, we identified the best-fit models by minimum  2 fit.

Membership and Multiplicity
From an initial list of 238 stars with complete kinematic information we identified 148 as members of BPMG with a Bayesian probability P > 0.9.Fifteen members were already identified in the literature as binaries.30 systems were identified as resolved binary systems in Gaia DR3 and of these, 14 new companions were independently confirmed in BANYAN.16 binary companions were resolved in AO images; 4 of which were identified in the Gaia search.Left: Best-fit isochrones by minimum  2 .For the SPOTS isochrones, f ′ denotes the revised spot fraction (Eq. 1) from YBC corrections.In both panels, we show a best-fit spline to the main BPMG locus in red.Right: Each of the nine best-fit isochrone curves relative to the spline fit to better visualize the regions of the CMD where the models tend to fit or fail.The  2 values for each fit are given in Table 1.
proportion of binary to single systems in this sample, consistent with previous estimates for pre-MS stars in the literature (∼50-58%, Ghez et al. (1997)).The upper bound of the multiplicity fraction fraction assumes the AO-resolved systems within 4" are bound, as is supported statistically (Kraus & Hillenbrand 2008).We also include all over-luminous sources (candidate binaries identified by photometry) in this estimate, as a majority of them were additionally suggested to be binaries by Gaia RUWE or proper motion matching.Excluding these systems would yield a conservative multiplicity lower-limit of 38%.We compare our results to other YMGs and the field in Fig. 5 et al. 2017;Zúñiga-Fernández et al. 2021).This is possibly because the binary analysis in this work was more thorough than those of those other, more distant groups and clusters.

Age Estimates
For our age analysis, we excluded 47 confirmed members systems due to unresolved binarity; 29 of these were based on the Gaia RUWE criterion, 21 based on RV variability, and 18 based on photometry; with some systems flagged by multiple methods.We show the best-fit isochrones resulting from the CMD fits in Fig. 4. The LDB fits are shown in Fig. 6.We tabulate the best-fit values from each model and fit method in Table 1, and these are shown graphically in Fig. 7.For the Dartmouth Standard and Magnetic models, from which the lowest fit residuals for each method were derived, we assessed the concordance or internal agreement among the ages derived from the three fit methods (Gaia-CMD, Gaia-2MASS-CMD, and LDB).We performed an ANOVA analysis to assess the variance between age estimates.To do this, we computed an F-statistic as the ratio of inter-method to intra-method  2 .This analysis should be considered a rough representation of the method concordance, since we do not compute individual ages for every star, but rather a single age fit to the group locus.The measurements for each method are therefore independent of each other, and their variances cannot necessarily be compared in the way an ANOVA intends.The Dartmouth Magnetic models are the best-fitting in both the Gaia   vs.  P −  P CMD and LDB.The method concordance (F=0.34;p=0.97;  ℎ=0.05)suggests, however, that the internal agreement between ages is not significant, and the variance is likely attributed to the differences in the models as well as the scatter in the data.We find a similar result for the Dartmouth Standard model fits (F=0.29;p=0.97;  ℎ=0.05).

DISCUSSION
We have constructed an updated catalog of the Beta Pictoris Moving Group using the most precise available photometry and astrometry and utilizing complete kinematics for membership validation.We further refined this sample by quantitatively identifying older MS interlopers aligned in kinematic space with the BPMG, but exhibiting older ages on a color-magnitude diagram.Our catalog contains a total of 106 single and resolved companions, and 47 unresolved binaries.To date, this is the most thorough vetting of stellar multiplicity in the BPMG sample, and improves upon the previous samples of Gagné & Faherty (2018) and subsequent literature by newly confirming 37 previously identified candidates and rejecting 90 others.The multiplicity fraction of 38-51±7% (excluding/including pho-  1 shown graphically.For CMD fits (circles, squares) the y-axis is  2 scaled on the left, while for LDB fits (diamonds), the y-axis is  2 scaled on the right.The two y-axis scales are not related.
tometric binaries) reported here is in agreement with the literature (e.g., Shan et al. 2017;Zúñiga-Fernández et al. 2021) within errors (Fig. 5).The slightly higher fraction reflects improved completeness of the binary search, as previous works have tended to utilize a subset of the binary identification methods used in this work (i.e., RVs, AO imaging).The binary fractions from Kounkel et al. (2019) and Jaehnig et al. (2017) in Fig. 5 are the ratio of confirmed spectroscopic binaries (SBs) to the total of their sample.Zúñiga-Fernández et al. ( 2021) also report SB fractions, but extrapolate to an expected occurrence rate.Since we have incompletely investigated SBs in this work (SB1s from RVs, and SB2s from literature), our multiplicity fraction is likely an underestimate.This sample refinement appears to affect the age estimates quite variably among the different models.Curiously, with the refined single-star fit sample, one would expect the isochrone-based ages to bias older (Sullivan & Kraus 2021), but the BHAC and Dartmouth Standard models consistently produce younger ages (∼11-20 Myr) compared to the literature (∼15-28 Myr) employing the same models (e.g., Binks & Jeffries 2014;Malo et al. 2014b).One key difference between the current and previous studies is the availability of the latest Gaia parallaxes to compute absolute magnitudes.There is notable internal disagreement between the models used in this study, suggesting that such age determinations are largely model-(and method-) dependent.We show in Fig. 1 the results of the present work in comparison to previous literature values.The new age estimates from the non-standard models still lie within the uncertainty spread of the literature, but tend to lie at the older end of the range.The overall best fits (by minimum  2 or  2 , producing the lowest fit residuals) from each method (the two CMDs and LDB) from the 27 estimates are in boldface in Table 1.
The model age estimate with the lowest overall residual is given by the Dartmouth magnetic models, which yield the lowest  2 value in a CMD.The Dartmouth magnetic models also result in the lowest fit residuals in two of the three fit methods.The Gaia G vs.  P −  P CMD places the BPMG at 33 Myr according to the Gaia CMD, though the best-fit LDB suggests 23 Myr.We assign an uncertainty on the BPMG age due to model selection of about ±10 Myr (see Fig. 1).The overall best-fit Gaia-2MASS isochrone is from the Dartmouth standard model and suggests a much younger age of 11 Myr, but the magnetic model returns a slightly lower  2 value.These two isochrone curves (Dartmouth magnetic 33 Myr and standard 11 Myr) are nearly identical, highlighted by the similarity in their offsets from the BPMG locus in the upper right panel of Fig. 4 (teal and cyan curves).The magnetic model should be more appropriate, as it should better approximate the actual physics in these highly active stars (Morris 2020).
One caveat of the isochrone fitting is that not all models are appropriate for all stars.We fit the isochrone models for BPMG members spanning a large mass range for homogeneity, however, it is expected that different models are better suited for age determination in different mass ranges.While the lower mass members should be better described by the magnetic and spot fraction models, the higher mass and hotter stars are expected to have fewer spots (Morris 2020).
In the Gaia   vs. B P -R P CMD, the standard models (BHAC and Dartmouth standard) produce the shape of the BPMG locus well, but slightly overestimate the luminosity of the faintest objects.The standard models do not reproduce the shape of the LDB observed in BPMG.This, in conjunction with the poorer isochrone fits, suggests that the magnetic and spot fraction models are better suited for age-dating BPMG and similar YMGs and young associations.The 2MASS   vs. B P -R P CMD generally produced poorer isochrone fits across the 9 models.And in both CMDs, the models tend to overestimate the luminosity of the bluest stars.
With increasing spot fraction, the SPOTS models tend to better fit the hottest stars while underestimating the luminosity or color of the coolest members.As mentioned in Sec.3.1, the SPOTS models are not well calibrated in this region of the CMD (Somers et al. 2020).It is possible that the YBC-interpolated isochrones used here are neglecting to properly account for extinction and reddening, though this should not be significant for the majority of this sample.Another possibility is that the revised luminosity spot fractions (as opposed to the surface spot fraction) overestimate the true spot fractions of the low mass members, but not the higher mass members.Hence, the spot fractions should not be considered constant with stellar mass (and T eff ), as was assumed when the isochrones were generated.Among the SPOTS models, the overall best-performing is the f ′ = 0.837 model, which suggests and age range of 28-45 Myr.However, very high spot fractions are likely unrealistic since they have not been observed among young stars (Cao et al. 2023), and would imply lower photometric variability than observed.
Fitting to the CMD vs. the LDB returned marginally discrepant ages of 33 +11 −9 and 23 ± 8 Myr, respectively.One explanation for this divergence, compared to, e.g., the consistent CMD and LDB ages in Malo et al. (2014b), is the use of revised solar abundances from Grevesse et al. (2007) in the stellar models over previous estimates (Grevesse & Sauval 1998).It is possible that the Grevesse et al. (2007) abundances used in the current models bias CMD ages older compared to those from (Grevesse & Sauval 1998).However, the precise effects of solar abundance estimates on model ages for young stars, and which abundances are to be adopted, are yet to be determined (e.g., Asplund et al. 2021;Serenelli et al. 2009).Another explanation could be a non-solar metallicity for the BPMG, but the deployment of solar-metallicity models.Viana Almeida et al. (2009) estimated the metallicity of the BPMG as [Fe/H] = −0.01 ± 0.08, but only from a single star.Preliminary tests with the Dartmouth magnetic models (G.Feiden, pers. comm.)show that a super-solar metallicity like that found in nearby young clusters (Brandt & Huang 2015;Dutra-Ferreira et al. 2016) shifts the lithium depletion pattern to higher masses, but lower  eff , at a given age.This could decrease or increase the age estimate depending on the relative sampling of the warmer vs. cooler sides of the boundary (Fig. 6).These models also show that a higher metallicity produces a systematic shift to older ages in a CMD of these fully convective stars.Clearly, a robust estimate for the metallicity of the BPMG is needed and, if necessary, fitting with additional models appropriate to that value.
It is important to also consider the connections between multiplicity, rotation, and age estimates (Messina et al. 2016;Messina 2019).Tidal locking in close binaries can lead to spin-up as angular momentum is transferred from the orbit to rotation, while widely separated companions can enable spin-up by truncating and accelerating the dissipation of the circumstellar disk, halting angular momentum transfer from the star (Rosotti & Clarke 2017).Both cases lead to more rapid rotation, but past studies of the connections between binary separation, disk dissipation, and rotation have been ambiguous or contradictory (e.g., Messina et al. 2016;Allen et al. 2017;Kuruwita et al. 2018).Rotation affects the magnetic and surface activity, complicating isochrone fitting, but also affects interior mixing and lithium depletion (Somers & Pinsonneault 2015;Messina et al. 2016;Bouvier et al. 2018).Messina et al. (2016) employed the same Dartmouth models for age-dating BPMG and reported that the results significantly improve when Li EW is de-trended with rotation.
If an older age for BPMG is borne out, there are several implications for our understanding of disk and planet evolution.If the BPMG is indeed older than previously accepted, this is in agreement with the expected typical disk lifetimes of a few to tens of Myr (Richert et al. 2018), since only debris disks remain in BPMG.Additionally, any wide-orbit, directly detected substellar and planetary companions of BPMG members, e.g.51 Eridani, and other well-studied systems mentioned in Sec. 1, would be more massive than previously determined based on their luminosity (Sullivan & Kraus 2021).
The relationship between binary separation, formation, and pre-MS evolution is not well known; this sample of young pre-MS binaries at varying separations paves the way for a detailed study of how the presence of a companion may affect pre-MS spin-down, especially since much of BPMG has been observed by continuous photometric surveys such as TESS.

CONCLUSIONS
In this study we presented an updated membership census and binary catalog for the Beta Pictoris Moving Group, the nearest and youngest co-moving stellar group in the solar neighborhood.We utilized the latest and most precise space-based astrometry, photometry, and radial velocity data from Gaia DR3 to compute precise membership probabilities and identify binaries contaminating the samples of previous photometric age-dating studies.The Gaia data was supplemented with ground-based RVs from surveys and literature, and photometry from the 2MASS survey.
With the Gaia and RV survey data, we were able to confirm 37 new members of BPMG and reject 90 candidates (Table 3) previously identified in the literature which lacked full kinematic data.We found multiplicity fraction in the range 38-51±7% using multiple methods to identify resolved and unresolved binaries among the confirmed BPMG members.The single-star and resolved companion sample was used in age-dating analyses to minimize the expected age-biasing effects of unidentified binaries in a stellar sample.
We used 9 suites of evolutionary models: the BHAC standard models, Dartmouth standard and magnetic models, and the SPOTS models at 6 discrete spot fractions.We found that the difference in isochrone age between our refined sample and earlier works is smaller than the variation among isochrone ages for the different models used here.From isochrone fitting in two CMDs and from LDB fitting, we found the BPMG age ranges from 10-50 Myr.This result highlights the sensitivity of age estimates to different magnetic and spot fraction models.However, the same standard models employed in previous works still give much younger ages of 10-12 Myr, which indicates that unresolved binaries have not had a major effect on previous determinations of BPMG's age.This study shows that isochrone fitting for young stellar populations is largely sensitive to the models themselves, and that, overall, the surface spot and magnetic models out-perform the previous standards.

Figure 2 .
Figure 2. left:Color-magnitude diagram (CMD) based on Gaia absolute  mangitudes and  P −  P colors of confirmed BPMG members.right: CMD using absolute 2MASS   magnitude in place of   .All confirmed BPMG members are shown as grey dots, while those passing the binary exclusion and photometric criteria and used for age estimation (outlined in black) are used in the isochrone fitting described in Sec.3.1.The solid green lines are the best-fit lines in each diagram (5th degree polynomial) and the dashed green lines are the equal-mass binary loci (displaced by −2.5 log 2 ≈ −0.75 mag from the best-fit loci).The blue dashed line is the equal-mass binary locus constructed from the empirical main sequence (blue dotted line) ofPecaut & Mamajek (2013).

Figure 3 .
Figure3.Left: Separation vs. magnitude contrast for binary components in the BPMG primary sample from the present Gaia cone search astrometry comparison and AO image inspection, as well as literature-reported binaries.Unresolved photometric, spectroscopic, and RUWE-identified binaries are not included since these parameters cannot be unambiguously determined.Also plotted are 50% recovery contours for Gaia RUWE binaries fromWood et al. (2021), Gaia -resolved binaries fromZiegler et al. (2018) and for AO from our injection-and-recovery tests.The estimated RV detection limit, based on an assumed primary mass of 1  ⊙ and average RV amplitude of 20 km s −1 , is shown in light green.Right: Results of the AO injection-and-recovery analysis.Simulated binaries that were not recovered by eye (false negatives) are plotted as open points, and successful recoveries are plotted as filled points.The contour represents a 50% detection efficiency.

Figure 4 .
Figure 4.Gaia (top)  and 2MASS-Gaia (bottom) color-magnitude diagrams of the single star BPMG sample.All BPMG members confirmed in this study are shown as grey dots, while those passing the binary exclusion and photometric criteria are outlined in black.Only those outlined are used in the isochrone fitting.Left: Best-fit isochrones by minimum  2 .For the SPOTS isochrones, f ′ denotes the revised spot fraction (Eq. 1) from YBC corrections.In both panels, we show a best-fit spline to the main BPMG locus in red.Right: Each of the nine best-fit isochrone curves relative to the spline fit to better visualize the regions of the CMD where the models tend to fit or fail.The  2 values for each fit are given in Table1.

Figure 5 .
Figure 5. Binary fractions from this work compared to those from previous BPMG studies and different YMGs.Estimates for BPMG are outlined in black, and the dark grey shaded region is the result from this work.The lower bound represents the conservative binary fraction that does not include photometric binaries, and the upper bound is the fraction including photometric binaries.The Kounkel et al. (2019) (light grey shaded) region is an estimate of from RVs and the double-lined spectroscopic binary (SB2) fraction.Jaehnig et al. (2017) give spectroscopic binary (SB) fractions in star-forming regions (sky blue points).The open and closed points represent the raw observed binary fractions and modeled binary fractions correcting for completeness, respectively.Shan et al. (2017) report multiplicity fractions by association (turquoise points) as well as an overall fraction (turquoise shaded region) based on AO imaging.Zúñiga-Fernández et al. (2021) also report a SB fraction (dark blue points).Susemiehl & Meyer (2022) report a multiplicity fraction for M dwarfs, considering systems with mass ratios from 0.1 to 1, and separations from 0 to ∞ (purple shaded region).
Binks & Jeffries (2016) all Gaia DR2 sources with precise parallaxes within the nearest 100 pc for YMG membership.We screened additional confirmed and candidate members from Binks & Jeffries (2014),Binks & Jeffries (2016), Mamajek & Bell (2014) the age of the BPMG in the literature, i.e.Table 6 of Miret-Roig et al. (2020), updated from Table1ofMamajek & Bell (2014), along with the CMD-and LDB-based ages from this study.The best-fit ages for each of three methods used in this work are circled in red.Estimated ages of other young stellar associations Taurus (TAU), Upper Scorpius (USCO), Carina (CAR), and Columba (COL) from Table Table 2 catalogs confirmed BPMG members in 132 systems, including 64 single stars, 42 resolved binaries, and 47 unresolved binaries.Members of binary systems and those for which the 2MASS Ks magnitudes were recomputed based on Gaia photometry are flagged in Table 2.We found that 38-51% of BPMG systems are binaries, the Baraffe et al. (2015)alized relative to an assumed primordial value () = 3.2) vs. Gaia  P −  P color for BPMG members with 6708Å doublet equivalent width measurements from the literature.Points plotted as triangles represent upper limits.The typical measurement uncertainty is represented as the error bar in the bottom left.For each model, the best (minumum  2 ) fit model LDB is plotted.For the incomplete LDB curves from theBaraffe et al. (2015)(BHAC) models, the fitting sample includes only those within the color ranges covered by the models.
Figure 7.The results in Table

Table 1 .
BPMG age estimates for each of the 9 evolutionary model sets and 3 age-dating methods (2 CMD and 1 LDB).The 1 errors are given in the Age column for the best-fitting models (Dartmouth Standard and Magnetic; see Sec. 4.2).Minimum  2 values are given for CMD best-fit isochrones (seeSec.3.1).Additionally, the minimum  2 value from the best-fit LDB ages.The best-fit age results with lowest overall residuals from each fit method (Gaia CMD, Gaia-2MASS CMD, and LDB) are boldface.

Table 2 .
Confirmed members of the Beta Pictoris Moving Group 1

Table 3 .
Rejected Beta Pictoris Moving Group candidates.The probability of membership to BPMG based on proper motions, parallax, and RV is listed in column log    .:mm:ss.ss)Dec (dd:mm:ss.s)log P