Ultracool dwarf benchmarks with Gaia primaries

MNRAS 470, 4885–4907 (2017) doi:10.1093/mnras/stx1500 Advance Access publication 2017 June 16 Ultracool dwarf benchmarks with Gaia primaries F. Marocco,1‹ D. J. Pinfield,1 N. J. Cook,1,2 M. R. Zapatero Osorio,3 D. Montes,4 J. A. Caballero,3,5 M. C. Gálvez-Ortiz,3 M. Gromadzki,6 H. R. A. Jones,1 R. Kurtev,7,8 R. L. Smart,1,9 Z. Zhang,10,11 A. L. Cabrera Lavers,10,12 D. Garcı́a Álvarez,10,12 Z. X. Qi,13 M. J. Rickard1 and L. Dover1 1 Centre for Astrophysics Research, School of Physics, Astronomy and Mathematics, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK of Science, York University, 4700 Keele Street, Toronto, ON M3J 1P3, Canada 3 Centro de Astrobiologı́a (CSIC-INTA), Carretera de Ajalvir km 4, E-28850 Torrejón de Ardoz, Madrid, Spain 4 Dpto. Astrofı́sica, Facultad de CC. Fı́sicas, Universidad Complutense de Madrid, E-28040 Madrid, Spain 5 Landessternwarte, Zentrum für Astronomie der Universität Heidelberg, Königstuhl 12, D-69117 Heidelberg, Germany 6 Warsaw University Astronomical Observatory, Al. Ujazdowskie 4, PL-00-478 Warszawa, Poland 7 Instituto de Fı́sica y Astronomı́a, Universidad de Valparaı́so, Av. Gran Bretaña 1111, Playa Ancha, Casilla 5030, Valparaı́so, Chile 8 Millennium Institute of Astrophysics, Santiago, Chile 9 Istituto Nazionale di Astrofisica, Osservatorio Astronomico di Torino, Strada Osservatorio 20, I-10025 Pino Torinese, Italy 10 Instituto de Astrofı́sica de Canarias, E-38205 La Laguna, Tenerife, Spain 11 Department Astrofı́sica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain 12 GTC Project Office, E-38205 La Laguna, Tenerife, Spain 13 Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, China 2 Faculty ABSTRACT We explore the potential of Gaia for the field of benchmark ultracool/brown dwarf companions, and present the results of an initial search for metal-rich/metal-poor systems. A simulated population of resolved ultracool dwarf companions to Gaia primary stars is generated and assessed. Of the order of ∼24 000 companions should be identifiable outside of the Galactic plane (|b| > 10 deg) with large-scale ground- and space-based surveys including late M, L, T and Y types. Our simulated companion parameter space covers 0.02 ≤ M/M⊙ ≤ 0.1, 0.1 ≤ age/Gyr ≤ 14 and −2.5 ≤ [Fe/H] ≤ 0.5, with systems required to have a false alarm probability <10−4 , based on projected separation and expected constraints on common distance, common proper motion and/or common radial velocity. Within this bulk population, we identify smaller target subsets of rarer systems whose collective properties still span the full parameter space of the population, as well as systems containing primary stars that are good age calibrators. Our simulation analysis leads to a series of recommendations for candidate selection and observational follow-up that could identify ∼500 diverse Gaia benchmarks. As a test of the veracity of our methodology and simulations, our initial search uses UKIRT Infrared Deep Sky Survey and Sloan Digital Sky Survey to select secondaries, with the parameters of primaries taken from Tycho-2, Radial Velocity Experiment, Large sky Area Multi-Object fibre Spectroscopic Telescope and Tycho–Gaia Astrometric Solution. We identify and follow up 13 new benchmarks. These include M8–L2 companions, with metallicity constraints ranging in quality, but robust in the range −0.39 ≤ [Fe/H] ≤ +0.36, and with projected physical separation in the range 0.6 < s/kau < 76. Going forward, Gaia offers a very high yield of benchmark systems, from which diverse subsamples may be able to calibrate a range of foundational ultracool/sub-stellar theory and observation. Key words: binaries: visual – brown dwarfs – stars: late type. 1 I N T RO D U C T I O N ⋆ E-mail: [email protected] Ultracool dwarfs are a mixture of sub-stellar objects that do not burn hydrogen, and the lowest mass hydrogen fusing stars. While 2017 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society  C Downloaded from https://academic.oup.com/mnras/article/470/4/4885/3869256 by guest on 09 December 2021 Accepted 2017 June 14. Received 2017 June 14; in original form 2017 January 20 4886 F. Marocco et al. MNRAS 470, 4885–4907 (2017) systems. Indeed, to take full advantage of the Gaia benchmark population within a reasonable programme of follow-up study, we aim to identify a subset with a focus on covering the full range of UCD properties (i.e. biased towards outlier properties). This sample should reveal the nature of UCDs extending into rare parameter space, i.e. high and low metallicity, youthful and ancient, and the coolest UCDs. To access this outlier benchmark population, it is crucial to identify systems in a very large volume. Wide companions can be confirmed through an assessment of their false alarm probability, using a variety of parameters. System components should have an approximate common distance, since the orbital separation is much less than the system distance, as well as common proper motion and approximately common radial velocity (RV), since the orbital motion should be small compared to system motion. Previous studies have focused on common distance and proper motion, but across the full Gaia benchmark population we may use a different compliment of parameters. In particular, for more distant systems, proper motion will be smaller and RV may be more useful. Once discovered, benchmark systems need to be characterized via detailed spectroscopic studies of both the primary stars and their sub-stellar companions. While Gaia will provide (in addition to astrometry) RV and atmospheric parameter estimates for the primaries, most UCDs will be too faint to be detected by the ESA satellite. Even those UCDs that are bright enough to be astrometrically observed by Gaia will be too faint for its Radial Velocity Spectrometer. So further study of the UCDs will be needed to determine proper motions, RVs, spectral indices and metallicity/age indicators necessary to fully exploit these benchmarks. In this paper, we explore the full scope of the expected Gaia benchmark population, and then present discoveries from our first selection within a portion of the potential parameter space. Sections 2 and 3 describe a simulation we performed of the local Galactic disc, containing wide UCD companions to Gaia stars as well as a population of field UCDs. We simulate constraints on benchmark candidates (using appropriate limits set by Gaia and available deep large-scale infrared surveys), and calculate false alarm probabilities (based on a range of expected follow-up measurements) thus identifying the full benchmark yield within our simulation. We assess the properties of this population, address a series of pertinent questions and determine how best to optimize a complete/efficient identification of the full population of benchmark systems for which Gaia information will be available. In Sections 4–7, we then outline and present discoveries from our initial search. We have targeted systems where the Gaia primary has metallicity constraints (from the literature), and where the UCD companion is a late M or L dwarf (detected by the UKIRT Infrared Deep Sky Survey, hereafter UKIDSS, or the Sloan Digital Sky Survey, hereafter SDSS) with an on-sky separation ≤3 arcmin from its primary. Conclusions and future work are discussed in Section 8. 2 S I M U L AT I O N Our simulation consists of both a field population of UCDs and a population of wide UCD companions to Gaia primary stars. This two-component population allowed us to simulate the calculation of false alarm probabilities (the likelihood of field UCDs mimicking wide companions by occupying the same observable parameter space), and thus identify simulated benchmark UCDs that we would expect to be able to robustly confirm through a programme of follow-up study. Downloaded from https://academic.oup.com/mnras/article/470/4/4885/3869256 by guest on 09 December 2021 most hydrogen-burning ultracool dwarfs (hereafter UCDs) stabilize on the stellar main sequence after approximately 1 Gyr, their substellar counterparts continuously cool down (since they lack an internal source of energy) and evolve towards later spectral types. Their atmospheric parameters are a strong function of age. The degeneracy between mass and age in the UCD regime does not affect higher mass objects (Burrows et al. 1997). Measuring directly the dynamical mass of a celestial body is possible only if the object is part of a multiple system, or via microlensing events. But so far the census of UCDs with measured dynamical masses is very limited (see e.g. Konopacky et al. 2010; Dupuy, Liu & Ireland 2014; Dupuy et al. 2015). Similarly, age indicators are poorly calibrated and, therefore, scarcely reliable, especially for typical field-star ages (> 1 Gyr). The spectra of UCDs are characterized by strong alkali absorption lines, as well as by broad molecular absorption bands (primarily due to water, hydrides and methane; see e.g. Kirkpatrick 2005). A number of these features have been shown to be sensitive to metallicity and surface gravity (both proxies for age), but the majority of studies have been so far purely qualitative (e.g. Lucas et al. 2001; Bihain et al. 2010; Kirkpatrick et al. 2010), and the quantitative attempts to calibrate these age indicators suffer from large scatter and limited sample size (e.g. Cruz, Kirkpatrick & Burgasser 2009; Allers & Liu 2013) or simply do not extend all the way down through the full UCD regime (e.g. Lépine, Rich & Shara 2007; Zhang et al. 2017). Moreover, the cooling tracks for sub-stellar objects are sensitive to the chemical composition of the photosphere, further complicating the scenario (Burrows et al. 1997). The metallicity influences the total opacity by quenching/enhancing the formation of complex molecules and dust grains, all believed to be key factors in shaping the observed spectra of sub-stellar objects. Although a number of absorption features are known to be sensitive to the total metallicity (e.g. Kirkpatrick et al. 2010; Pinfield et al. 2012), no robust calibration has so far been developed to determine the abundances of sub-stellar objects. A way to achieve more accurate, precise and robust calibrations is to study large samples of benchmark UCD objects for which properties such as mass, age and composition may be determined/constrained in independent ways. Benchmark systems come in a variety of forms (e.g. Pinfield et al. 2006), but here we focus on UCDs as wide companions. Such benchmark UCDs (hereafter ‘benchmarks’) may be easily studied, are expected to be found over a wide range of composition and age (i.e. comparable to wide stellar binary populations), and are sufficiently common to offer large sample sizes out to reasonable distance in the Galactic disc (see Gomes et al. 2013). In general, system age constraints and chemical composition can be inferred from the main-sequence primaries (assuming the most likely scenario that the components formed together). This constrains the atmospheric properties of the UCD companions allowing calibration of their spectroscopic atmospheric parameter indicators. While a benchmark population has previously been found and characterized (see e.g. Day-Jones et al. 2011a; Deacon et al. 2014; Baron et al. 2015; Smith et al. 2015; Kirkpatrick et al. 2016; Gálvez-Ortiz et al. 2017), their number remains limited and the parameter space is therefore largely undersampled. The advent of the European Space Agency (ESA) cornerstone mission Gaia (Gaia Collaboration et al. 2016a) provides the potential to greatly expand the scope/scale of benchmark studies. Combined with the capabilities of deep wide-field infrared surveys optimizing sensitivity to distant UCD companions, Gaia will yield exquisite parallax distances and system property constraints (e.g. Bailer-Jones 2003) for an unprecedented sample of benchmark Gaia benchmarks 2.1 The field population benchmarks; Section 2.3) produces a population of ∼1700 000 UCDs. 2.2 The Gaia benchmark population We generated a model population of benchmark systems by selecting random field stars, and adding one UCD companion around a fraction of them. Primary stars were chosen randomly from GUMS, assuming the fraction of L dwarf companions to main-sequence stars, in the 30–10 000 au separation range, to be 0.33 per cent, as measured by Gomes et al. (2013). Note that Gomes et al. (2013) only measured the fraction of main-sequence stars hosting L dwarf companions, and here our simulation assumes the same system fraction for initially injected companions around all types of primaries, i.e. main-sequence stars, white dwarfs (WDs), giants and sub-giants. The fraction of stars hosting late M, T and Y dwarfs follows from the above normalization coupled with our other simulated characteristics. More details on the simulation of benchmark systems are given in the following subsections. 2.2.1 Primaries: GUMS The primary stars of our benchmark systems are selected from GUMS (Robin et al. 2012). The detailed description of GUMS can be found in Robin et al. (2012), and here we only briefly summarize the relevant facts. GUMS represents a snapshot of what Gaia should be able to see at an arbitrary given epoch. As such, it contains main-sequence stars, giants and sub-giants, WDs, as well as rare objects (Be stars, chemically peculiar stars, Wolf–Rayet stars, etc.), thus providing us with a diverse and reasonably complete sample of potential primaries. The stars were generated from a model based on the Besançon Galaxy Model (hereafter BGM; Robin et al. 2003). Since the BGM produces only single stars, binaries and multiple systems were added in with a probability that increases with the mass of the primary star, and orbital properties following the prescriptions of Arenou (2011), resulting in a fraction of binary systems within 10 pc of 24.4 ± 0.4 per cent (Arenou 2011). Exoplanets are added around dwarf stars, following the probabilities given by Fischer & Valenti (2005) and Sozzetti et al. (2009), and with mass and period distributions from Tabachnik & Tremaine (2002). GUMS does not include brown dwarfs. It is important to note here that GUMS generated stars in several age bins, following a constant formation rate over the 0–10 Gyr range, with the addition of three bursts of star formation at 10, 11 and 14 Gyr representing the bulge, thick disc and spheroid, respectively. We selected GUMS primaries within 1 kpc of the Sun and with |b| > 10 deg, and down to G < 20.7 mag (i.e. the limit for Gaia detection; Gaia Collaboration et al. 2016b; Smart et al. 2017). No other restrictions are placed on the selected objects. The full sample of GUMS potential primaries amounts to ∼ 63 000 000 stars. 2.2.2 Companions For our randomly assigned companion population, we assigned masses following the Chabrier (2005) IMF. Distance, age, metallicity, proper motion and RV values were set using the associated primary stars. Similarly to our simulated field population, Teff , log g, radius, and the absolute Mauna Kea Observatory (MKO), SDSS and WISE magnitudes were calculated for the companions using the BT-Settl models (Baraffe et al. 2003, 2015). MNRAS 470, 4885–4907 (2017) Downloaded from https://academic.oup.com/mnras/article/470/4/4885/3869256 by guest on 09 December 2021 We simulated the UCD field population within a maximum distance of 1 kpc, and over the mass range 0.001 < M/M⊙ < 0.12 (a parameter space that fully encompasses our detectable population of UCD benchmarks; see Sections 3.1.1 and 3.1.4). The overall source density is normalized to 0.0024 pc−3 in the 0.1–0.09 M⊙ mass range, following Deacon, Nelemans & Hambly (2008), and consistent with the values tabulated by Caballero, Burgasser & Klement (2008). We simulated the field population across the whole sky with the exception of low Galactic latitudes (i.e. |b| < 10 deg), since detecting UCDs in the Galactic plane is challenging due to high reddening and confusion (see e.g. Folkes et al. 2012; Kurtev et al. 2017). Each UCD is assigned a mass and an age following the Chabrier (2005) lognormal initial mass function (hereafter IMF) and a constant formation rate. In Peña Ramı́rez et al. (2012), it is shown that the Chabrier (2005) IMF describes very well the σ Orionis observed mass function, except for the very low mass domain (M  0.01 M⊙ ), where the discrepancy becomes increasingly large. However, the difference is at very low masses, where the number of expected detections is low given the observational constraints (see Section 3). The observable properties of the UCD are determined using the latest version of the BT-Settl models (Baraffe et al. 2003 isochrones in the 0.001 < M < 0.01 M⊙ mass regime, and Baraffe et al. 2015 isochrones in the 0.01 < M < 0.12 M⊙ mass regime). The isochrones are interpolated to determine Teff , log g, radius, and the absolute UKIDSS, SDSS and Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) magnitudes. The UCDs are then placed in the Galaxy by generating a set of XYZ Cartesian heliocentric coordinates in the same directions as UVW Galactic space motions (X positive towards the Galactic Centre, Y positive in the direction of Galactic rotation and Z positive towards the North Galactic Pole). We assume a homogeneous distribution in X and Y (similar to previous work; e.g. Deacon & Hambly 2006). Although the nearest spiral arm is located at ∼800 pc (Sagittarius–Carina spiral arm; see Camargo, Bonatto & Bica 2015), the most distant of our simulated benchmark population are actually at ∼550 pc (see Section 3.1.4), and our assumption should thus be reasonable. The distribution in Z follows the density laws adopted by the Gaia Universe Model Snapshot (GUMS; see table 2 in Robin et al. 2012). The XYZ coordinates are then converted to right ascension (α), declination (δ) and distance using standard transformations. We assigned to each UCD the UVW components of its velocity by drawing them from a Gaussian distribution centred on zero. The velocity dispersions (σ U , σ V and σ W , respectively) depend on the age of the UCD and are taken, for consistency, from Robin et al. (2012, table 7). V was corrected for the asymmetric drift, also following Robin et al. (2012). UVW are then converted to proper motion and RV using standard transformations. Apparent magnitudes for our simulated objects were calculated by applying the distance modulus ignoring reddening and extinction, which should be low level since our simulated objects are not at low Galactic latitude and are within the local volume. We included unresolved binaries within our sample by assuming a 30 per cent binary fraction (e.g. Marocco et al. 2015) and that all unresolved binaries are equal mass (a reasonable approximation according to e.g. Burgasser et al. 2007). Limiting this field simulation to J < 19 mag or W2 < 15.95 mag (i.e. the same photometric limits we apply to our simulated 4887 4888 F. Marocco et al. Table 1. Properties of the confirmable Gaia benchmarks (CGBs) generated by our simulation. The first column indicates the type of primary, and the second column indicates the type of companion. In the third column, we present the number of companions generated by our simulation when assuming a sloping separation distribution (minimum of range) and a log-flat separation distribution (maximum of range). The last column indicates the number of companions that could be expected to be unresolved binaries themselves, estimated assuming a binary fraction of 32 per cent (see Section 3.3 for details). Primary Single (min–max) Doubles (min–max) Main-sequence stars M dwarfs L dwarfs T dwarfs 16 462–20 842 2392–2948 121–159 5268–6669 765–943 39–51 Sub-giants M dwarfs L dwarfs T dwarfs 76–74 5–1 0 24–24 2–0 0 WDs M dwarfs L dwarfs T dwarfs 96–132 25–38 0–2 31–42 8–12 0–1 Metal-rich stars ([Fe/H] > 0.2 dex) M dwarfs L dwarfs T dwarfs 539–735 84–104 7–7 172–235 27–34 2–2 Metal-poor stars ([Fe/H] < − 0.3 dex) M dwarfs L dwarfs T dwarfs 1717–2209 169–154 10–10 549–707 54–49 3–3 Young stars (< 500 Myr) M dwarfs L dwarfs T dwarfs 35–52 9–15 0 11–17 3–5 0 Thick disc stars M dwarfs L dwarfs T dwarfs 1251–1635 111–110 8–3 400–523 35–35 3–1 Halo stars M dwarfs L dwarfs T dwarfs 28–37 2–2 1–0 9–12 1–1 0 Any Y dwarfs 0(1)–1(3)a 0(0)–0(1)a a The numbers in parentheses assume an extension to the main simulation reaching W2 ∼ 17 mag through a shift-and-stack approach with multiple NEOWISE scans of the sky. For our main simulation, the frequency of UCD companions was assumed to be flat with the logarithm of projected separation (s). This is reasonably consistent with observations over the range where completeness is high (<10 kau; see e.g. Deacon et al. 2014). Companions are assigned out to s = 50 kau, with limited constraints by previous observations on the wider part of this range (see e.g. Caballero 2009). However, we account for the truncation of wide companions through dynamical interaction (see below), which provides a more physical means of shaping the frequency distribution of the widest benchmark companions. And note that we also carried out a ‘rerun’ simulation with the frequency of UCD companions declining linearly with the log of projected separation, more closely matching observations across the full range (see Section 3.3 and Table 1 ), although it is our main simulation that we discuss in detail in Section 3.3. The position of the UCD companion in the sky, relatively to the primary, is generated assuming a homogeneously distributed position angle. Dynamical interactions between stars are known to cause the disintegration of multiple systems (Weinberg, Shapiro & Wasserman 1987). This is particularly critical in the case of our simulated wide benchmarks. The chance of a system undergoing such disintegration increases as a function of time, and can be estimated, given the MNRAS 470, 4885–4907 (2017) system total mass and age, using the method of Dhital et al. (2010). The average lifetime τ of a binary system is given by τ ≃ 1.212 Mtot Gyr, a (1) where Mtot is the total mass of the binary in units of M⊙ and a is the semi-major axis in pc. We removed from our simulated sample all systems whose age is greater than their expected lifetime. To convert from s to a, we assumed a randomly distributed inclination angle between the true semi-major axis and the simulated projected separation. While this is obviously an approximation (a complete treatment would take into account the full set of orbital parameters), it leads to a median a/s = 1.411, closely matching the 1.40 ratio derived from theoretical considerations by Couteau (1960). Particular care was taken while treating UCD companions to WDs. In that case, the total mass of the system and separation change as the main-sequence progenitor evolves into a WD. Therefore, we first estimated the cooling age of the WD using its Teff and log g (given by GUMS) and the cooling tracks for DA WDs from Tremblay, Bergeron & Gianninas (2011, since all WDs in GUMS are assumed to be DAs; Robin et al. 2012). We estimated the mass of the WD progenitor using the initial-to-final mass relation of Catalán Downloaded from https://academic.oup.com/mnras/article/470/4/4885/3869256 by guest on 09 December 2021 Companion Gaia benchmarks et al. (2008). We assumed that the orbit of a companion around a star that becomes a WD will expand stably, such that aMS aWD = , (2) MWD MMS where MWD and MMS are the mass of the WD and of its progenitor, and aWD and aMS are the semi-major axes of the orbit in the WD and main-sequence phase, respectively. We then calculate a ‘disintegration probability’ for each stage as follows (3) pMS = τstar /τMS , (4) where τ cool is the cooling age of the WD, τ star is the ‘main-sequence age’ (the difference between the total age given by GUMS and τ cool ), τ WD is the expected lifetime in the ‘WD stage’ (i.e. assuming a = aWD and Mtot = MWD + MUCD in equation 1) and τ MS is the expected lifetime in the ‘main-sequence stage’ (i.e. assuming a = aMS and Mtot = MMS + MUCD in equation 1). We removed systems whose total disintegration probability pWD + pMS is greater than one. We note that the initial-to-final mass relation only holds in the mass range 0.5 < MWD /M⊙ < 1.1. For WDs more massive than 1.1 M⊙ , we simply assume the main-sequence lifetime of the progenitor to be negligible compared to the cooling age of the WD, since MMS would be greater than 6 M⊙ . For WDs less massive than 0.5 M⊙ , we assume the mass-loss during the post-main-sequence evolution to be negligible, hence MMS ≃ MWD . 2.3 Simulating constraints on candidate selection and follow-up After generating the field and benchmark populations, we simulated limitations on UCD detection within infrared surveys, as well as the accuracy of observational follow-up. This involves imposing magnitude and minimum separation cuts, and generating realistic uncertainties (typically achieved) on the observables (magnitudes, distance, RV and proper motion). The first step is to set detection limits for our simulated UCDs. Current near- and mid-infrared (hereafter NIR and MIR) surveys probe the sky at different depths and with different levels of multiband coverage, but rather than try to simulate all these different surveys (which would be convoluted, and may change in the future) we took a somewhat simplified approach. In the NIR, we chose a depth limit of J ≤ 19 mag, which can be achieved in a variety of ways. The ongoing Visible and Infrared Survey Telescope for Astronomy (hereafter VISTA) Hemisphere Survey (VHS; McMahon et al. 2013) is scanning the Southern hemisphere down to J = 21.2 mag and Ks = 20.0 mag, allowing for the detection and selection of UCD candidates, e.g. via J − Ks colour criteria. In the Northern hemisphere, the combination of the UKIDSS Large Area Survey (ULAS; with limiting magnitude J = 19.5 mag), SDSS, UKIDSS Hemisphere Survey (with J depth similar to ULAS) and Pan-STARRS 1 (Chambers et al. 2016) will allow the effective selection of UCD candidates (using e.g. z − J criteria) across the full hemisphere. While NIR surveys should be ideal to select most UCDs, some with very red NIR–MIR colours (particularly Y dwarfs) will be best detected in the MIR. WISE is scanning the whole sky down to a 5σ limit of W2 = 15.95 mag. We can therefore expect to identify UCDs down to this limit, by selecting W2-only detections or objects with very red W1 − W2 colours. These objects would be much fainter in the NIR bands; however, the spectroscopic follow-up of WISEselected targets (down to J ∼ 21–22 mag) is routinely achieved with the aid of the latest generation of 6–8 m class telescopes (e.g. Cushing et al. 2011; Kirkpatrick et al. 2012, 2013; Pinfield et al. 2014). Any simulated object (either in the field or part of a benchmark system) fainter than J = 19.0 mag and W2 = 15.95 mag is therefore considered undetectable and removed from our simulated population. Since we are only targeting resolved star+UCD systems, we need to remove all unresolved companions. The angular resolution of existing NIR and MIR surveys is rather patchy, varying from 1–2 arcsec in the best cases (e.g. SDSS, VISTA, UKIDSS) to ∼6 arcsec for WISE. Additionally, large area surveys are known to have issues identifying and cataloguing sources around bright stars, pushing the detection limit for faint companions out to larger separations. We chose to adopt an ‘avoidance radius’ dependent on brightness, i.e. an area of sky around a star where faint UCDs will go undetected. Examination of a range of example stars in SDSS (where this effect is quite clear) led us to set the following values: (i) 15 arcmin for stars with V ≤ 4 mag; (ii) 10 arcmin for stars with 4 < V ≤ 6 mag; (iii) 5 arcmin for stars with 6 < V ≤ 8 mag; (iv) 4 arcsec for stars with V > 8 mag (if J ≤ 19 mag); (v) 20 arcsec for stars with V > 8 mag (if J > 19 mag and W2 ≤ 15.95 mag). Simulating realistic uncertainties on distance, proper motion and RV is a complicated exercise, since it depends not only on the brightness of the UCD, but also on the type of follow-up assumed. For instance, dedicated astrometric campaigns can achieve a high level of precision on parallax and proper motion down to very faint magnitudes (∼1 mas down to J ∼ 20 mag; e.g. Dupuy & Kraus 2013; Smart et al. 2013), but are time-consuming and limited to a relatively small number of objects. However, we take a simplified approach since it is more common to measure proper motion for UCDs using just two epochs, i.e. following up the original discovery images at a later epoch allowing a long enough time baseline. With secondepoch images often obtained with a different telescope/filter, the precision of such measurements is limited. With medium- to highresolution spectroscopy, one can obtain RVs down to a precision of a few km s−1 or less (e.g. Zapatero Osorio et al. 2007; Blake, Charbonneau & White 2010) but these observations are limited to the brighter objects only. Following these considerations, we adopted a 10 mas yr−1 uncertainty for proper motions (assuming a two-epoch measurement and a ∼5 yr baseline), and a precision of 2 km s−1 for RVs down to J = 18 mag (Marocco et al. 2015). For fainter objects, we consider a RV measurement to be currently unfeasible, and we therefore assume their RV to be unconstrained. To simulate distance uncertainties, we considered the current spectrophotometric distance calibrations. Although based on an increasing number of UCDs with measured parallaxes (see e.g. Marocco et al. 2010; Dupuy & Liu 2012), these calibrations are limited by the intrinsic scatter in the UCD population, primarily due to age and composition differences among objects of similar spectral type. The typical scatter around the polynomial spectrophotometric distance relations is ∼0.4 mag (Dupuy & Liu 2012). We thus adopted distance modulus uncertainties of 0.4 mag. No systematic uncertainty is considered for e.g. young or peculiar objects. For unresolved binaries (in our simulated field population), where a spectrophotometric relation would lead to an incorrect distance estimate, we assume their observed distance to be 30 per cent closer than their real distance, with distance modulus uncertainties of 0.4 mag. MNRAS 470, 4885–4907 (2017) Downloaded from https://academic.oup.com/mnras/article/470/4/4885/3869256 by guest on 09 December 2021 pWD = τcool /τWD 4889 4890 F. Marocco et al. 2.4 Companionship probabilities and ‘confirmable Gaia benchmarks’ 3.1 Intrinsic properties of CGBs 3.1.1 Mass, age and spectral type Fig. 1 shows the mass–age distribution of the CGBs. For each of the GUMS age bins, we have introduced a random scatter (across the bin) so as to obtain a more natural continuum of ages and make the plot easier to view. Note that this leads to ‘step-like’ behaviour as one moves across the age bins. The M–L and L– T spectral type transitions can be seen separating the grey–red and red–blue plotting colours. The large majority of CGBs are low-mass stars, with 82 per cent having masses above the sub-stellar limit and 18 per cent being brown dwarfs. The lowest mass CGBs are found in the youngest age bin with masses down to ∼0.02 M⊙ . Most CGBs (87 per cent) are ultracool M types, with about 13 per cent having L or T spectral type. The predominance of late M-type CGBs is due to a combination of three factors: (i) M dwarfs are brighter and therefore can be seen out to larger distance given our adopted magnitude limits; (ii) M dwarfs are intrinsically more numerous given the adopted IMF; (iii) M dwarf companions are generally more massive than L and T dwarfs, and are therefore more likely to survive dynamical disruption (a less significant factor, but not negligible). At the oldest extremes, the large number of M-type CGBs includes 37 halo systems (with ages > 13 Gyr). Of the 2987 L-type CGBs, about two-thirds are young disc, and ∼30 per cent old disc. There are also 110 thick disc L-type CGBs, but a very limited number of halo L types (just 2 simulated CGBs). Most of the 160 T-type CGBs are nearly evenly split between the young and old disc populations (43 and 52 per cent, respectively), with a small but potentially interesting collection of 3 T-type CGBs in the thick disc. Our simulation does not predict any T-type CGBs in the halo. At young ages, there are 67 late M- and L-type objects < 500 Myr, but no T-type objects in this age range. Our most youthful age bin (< 100 Myr) contains 52 M types and 15 L dwarfs. Our main simulation does contain one Y dwarf (with Teff = 490 K). Although WD 0806-661 B (Luhman, Burgasser & Bochanski 2011) is a known wide Y dwarf companion to a WD (discovered in Spitzer data), it lies beyond our all-sky photometric 3 S I M U L AT I O N R E S U LT S A N D D I S C U S S I O N We now discuss the results of our simulated population of CGBs. Primarily, we consider our main simulation (resulting from a projected separation distribution that is flat in log s out to 50 kau), which likely represents an upper bound on the overall population size. However, at the end of this discussion, we present CGB subset sizes for both flat and sloping separation distributions, where the sloping distribution is a closer match to observations (albeit with biases and selection effects), and thus provides a likely lower bound. The output of our main simulation consists of 36 559 ultracool companions (with Teff < 2800 K) with J ≤ 19 mag or W2 ≤ 15.95 mag. When we consider our statistical requirements to confirm companionship, this reduces to 24 196 CGBs. In the following subsections, we discuss the distribution of this CGB sample within intrinsic and observable parameter space, and then consider MNRAS 470, 4885–4907 (2017) Figure 1. The mass–age distribution of our simulated CGBs, with spectral types M7–M9, L and T plotted in grey, red and blue, respectively. Downloaded from https://academic.oup.com/mnras/article/470/4/4885/3869256 by guest on 09 December 2021 We determined a series of companionship probabilities for each simulated benchmark system, appropriate for the observable properties that would be available at each stage of a search-and-follow-up programme. This companion confirmation programme was represented through the following stages: (i) cross-matching the GUMS primary with the simulated field + benchmark UCDs out to the separation of the simulated companion, to account for both crosscontamination (i.e. a UCD companion to star A being erroneously associated with nearby star B), and potential companion mimics whose spectrophotometric distance is consistent with the parallax distance of the primary (within 2σ ); (ii) obtaining the proper motions of the candidate primary and companion and ensuring that they are consistent (within 2σ ); (iii) obtaining the RVs of the candidate primary and companion and ensuring that these too are consistent (within 2σ ). While the 2σ distance criteria are prone to contamination from background unresolved binary UCDs (whose underestimated spectrophotometric distance may fall within the 2σ distance range of the primary), it also makes it unlikely that unresolved binary UCD companions will be mistakenly rejected. This is because such unresolved binaries are overluminous by no more than 0.75 mag (the equalmass limit), which is within two times our adopted σ = 0.4 mag uncertainty. For each benchmark UCD companion, ‘mimics’ were sought in the field population (i.e. UCDs that meet the observational requirements for companionship). This simulate-and-search exercise was carried out 10 000 times following a Monte Carlo approach, for each search-and-follow-up stage. A ‘false alarm probability’ was then determined equal to the number of trials where at least one mimic was found divided by 10 000. And the companionship probability was set as one minus the false alarm probability. Our approach cannot accurately calculate very small (<0.01 per cent but non-zero) false alarm probabilities; however, it is effective at identifying systems with a strong companionship probability. We chose a minimum threshold for companionship probability of 99.99 per cent (close to 4σ confidence) for the confirmation of simulated benchmarks. Some benchmark systems were confirmed after early stages of our search-and-follow-up (see discussion in the next section), but we made our full selection of confirmed simulated systems by applying the threshold at the final stage. We refer to this full sample as ‘confirmable Gaia benchmarks’ (CGBs). prioritized CGB subsets and an optimized search-and-follow-up approach. Gaia benchmarks Figure 3. The projected separation s versus age distribution of our simulated CGBs. UCD spectral types are coloured as in Fig. 1. limits. Our one simulated Y CGB would be within the WISE AllSky Survey, and would be bright enough for spectroscopic follow-up with current facilities (Section 3.3 provides further discussion on the potential for Y dwarf CGBs). These results are summarized in Table 1, and overall suggest a potentially very large CGB population. Although dominated by low-mass stars and late M dwarfs, there should be numerically substantial samples of L and T CGBs across a wide range of age and kinematic population. 3.1.2 Binary constituents Fig. 2 shows the distribution of secondary versus primary Teff for the CGB population. About 60 per cent of CGBs have M dwarf primaries, with primary Teff down to ∼3000 K (∼M5). Most of the remainder have FGK primaries, with just 76 CGBs containing hotter BA-type primaries. In addition to the main-sequence primaries, there are 75 CGBs with sub-giant primaries, and 172 with WD primaries. Below Teff ∼ 3000 K, we observe a sharp drop in the number of primaries. This is because the absolute magnitude sequence for M dwarfs is very steep for optical bands (dropping nearly 3 mag between M5 and M7 in the SDSS r band; Bochanski, Hawley & West 2011), and therefore the G < 20.7 mag limit results in a sharp cut-off in the population. As was discussed by Pinfield et al. (2006), sub-giants and WDs make very useful benchmark primaries. It is possible to constrain the metallicity and ages of subgiant stars quite accurately (as they evolve relatively quickly across the Hertzsprung-Russell diagram) using well-understood models. WD primaries provide lower limit system ages from their cooling age. Furthermore, higher mass WDs have higher mass shorter lived progenitors and the cooling ages will be a better proxy for total system age. 3.1.3 Projected separation Fig. 3 shows the projected separation versus age distribution of the CGBs, with UCD spectral types coloured as in Fig. 1. As we described in Section 2.2.2, our initial separation distribution is flat in Figure 4. The projected separation versus distance distribution of our simulated CGBs. UCD spectral types are coloured as in Fig. 1. The two dotted lines represent the 4 and 20 arcsec limits (see Section 2.3). log s, and is truncated at 50 kau. This truncation is seen in Fig. 3, as is the effect of dynamical break-up that removes CGBs if their age exceeds the (mass and separation sensitive) dynamical-interaction lifetime. This dynamical effect essentially leads to a reduced truncation across the old disc, thick disc and halo, but also thins the CGB population for separations greater than a few thousand au. It is interesting to note that the dynamical-interaction lifetime limits all thick disc CGBs to separations < 20 kau, and all halo CGBs to separations < 10 kau. While our input assumptions about the separation distribution have some inherent uncertainties, our simulation results provide some useful constraints on suitable limits for the separation of CGBs across a range of kinematic populations. 3.1.4 Distance Fig. 4 shows the distance versus projected separation distribution of the CGBs, with UCD spectral types coloured as in Fig. 1. Late M-, L- and T-type CGBs are available out to distances of ∼550, 400 and 70 pc, respectively, and with numbers very limited for distances < 20 pc. There is a fairly uniform increase in the number of Ltype CGBs over the distance range 50–250 pc, since the increase in space volume at larger distance is counteracted by the decrease MNRAS 470, 4885–4907 (2017) Downloaded from https://academic.oup.com/mnras/article/470/4/4885/3869256 by guest on 09 December 2021 Figure 2. The distribution of secondary versus primary Teff for our simulated CGB population. Spectral type divisions are indicated along the top and right axes. Different primary types are plotted in different shades (dark grey, light grey and black for main-sequence stars, sub-giants and WDs, respectively). We overplot density contours over the main-sequence samples to show the structure of the distribution in highly crowded regions. Contour labels are, from the outermost to the innermost, 10, 20, 50 and 100 objects per 25 × 100 K2 bin. 4891 4892 F. Marocco et al. in the range of L sub-types that are detectable at this distance (i.e. all L CGBs can be detected at ∼50 pc whereas only early L CGBs are detectable at ∼400 pc; see also Caballero et al. 2008). Dashed lines delineate regions where CGBs are undetectable in the NIR and MIR surveys because they are unresolved from their primaries (as dictated by our minimum angular separations limits of 4 and 20 arcsec, respectively; see Section 2.2.2). The NIR angular resolution limit has a much greater impact because the majority of CGBs are detectable in the J band (this will be discussed further in Section 3.2). 3.1.5 Metallicity Fig. 5 shows UCD Teff versus metallicity for the CGBs. UCD spectral type divisions and approximate metallicity class ranges (Lépine et al. 2007; Zhang et al. 2017) are indicated along the right and top axes, with sd standing for ‘sub-dwarf’, esd for ‘extreme subdwarf’ and usd for ‘ultra sub-dwarf’. There is a sizeable subset of 735 metal-rich ([Fe/H]> 0.2 dex) M-type CGBs, and a smaller but significant subset of 104 metal-rich L types (though there are very few metal-rich T types). Within metallicity classes, a large subset of 2098 sdM CGBs should be available, with smaller subsets of 99 esdM and 12 usdM CGBs. In addition, there is a subset of 149 sdL CGBs, as well as 4 esdL and 1 usdL types. Our simulation also contains 10 sdT CGBs. For these metal-rich/poor CGBs, about 53 per cent have M dwarf primaries and most of the remainder have FGK primaries (as was discussed in Section 3.1.2, and summarized in Table 1). CGBs with sub-giant or WD primaries are predominantly solar metallicity dwarfs. Most of the CGBs with sub-giant primaries are late M type (save for one L type). CGBs with WD primaries are mostly (∼80 per cent) late M type (cf. Day-Jones et al. 2008), with the remainder generally L type, except for two T-type CGBs (cf. Day-Jones et al. 2011b). MNRAS 470, 4885–4907 (2017) Figure 6. The Gaia G-band magnitudes of CGBs versus V-band primary magnitudes, with UCD spectral types coloured as in Fig. 1. Dashed lines indicate the Hipparcos limit (V ∼ 10 mag) for the primaries, and the Gaia G detection limit for the CGBs. 3.2 Observable properties of CGBs 3.2.1 Gaia combined with infrared surveys Fig. 6 shows the Gaia G-band magnitudes of CGBs versus V-band primary magnitudes, with UCD spectral types coloured as in Fig. 1. Dashed lines indicate the Hipparcos limit (V ∼ 10 mag) for the primaries, and the Gaia G detection limit for the CGBs. To determine how the use of Gaia primaries improves over samples with Hipparcos/Gliese primaries, we counted simulated CGBs in which the primaries have V < 10 mag or distance < 25 pc. This produced 583 Hipparcos/Gliese systems, including 129 L dwarfs and 18 T dwarfs. The entire simulated CGB sample thus represents a 40-fold increase over samples with Hipparcos/Gliese primaries. Also, to compare the approach of using infrared surveys (for CGB detection) to detecting UCDs with Gaia itself, we counted simulated CGBs with Gaia G < 20.7 mag giving 2960 systems. Most of these are late M dwarfs, with 125 L dwarfs and no T dwarfs. The infrared surveys thus improve CGB sample size by 10-fold for late M and 30-fold for L dwarfs compared to Gaia alone. They also provide sensitivity to > 150 T-type CGBs that are undetectable with Gaia. It is also interesting to compare the predictions of our simulation with the results of previous work to identify large samples of wide binaries