The $\textit{Gaia}$ Ultra-Cool Dwarf Sample $-$ III: Seven new multiple systems containing at least one $\textit{Gaia}$ DR2 ultra-cool dwarf

We present ten new ultra-cool dwarfs in seven wide binary systems discovered using $\textit{Gaia}$ DR2 data, identified as part of our $\textit{Gaia}$ Ultra-Cool Dwarf Sample project. The seven systems presented here include an L1 companion to the G5 IV star HD 164507, an L1: companion to the V478 Lyr AB system, an L2 companion to the metal-poor K5 V star CD-28 8692, an M9 V companion to the young variable K0 V star LT UMa, and three low-mass binaries consisting of late Ms and early Ls. The HD 164507, CD-28 8692, V478 Lyr, and LT UMa systems are particularly important benchmarks, because the primaries are well characterised and offer excellent constraints on the atmospheric parameters and ages of the companions. We find that the M8 V star 2MASS J23253550+4608163 is $\sim$2.5 mag overluminous compared to M dwarfs of similar spectral type, but at the same time it does not exhibit obvious peculiarities in its near-infrared spectrum. Its overluminosity cannot be explained by unresolved binarity alone. Finally, we present an L1+L2 system with a projected physical separation of 959 au, making this the widest L+L binary currently known.


INTRODUCTION
Ultra-cool dwarfs (UCDs, spectral type ≥M7) in binary systems with main sequence and post-main sequence stars are valuable benchmarks (Pinfield et al. 2006), providing robust tests of ultra-cool atmospheric and evolutionary models. Under the reasonable assumption of common origin, a bright and temperature with the gaseous giant planets in exosolar systems (e.g. Faherty et al. 2016), but can be studied without the additional complication of the planets' vicinity to a bright host star. A full understanding of ultra-cool atmospheres is therefore of vital importance if we wish to understand exoplanets and their formation and evolution.
The recent second data release from the ESA mission Gaia (Gaia Collaboration et al. 2016aCollaboration et al. , 2018 provides exquisite astrometry for ∼1.3 billion objects within our Galaxy (Lindegren et al. 2018), allowing access to a huge population of wide binaries consisting of an ultra-cool dwarf in a system with a star or white dwarf . In particular, the greatly increased volume probed by Gaia, and the resulting increased pool of potential primary stars, offers for the first time the possibility to map the full age-temperature-metallicity parameter space, large regions of which are currently undersampled or completely unexplored (see e.g. Day-Jones et al. 2011;Deacon et al. 2014;Marocco et al. 2017). While the advent of Gaia expands the pool of potential primaries, existing optical and near-infrared surveys, and the astrometric catalogues that spawned from them (e.g. ULAS, Smith et al. 2014;VIRAC, Smith et al. 2018;CatWISE, Eisenhardt et al. 2019), grant access to a vast population of ultra-cool dwarfs across spectral types M, L, T, and Y. We have therefore set out to complete the nearby census of these objects, to fully explore and characterise ultra-cool atmospheres.
In this paper we present seven new multiple systems containing at least one Gaia DR2 ultra-cool dwarf component.
In Section 2 we describe our candidate selection; in Section 3 we summarise observing and data reduction procedures; in Section 4 we discuss in more detail the newly discovered systems; in Section 5 we compare the main features in the spectra of the new UCDs; finally in Section 6 we summarise our findings and discuss future work.

CANDIDATE SELECTION
We identified an initial list of 8013 UCD candidates from the Gaia DR2 catalogue as follows. First, we queried the catalogue for Gaia sources fainter than the maximum brightness that an UCD at the measured parallax could have, as predicted by the BT-Settl models (Allard et al. 2012a(Allard et al. , 2013. The maximum distance is 373 pc, the distance at which the brightest, hottest UCD (T eff ∼ 2500 K) would be fainter than the Gaia limiting magnitude (G = 20.7 mag). We required the G − G RP colour to be redder than 1.4 mag (since UCDs are typically redder than that; Smart et al. 2017Smart et al. , 2019. To minimise the number of sources with spurious astrometric measurements, we removed candidates within 5 degrees of the Galactic plane and inside an ellipse centred at the Galactic centre with semi-major axis along the Galactic longitude axis of 50 • , and 8 • along the Galactic latitude axis. To retain only reliable astrometric measurements, we required sources to have more than six visibility periods and astrometric excess noise lower than 5 mas. We computed posterior probability densities of the distance given the parallax measurements and associated uncertainties using an exponentially decreasing constant volume density prior, and selected sources with a posterior probability to be within 373 pc greater than 0.5. We then fit a principal curve (Hastie & Stuetzle 1989) in the M G versus G − G RP plane to the values of the resulting set, and calibrated the curve in effective temperature using the spectral types of sources in the Gaia Ultra-Cool Dwarfs Sample (Smart et al. 2017(Smart et al. , 2019 and the Stephens et al. (2009) conversion between spectral types and effective temperatures. Finally, we computed the projections of the UCD candidate positions in the M G versus G RP plane along the principal curve and assigned effective temperatures accordingly. A cut at T eff ≤ 2500 K resulted in the 8013 candidates mentioned above.
We searched for binaries among these 8013 UCD candidates using the criteria defined in Smart et al. (2019, hereafter GUCDS II): where ρ is the separation on the sky in arcseconds, ∆ is the difference between the candidate UCD and primary parallax, and σ are the parallax and parallax uncertainty for the UCD (in mas), ∆µ is the difference of the total proper motions, and ∆θ is the difference of the position angles. The maximum ρ was chosen to correspond to 100,000 au as a conservative upper limit for the projected physical separation (s). This separation meets the binding energy criterion of |U * g | > 10 33 J as developed by Caballero (2009) for a system of a 0.1 and a 2 M objects. The parallax criterion is a compromise between a standard 3σ criterion, and a more conservative 1.0 mas difference to allow for parallaxes that had unrealistically low errors. For the proper motion, using a standard 3σ criterion would remove nearby objects with significant orbital motion, so we choose a conservative 10% agreement, which is large enough to accommodate most orbital motions but small enough to reduce false positives. As discussed in GUCDS II these criteria fail for the nearby binary systems GJ 1048 A/B and G 239-25 A/B (in both cases because the modulus of the proper motions differs by more than 10%). Therefore, our catalogue of binary candidates should not be regarded as complete.
Of the 8013 UCD candidates, 840 have a possible companion according to the criteria above. The seven systems presented here are those that we could observe during our observing nights at the Palomar Observatory. We present their astrometric properties and spectral types in Tables 1 and 2. We collected optical and near-infrared photometry for both components of our newly discovered systems from Gaia DR2, 2MASS (Skrutskie et al. 2006), PanSTARRS DR1 (Chambers et al. 2016), and AllWISE (Cutri et al. 2013). The photometry is also presented in Tables 1 and 2. In Figure 1 we show a colour-magnitude diagram based on Gaia colours and astrometry. The small grey points are objects in Gaia DR2 nominally within 50 pc, selected using Equation C.1 and C.2 from Lindegren et al. (2018). Red points are UCDs identified in Gaia DR2 by GUCDS II. The position of the seven systems presented here is highlighted with different symbols, with the primary plotted in blue and the companion in green. Two objects stand out at first glance: HD 164507 B, and 2MASS J23253550+4608163. We will discuss their properties in Sections 4.1 and 4.4.  Figure 1. Colour-magnitude diagrams depicting the full stellar sequence (top) and a zoom into the ultra-cool dwarfs region (bottom). The small grey points are stars in Gaia DR2 nominally within 50 pc, selected using the criteria described in Appendix C of Lindegren et al. (2018). Red points are the ultra-cool dwarfs identified in Gaia DR2 by GUCDS II. The seven systems presented here are plotted with different symbols, with the primary in each system plotted in blue and the companion in green. Vertical error bars are typically smaller than the symbols. Detailed analysis of individual systems can be found in Section 4.

OBSERVATIONS
We obtained near-infrared spectra for the ultra-cool dwarfs in our newly discovered binary systems using TripleSpec on the 200" telescope at the Palomar Observatory on 2018 April 27-29, 2018 October 16 and 18, and 2019 April 16 (Proposals 2018A J12, 2018B J08, and 2019A J14; PI: Mamajek; see Appendix A). TripleSpec is a near-infrared echelle spectrograph, that delivers a resolution of 2500-2700 over the wavelength range 1.0 − 2.4µm (Herter et al. 2008). Targets were observed following a standard ABBA nodding pattern with a nod throw of 11 arcsec. The slit was aligned to the parallactic angle to minimise atmospheric distortion, with the exceptions of HD 164507 B, V478 Lyr C, and 2MASS J232535.09+460809.3, for which we rotated the slit to avoid the bright primary. We observed an A0 V star (selected using the Gemini Telluric Standard Search on-line tool 1 ) for telluric correction after each target, matching the airmass of observation as closely as possible.
The data were reduced using a modified version of the IDL package Spextool (Cushing et al. 2004). The program applies basic calibration (dark subtraction and flat fielding), then pair-wise subtracts the images to remove sky background. The individual orders of the echelle spectra are traced and extracted, and wavelength calibration is achieved using the numerous OH sky lines. The individual orders are corrected for telluric absorption and flux calibrated using the observed telluric standard star, chosen to match the Vega spectrum used as template in Spextool. The individual orders are then merged, using their overlap to determine flux adjustments when needed. The reduced spectra are presented in Figures 2, 3, and 4.
We assigned a spectral type to our targets via standard template-matching using the classifyByStandard routine in the Python package SPLAT 2 (Burgasser et al. 2016). The code interpolates the templates to the same wavelength grid of the observed spectra, and then minimises the χ 2 of the fit, treating the scaling between the flux-calibrated target and the normalised templates as a free parameter of the fit. The classifyByStandard routine offers the possibility to classify objects by fitting the full spectrum, as well as by fitting only the J band, following the prescriptions of Kirkpatrick et al. (2010). The spectral types obtained with the two methods agree to within ±1 subtype, with the exception of CD-28 8692 B and 2MASS J0139+8110 A. We discuss the discrepancies and our adopted classification in Section 4.3 and 4.5. We used the standard M, L, and T templates defined in Burgasser et al. (2006) and Kirkpatrick et al. (2010), except in the case of V478 Lyr C, where standard templates gave poor fits. Further details on the spectral typing for this source are given in Section 4.2. The results from template matching are presented in Figures 2, 3, and 4, and the assigned spectral types are listed in Tables 1 and 2. 1 https://www.gemini.edu/sciops/instruments/ nearir-resources/spectroscopic-standards-/ telluric-standard-search 2 http://pono.ucsd.edu/~adam/browndwarfs/splat

HD 164507 AB
The primary is a very well characterised G5 IV star that is included in the catalogue of radial velocity standards for Gaia (Soubiran et al. 2013). Several independent estimates of the atmospheric and evolutionary parameters for this subgiant can be found in the literature, and here we briefly summarise those based on high-resolution spectroscopy only. Valenti & Fischer (2005) obtained R ∼70,000 spectroscopy for HD 164507 using the High-Resolution Echelle Spectrometer (HIRES) on the 10 m telescope at Keck Observatory (Vogt et al. 1994). They derived atmospheric parameters using version 2.1 of the software package Spectroscopy Made Easy (SME; Valenti & Piskunov 1996) and the atmospheric models by Kurucz (1992). Mass and age for the star were then derived using the Y 2 isochrones (Demarque et al. 2004). Takeda et al. (2007) and Maldonado et al. (2013) derived independent age and mass using the atmospheric parameters from Valenti & Fischer (2005). Takeda et al. (2007) employed the Yale Rotational Evolution Code (YREC) in its non-rotating mode (Demarque et al. 2008) to generate their set of isochrones, while Maldonado et al. (2013) used the Valenti & Fischer (2005) spectroscopic T eff and metallicity together with Hipparcos data as inputs for PARAM 3 (da Silva et al. 2006) to derive age and mass for HD 164507. Jofré et al. (2015) used high-resolution spectroscopy from SOPHIE on the 1.93 m telescope at the Observatoire de Haute-Provence (Perruchot et al. 2008). The fundamental stellar parameters (T eff , log g, [Fe/H], ξ t ) were computed homogeneously using the FUNDPAR code (Saffe 2011). The chemical abundances of 14 elements (Na, Mg, Al, Si, Ca, Sc, Ti, V, Cr, Mn, Co, Ni, Zn, and Ba) were obtained using the 2009 version of the MOOG 4 code (Sneden 1973). Rotational velocities were derived from the full width at half maximum of isolated Fe lines. Again, mass and age were derived using PARAM. Niedzielski et al. (2016) used the High Resolution Spectrograph (Tull 1998) on the Hobby-Eberly Telescope. The T eff , log g, ξ t , and [Fe/H] were obtained from the measured equivalent width of neutral and ionised iron absorption lines, with the TGVIT code (Takeda et al. 2002(Takeda et al. , 2005. The stellar mass and age were determined using a Bayesian method described in Adamczyk et al. (2016), with theoretical stellar models from Bressan et al. (2012) et al. 1993). Abundances and ξ t were calculated using measured equivalent widths and plane-parallel MARCS model atmospheres (Gustafsson et al. 2008), while T eff and log g were computed from broad-band photometry and the photometric calibration of Casagrande et al. (2010). Finally, Luck (2017) determined mass and age using various sets of isochrones from Bertelli et al. (1994), Demarque et al. (2004), Dotter et al. (2008), and the 2016 version of the BaSTI isochrones (Pietrinferni et al. 2004). . Spectral types are assigned using SPLAT (see Section 3), except for LT UMa A, whose spectral types is taken from Strassmeier et al. (2000). Notes on photometry: a contaminated by bright star halo; b saturated; c contaminated by bright star; d contaminated by diffraction spike.
Finally, Gaia DR2 quotes T eff = 5560 +115 −62 K (see Andrae et al. 2018, for details on how Gaia DR2 atmospheric parameters are derived), and the best-fit template used for radial velocity measurement has T eff = 5500 K, log g=3.5, and [Fe/H]=+0.2 (Sartoretti et al. 2018), all in good agreement with the literature values.
The atmospheric parameters discussed above are listed in Table 3. The values derived are in general agreement with each other, and in particular point towards a slightly super solar metallicity ([Fe/H]=0.03-0.19 dex), and an age for the system in the range 3.0 − 5.9 Gyr.
More accurate age constraints on this star will be provided by TESS (Ricker et al. 2015) via gyrochronology, making this system an exquisite benchmark for UCD models and retrieval codes testing (Line et al. 2015;Burningham et al. 2017Burningham et al. , 2013. The L1 companion, HD 164507 B, is an outlier in the colour-magnitude diagram of Figure 1. With a G − G RP colour of 2.028 ± 0.067 mag, it is among the reddest UCDs in the Gaia sample. Objects with similar G − G RP colour are found in GUCDS II to be either tight binaries or suspect tight binaries. The red G − G RP colour in this case would be due to the fact that G RP (and G BP ) magnitudes are determined by integrating the G RP fluxes in a 3.5 × 2.1 arcsec 2 window, and there is currently no treatment of multiple sources in the same window in Gaia DR2 (Evans et al. 2018). Therefore, an excess in G RP for close binary systems is expected. However, there is no evidence for binarity of HD 164507 B. The source is not resolved by Gaia, and the goodness-of-fit and astrometric excess noise reported in Gaia DR2 (2.7359 and 2.108 mas, respectively) are both con-sistent with the mean values for UCDs found in GUCDS II (5.2 ± 2.6 and 2.2 ± 1.2 mas, respectively). The primary has higher-than-solar metallicity (0.03 < [Fe/H] < 0.19 dex, see Section 4.1), and higher metallicity UCDs are expected to have redder than average colours because of the enhanced dust content in their photosphere (e.g. Looper et al. 2008;Marocco et al. 2014). However the near-infrared spectrum of HD 164507 B does not show obvious peculiarities (see Figure 2). Finally, youth is also typically associated with redder-than-usual colours (see e.g. Faherty et al. 2016), but young and suspected young objects in GUCDS II form a relatively tight sequence with 1.6 G − G RP 1.8 mag, and the age of the system rules out youth as a cause. Optical spectroscopy for this UCD is desirable to shed light on its nature.
We derive T eff for the companion using the Filippazzo et al. (2015) spectral type to T eff polynomial relation, and obtain T eff = 2100 ± 29 K. Linear interpolation of the BT-Settl isochrones 5 for solar and super-solar metallicity ([Fe/H] = +0.5 dex) in the age range 3.0-5.9 Gyr, and for T eff = 2100 ± 29 K, gives a mass for the companion in the range 50-77 M Jup , at or below the hydrogen burning limit.
To compute the bolometric luminosity (L bol ) we need to determine a bolometric correction, since our TripleSpec spectrum only covers the 1.0 < λ < 2.4 µm range. We did this by fitting the TripleSpec spectrum with the BT-Settl atmospheric models (Allard et al. 2012b) with the fitting technique developed by Cushing et al. (2008). The models cover the T eff space in steps of 50 K, the log g space in steps of 0.5 dex, and the [Fe/H] space in steps of 0.5 dex.
We flux calibrated the target's spectrum using the measured 2MASS J-band magnitude, and then allowed the scaling factor between the flux-calibrated spectrum and the models to be a parameter of the fit. The best fit scaling factor gave us a measurement of the radius (R) of the target via the simple geometric dilution factor (R/d) 2 . We restricted the range of models to be considered for fitting to the ± 200 K range around the predicted T eff of 2100 K and the metallicity to be within ±0.5 dex of the metallicity of the primary, for which we chose the mid point of the values quoted in the literature, i.e. 0.11 dex.
We used the scaled best-fit atmospheric model to complete the TripleSpec spectrum at long and short wavelength (λ < 1µm and λ > 2.4µm). L bol was then computed by summing the flux density over the full model+TripleSpec spectrum, and multiplying it by 4π d 2 . The uncertainty on L bol was computed by propagating the uncertainty on the measured spectrum, as well as the uncertainty on the 2MASS magnitude used for flux calibration, and the uncertainty on the distance. The best fit model for HD 164507 B has T eff = 2300 K, log g = 5.0, and [Fe/H] = +0.5 dex. The radius corresponding to the best fit scale factor is 0.88 R Jup , and the bolometric luminosity is log 10 (L bol /L ) = −3.144 +0.039 −0.043 . Approximately 17% of the bolometric luminosity reported here is outside of the TripleSpec wavelength range (1.0-2.4 µm). This fraction decreases with spectral type, as the contribution from the optical portion of the spectral energy distribution collapses, while the longer wavelength flux does not increase significantly. The model-dependent fraction of L bol approaches ∼ 40% for the late-Ms in our sample, and decreases down to ∼ 8% for the L2s. The best fit model for HD 164507 B is shown in Figure 5. The overall fit is poor: (i) the model has a triangular H band, while the target has a much flatter H-band spectrum; (ii) the alkali lines in the J band are much too shallow in the model compared to the observed ones; (iii) the K-band spectrum in the model is too flat, and (iv) the overall spectrum is too blue compared to our target. The best fit T eff is 200 K warmer than the prediction from the Filippazzo et al. (2015) polynomial.

V478 Lyr ABC
The primary is a chromospherically active G8 V single-lined spectroscopic binary with a period of about 2.13 d (Fekel 1988). This star was found to have strong ultraviolet emission features and a filled-in Hα absorption line that is variable in strength. Therefore, Fekel (1988) classified it as an early-type BY Draconis system. The secondary had its mass estimated to be about 0.3 M and to be probably an M2-M3 dwarf. The inclination of the system was measured to be 67 ± 12 deg. The lithium abundance of the G8 dwarf, estimated from the equivalent width of the Li i 6707.8Å line (47 mÅ), led Fekel (1988) to propose an age for the system that is somewhat less than that of the Hyades cluster (680 Myr; Gossage et al. 2018).
Using the BANYAN Σ on-line tool, the Gaia DR2 astrometry, and the mean radial velocity from Nordström et al. (2004), we find a probability of 0% for the object to be a member of any of the young moving groups considered in BANYAN Σ (including the Hyades).
Nevertheless, the UCD companion, dubbed V478 Lyr C, shows a somewhat triangular H-band spectrum, a feature previously associated with youth (Lucas et al. 2001;Allers & Liu 2013). Gravity-sensitive spectral indices and pseudoequivalent width defined in Allers & Liu (2013) however lead to a L1 field surface gravity (FLD-G) classification for the companion. Intermediate surface gravity (INT-G) and verylow surface gravity (VL-G) objects in the Allers & Liu (2013) sample have typical age <200 Myr and, according to a more recent study conducted by Martin et al. (2017), the reliability of the gravity classification drops significantly for objects with age >100 Myr. On the other hand, the L1 companion to the young A3V star β Circini has a flat H-band spectrum (and no low gravity features, see Smith et al. 2015). The age of the β Circini system has been estimated to be in the 370-500 Myr range. We would therefore expect the V478 Lyr system to be somewhat younger than the β Circini system, but likely older than ∼100 Myr.
Spectral typing via standard template matching leads to an L2 type. However the fit in Figure 4 is poor, with the standard not only failing to match the H band shape, but also underestimating the flux at the blue end of the spectrum (up to ∼ 1.2µm). Using the Kirkpatrick et al. (2010) method, i.e. fitting only the 0.9-1.4 µm range, the best fit template is the L1 standard, but the target shows flux excess at the longer wavelength, as expected for a low surface gravity object.
If we fit V478 Lyr C with the low gravity templates defined in Allers & Liu (2013), the best fit is the L0β standard. The fit to the H band is much more accurate, and the flux in the J band is less underestimated, but at the same time the fit to the H 2 O band at ∼ 1.4 µm is poorer. Given all of the above, we assign V478 Lyr C a spectral type of L1:. Filippazzo et al. (2015) derived a M H to T eff polynomial relation for young objects, but the available nearinfrared photometry for V478 Lyr C is heavily contaminated by the parent star (at ρ ∼ 17 arcsec). We computed a synthetic H-band magnitude using our flux-calibrated Triple-Spec spectrum and the 2MASS H-band response curve (Cohen et al. 2003). We estimated the accuracy of our synthetic H magnitude by comparing the synthetic magnitudes obtained for the other objects observed as part of our Triple-Spec run, against their measured 2MASS H (for all except HD 164507 B, since its photometry is also contaminated). The mean offset between our synthetic magnitudes and the measured ones is -0.007 mag and the 1σ dispersion around the mean is 0.44 mag. We therefore adopted 13.74±0.44 mag as our synthetic measurement, and obtain T eff = 1740 ± 130 K for V478 Lyr C as a result. Linear interpolation of the BT-Settl isochrones for solar metallicity in the age range 0.10 − 0.37 Gyr gives a mass for this object in the range 10-28 M Jup , straddling the deuterium fusion mass limit.
We determined L bol for V478 Lyr C following the same procedure described in Section 4.1. The best fit model has T eff = 1800 K, log g = 5.0, and solar metallicity. The log g = 5.0 is somewhat higher than one might expect, given the age of the system, and the fact that this object shows signs of youth. The radius resulting from the best-fit scaling factor is 1.31 R Jup and the bolometric luminosity is  (2005)   log 10 (L bol /L ) = −3.33 +0.26 −0.78 . The best fit model can be seen in Figure 5. The overall fit is good, with the model only slightly underpredicting the flux at the shortest wavelength (λ < 1.2 µm) but that is the region of lowest signal-to-noiseratio. Oh et al. (2017) found the SB1 primary to form a very wide co-moving pair with the G6V HD 171067, with a projected separation of ∼8 pc. The Oh et al. (2017) analysis however did not take into account radial velocity (RV). The measured system RV for V478 Lyr AB is −25.2 ± 4.8 km s −1 (Nordström et al. 2004), and is discrepant from the RV of HD 171067 (−46.197 ± 0.002 km s −1 ; Soubiran et al. 2013). As a result, the G6V is unlikely to be associated with the V478 Lyr triple system. V478 Lyr ABC joins the rank of triple systems consisting of a spectroscopic binary with a wide, low-mass tertiary component (see Allen et al. 2012, and references therein). These systems are precious for testing formation simulations of very close separation binaries, which require a mechanism to draw angular momentum away from an already close pair of objects. One proposed mechanism is through three-body interactions with cool dwarfs (see e.g. Sterzik & Durisen 2003;Delgado-Donate et al. 2004;Umbreit et al. 2005), and a key observable to test such scenario is the fraction of tight spectroscopic binaries that have a wide additional companion. Towards this goal, V478 Lyr AB was among the stars targeted by Allen et al. (2012), who conducted a deep nearinfrared survey looking for low-mass tertiary components around 118 known spectroscopic binaries within 30 pc of the Sun. However, V478 Lyr C was missed probably because of the combination of its tight angular separation from the binary (17.05 arcsec, close to the Allen et al. 2012 survey limit of 10-15 arcsec), the large magnitude difference between SB1 primary and L dwarf companion, and the large contamination by reddened background sources resulting from its proximity to the Galactic plane (b = 10.1 deg).
Finally, the estimated orbital period for this system is 8,000 yr, despite this being the most favorable configuration among the seven systems presented here -i.e. a relatively massive primary, with a relatively tight separation, and assuming a face-on circular orbit. If instead we assume the wide L1: companion is coplanar with the SB1, i.e. that the inclination angle is 67 ± 12 deg, then the orbital period would be ∼9,700 yr. In either case, no dynamical mass measurement is possible for the UCD. The other systems presented here have even longer estimated orbital periods.

CD−28 8692 AB
The primary is a slightly metal poor K5 V star. It has been monitored with HARPS for planets by Sousa et al. (2011), who found no evidence for RV variations. Sousa et al. (2011) also used the HARPS spectra to determine atmospheric parameters, and obtained T eff = 4799±90 K, log g = 4.43±0.18, and [Fe/H] = −0.22 ± 0.06 dex. They then estimated a mass of 0.715 ± 0.014M for the star using the measured atmospheric parameters and the Padova isochrones. Adibekyan et al. (2012) used the atmospheric parameters estimated by Sousa et al. (2011) and the HARPS spectra to measure detailed abundances of 12 chemical species, with typical precision in the 0.035 − 0.260 dex range.
Later, Delgado Mena et al. (2015) used the HARPS data to estimate atmospheric parameters and combined them with the Li i abundance to infer an age of 4.48 Gyr for this star.
The Gaia DR2 effective temperature for this star is 4742 +138 −116 K (Andrae et al. 2018), while the best-fit template used for radial velocity measurement has T eff = 4750 K, log g=4.5, [Fe/H]=-0.2 dex (Sartoretti et al. 2018). All Gaia DR2 values are in good agreement with the literature measurements.
The companion presented here is classified as L2, with a projected separation of 2026 au (50.91 arcsec). The L2 template is a good fit to the spectrum of the target, with the exception of a slightly suppressed K band (typical of metalpoor and high surface gravity dwarfs; Burgasser et al. 2008;Zhang et al. 2017), and a flux excess at ∼ 1.3 µm. Scatter in the strength of the ∼ 1.3 µm peak among objects of a given spectral type has been observed before (Cruz et al. 2018). The Kirkpatrick et al. (2010) method yielded a very different classification of L6. While the L6 template does indeed provide a slightly better fit to the J band reducing the overluminosity at ∼ 1.3µm, the target is much bluer than the L6 standard at longer wavelength. Low metallicity L dwarfs are indeed slightly bluer compared to their solar metallicity counterparts, but this system is only slightly metal poor, and therefore a large suppression of the H -and K -band flux is unlikely. Moreover, the absolute G magnitude for CD-28 8692 B is 17.406±0.004 mag, which is consistent with the median value for L2s (17.24±0.41 mag; Smart et al. 2019), but nearly two magnitudes overluminous compared to typical L6s (19.25±0.60 mag; Smart et al. 2019). Therefore, we retain a classification of L2 for this object.
Somewhat counterintuitively, the spectral indices for CD-28 8692 B are consistent with an INT-G classification. This is unexpected, since a relatively old, metal-poor object should exhibit surface gravity typical of standard field L dwarfs, or at most slightly higher. The transition between INT-G and FLD-G however is not very sharp, and scatter around the dividing line has been previously noted (Martin et al. 2017). The unusual metallicity of the CD-28 8692 AB system further affects the reliability of the gravity classification, as first noticed by Aganze et al. (2016) for the M9.5 companion to the metal poor M1 V GJ 660.1A ([Fe/H] = -0.63±0.06 dex). We therefore conclude that our INT-G classification for CD-28 8692 B is incorrect. The solar metallicity BT-Settl isochrones at T eff = 1960± 29 K (as given by the Filippazzo et al. 2015 polynomial relations) and age = 4.48 Gyr gives a mass of ∼70 M Jup . Although the system is slightly metal poor, we cannot use the publicly available BT-Settl isochrones for low metallicity ([Fe/H] = -0.5 dex), since they do not extend below 75 M Jup and T eff ∼ 3000 K.
We determined L bol for CD-28 8692 B following the same procedure described in Section 4.1. The best fit model has T eff = 1800 K, log g = 5.0, and solar metallicity. We determine a radius of 0.87 R Jup , and log 10 (L bol /L ) = −3.688 +0.047 −0.053 . The best fit model is shown in Figure 5. The model fit is of good quality, the main discrepancies being in the blue wing of the H band (the model underpredicting the observed flux) and at ∼1.3 µm, where the model does not correctly reproduce the sharp observed peak (see above).

2MASS J23253550+4608163 + 2MASS J23253519+4608098
2MASS J23253550+4608163 is overluminous compared to objects of similar G − G RP colour and spectral type. Typical M8 dwarfs have M G = 15.24 ± 0.63 mag (see GUCDS II), while our target has M G = 12.850 ± 0.004 mag 6 . The overluminosity cannot be explained by unresolved binarity alone, since an equal-mass binary would at most be 0.75 mag overluminous, while the target is almost 2.4 mag overluminous. Young objects can also be redder and overluminous compared to field-age objects. However, 2MASS J23253550+4608163 does not show any indication of youth in its near-infrared spectrum (see Figure 3, middle panel) and its kinematics are inconsistent with membership to any of the young moving groups using the BANYAN Σ on-line tool 7 (Gagné et al. 2018). Contamination by a background object could be another possibility, and this source is indeed flagged as duplicate (duplicated_source = 1), however the background object would need to have the G − G RP colour of a late-M dwarf, since the G − G RP of 2MASS J23253550+4608163 is in line with the median colour of M8 dwarfs (1.592 ± 0.005 mag vs. 1.61 ± 0.95 mag, see GUCDS II). External photometry from 2MASS, PanSTARRS-1 and AllWISE does not show evidence of contamination nor peculiar colours, but all absolute magnitudes are similarly overluminous when compared with M8 dwarfs.
An indication of possible problems is the relatively large goodness-of-fit (astrometric_gof_al) of 132 (cf. the mean value of 5.2 ± 2.6 for objects in GUCDS II), which may indicate that the parallax for this source is spuriously large. The companion, 2MASS J23253519+4608098, does not show any sign of peculiarity, neither photometric nor spectroscopic. This could therefore be an unfortunate case of chance alignment, with 2MASS J23253550+4608163 being a background M dwarf whose spurious astrometry is consistent, by chance, with being a companion to 2MASS J23253519+460898. The astrometry for 2MASS J23253519+46089 would instead be correct. The chance of such an unfortunate alignment is however extremely low, given the tight separation of the pair on the sky (7.24"). We therefore have no conclusive explanation for the overluminosity of this object.
We determined L bol for both components of this system following the same procedure described in Section 4.1. The best fit model for the A component has T eff = 2400 K, log g = 5.0, and [Fe/H] = +0.5 dex. The radius is 3.14 R Jup (c.f. model-predicted value of 2.33 R Jup ), which is unusually large for an ultra-cool dwarf, but probably a consequence of the overluminosity discussed above. The result is log 10 (L bol /L ) = −1.928 +0.026 −0.027 . The best fit model for the B component has T eff = 1800 K, log g = 5.5, and solar metallicity. The radius corresponding to the best-fit scaling factor is R = 1.39 R Jup which is somewhat large for an object with this temperature and surface gravity (the BT-Settl models predict R ∼ 0.9 R Jup ). The bolometric luminosity is log 10 (L bol /L ) = −3.265 +0.053 −0.060 . The best fit models for both components are shown in Figure 6. The fit to the spectrum of 2MASS J23253550+4608163 is overall poor. The model appears too blue compared to the observed spectrum with the flux at λ < 1.3µm being overestimated and the flux in the K band being underestimated. The shape of the H band is also poorly reproduced, with the model having a more pronounced peak, while the observed spectrum appears flatter. The fit to the L dwarf component, 2MASS J23253519+4608098, is good, with the model only slightly underpredicting the flux at λ < 1.25µm.

2MASS J01390902+8110003 + 2MASS J01385969+8110084
With a projected separation of 959 au, this system is to our knowledge the widest L+L dwarf binary known to date. The primary is an L1 based on the template fitting to the whole spectrum, while a fit to the J band alone results in a significantly earlier spectral type, M8. The discrepancy is mostly driven by the slightly overluminous blue end of the TripleSpec spectrum (λ < 1.1µm, see Figure 2). The L1 standard gives a good fit to the overall spectrum except for this wavelength range, which is however also the lowest signal-tonoise-ratio portion of the spectrum. On the other hand, the M8 template reproduces better this part of the spectrum, but starts to diverge from the observations at wavelength longer than ∼ 1.3µm, with the target being overall redder than the template. While in principle this could be evidence of youth, the morphology of the H band, and the depth of the Na i and K i absorption lines suggest that the object is not particularly young. We assume a spectral type of L1 for this object in the rest of the analysis. The companion is classified as L2 by both methods.
Various authors have focused on the identification of wide low-mass binaries. Recent examples include SLoW-PoKES (Sloan Low-mass Wide Pairs Of Kinematically Equivalent Stars; Dhital et al. 2010), Baron et al. (2015), and Gálvez-Ortiz et al. (2017). Extremely wide low mass binaries do exist, with separations out to tens of thousands of au, and are found in young clusters and moving groups (see e.g. GUCDS II, Alonso-Floriano et al. 2015) as well as in the field (Caballero 2012;Caballero & Montes 2012;Dhital et al. 2010). These systems are rare, with an estimated fraction of wide low-mass binaries in the field of 1-2% (Burgasser et al. 2009). Their paucity may be explained via Galactic dynamical evolution, with subsequent stellar encounters in the Galactic disk progressively increasing the separation between the low-mass binary components, eventually leading to its dissolution (Weinberg et al. 1987). This sets a hard lower limit on the binding energy (see e.g. Burgasser et al. 2003;Caballero 2009).
However rare, these systems pose a challenge to the formation models of low-mass stars and brown dwarfs. In particular, Kouwenhoven et al. (2010Kouwenhoven et al. ( , 2011 argued that systems with separation > 1000 au are unlikely to have been formed as primordial binaries (since their orbital separation would be comparable to the size of an embedded cluster), but instead originated during the cluster dissolution process. Dhital et al. (2010) observed a bimodal binary separation (also observed by Kouwenhoven et al. 2010), suggesting the presence of two populations, one old and tightly bound, and another young and weakly bound, recently formed and unlikely to survive more than a few Gyr.
For us to determine how strongly bound this system is, we need to constrain the mass of the components. The spectra, presented in Figure 2, do not present any obvious peculiarity, and both give a good fit to the standard templates. We can therefore reasonably assume that these two L dwarfs are of solar metallicity, and with age > 0.37 Gyr (following the same reasoning used in Section 4.2). We estimate the effective temperature for the two components using the Filippazzo et al. (2015) polynomial relation, and obtain 2100±29 K and 1960±29 K for the L1 and L2, respectively. Given these temperatures, interpolation of the BT-Settl isochrones in the 0.37 − 13 Gyr range gives a mass of 44 − 82 M Jup and 42 − 80 M Jup , respectively, corresponding to a total system mass in the 0.08 − 0.15 M range. The corresponding binding energy for the pair is 3 × 10 33 < |U * g | < 1 × 10 34 J, just above the |U * g | > 10 33 J limit proposed by Caballero (2009).
We can finally estimate how long the 2MASS J01390902+8110003 + 2MASS J01385969+8110084 system is likely to survive stellar encounters in the Galactic disk, using the method described in Dhital et al. (2010). Rearranging their equation 18, and assuming a lower limit on the total mass for this system of 0.08 M , we find that the expected lifetime would be >22 Gyr. Alternatively, we can compute the maximum separation for a binary of given total mass to remain bound for at least 10 Gyr, rearranging equation 28 from Weinberg et al. (1987) and following their assumption of an average Galactic stellar density of 0.16 pc −3 , an average stellar mass of 0.7 M , and a relative velocity for the stellar encounters of ∼20 km s −1 . We find the maximum separation for a system of total mass > 0.08 M to be > 1.5 × 10 3 au. The system is therefore bound.
We determined L bol for both components of the system following the same procedure described in Section 4.1. The best fit model for component A has T eff = 1950 K, log g = 5.5, and solar metallicity. We determine a radius of 2.16 R Jup and log 10 (L bol /L ) = −2.753 +0.029 −0.031 . The radius is unusually large, and inconsistent with the model-predicted radius for an object of such atmospheric properties (1.14 R Jup ).
The best fit model for the B component has T eff = 1800 K, log g = 5.0, and solar metallicity. The resulting radius is 1.58 R Jup and log 10 (L bol /L ) = −3.158 +0.034 −0.037 . The radius is once again unusually large, and even more inconsistent with the model-predicted radius for an object of such atmospheric properties (0.91 R Jup ).The best fit models for both components can be seen in Figure 6. The fit to the L1 (2MASS J01390902+8110003) is overall good, while the fit to the L2 (2MASS J01385969+8110084) is of slightly lower quality. Main discrepacies are an overall underestimated flux in the blue wing of the H band, as well as in the K band and at ∼ 1.3µm.

2MASS J18392917+4424386 + 2MASS J18392740+4424510
This is a very wide (811 au) M+L binary, akin to the 2MASS J01390902+8110003 + 2MASS J01385969+8110084 system. The primary is the only previously known UCD discussed in this paper, and was classified M9 V using IRTF/SpeX spectroscopy in Bardalez Gagliuffi et al. (2014). The TripleSpec spectrum for the companion is presented in Figure 2, and does not present any obvious peculiarity. We classify it as L2 via template matching.
Following the same method described above, we estimate the T eff for the two components to be 2400 ± 29 K and 1960±29 K respectively, leading to masses of 49−88 M Jup and 42 − 80 M Jup . The binding energy of the system is therefore 4 × 10 33 < |U * g | < 1 × 10 34 J. The expected lifetime (computed using the same procedure described in Section 4.5) is > 28 Gyr and the separation limit > 1.6 × 10 3 au. The system is therefore bound.
We determined L bol for the L dwarf following the same procedure described in Section 4.1. The best fit model has T eff = 1800 K, log g = 5.5, and solar metallicity. We determine a radius of 1.19 R Jup and log 10 (L bol /L ) = −3.440 +0.042 −0.046 . The best fit model is presented in Figure 5.

LT UMa AB
LT UMa is a variable star of BY Dra type, with an amplitude of 0.03 mag (no period listed) in The International Variable Star Index 8 , based on 11 observations by Strassmeier et al. (2000).
The companion was first identified by Pinfield et al. (2006) based on motion and colour, but no spectroscopy was presented there. The Washington Double Star Catalog lists the pair as WDS J08448+5532. The spectral types are T eff for HD 164507 B and CD-28 8692 B are computed using the spectral type to T eff polynomial relations for field-age objects derived in Filippazzo et al. (2015), while for V478 Lyr C we used the M H to T eff polynomial relation for young objects presented in the same paper.
48 − 77 M Jup range. We determined L bol for the M dwarf following the same procedure described in Section 4.1. The best fit model has T eff = 2300 K, log g = 5.0, and [Fe/H] = +0.5 dex. The radius is 1.13 R Jup , in good agreement with the model-predicted radius (1.18 R Jup ). The bolometric luminosity is log 10 (L bol /L ) = −2.968 +0.028 −0.029 . The best fit model can be seen in Figure 5. The quality of the fit is poor. The model has a triangular-shaped H band that is not present in the target, which instead displays a flat H-band spectrum. The alkali lines in the J band are also weaker in the model compared to the observed ones.

COMPARISON OF L DWARF SPECTRAL FEATURES
Despite the relatively small sample size, it is nonetheless interesting to compare the spectroscopic features in our newly discovered L companions. In particular, V478 Lyr C, HD 164507 B, and CD-28 8692 B offer an interesting comparison set. With very similar spectral type (L1:, L1, and L2, respectively), but different ages and metallicity, these three objects can be used to qualitatively determine the dependence of spectral features on these parameters. Properties for these three UCDs relevant to this analysis are summarised in Table 5. Figure 8 shows the normalised IR spectra, centred around four of the main absorption features in the spectra of early L dwarfs: the Na i doublet at ∼ 1.139 µm, the K i doublets at ∼ 1.173 µm and ∼ 1.248 µm, and the CO band head at 2.30 µm.
The alkali lines in V478 Lyr C and HD 164507 B show remarkable similarity, while those in CD-28 8692 B are deeper and broader. FeH absorption in the 1.24-1.25µm range is also stronger in CD-28 8692 B, as expected from its age, which confirms the known trend of alkali lines and hydride bands with metallicity (see e.g. Kirkpatrick et al. 2010). Surprisingly, the CO band at 2.293 µm appears deeper in CD-28 8692 B as well, while the band at 2.322 µm is in a region of too low signal-to-noise ratio. While the strength of this CO band is relatively insensitive to changes in effective temperature in the L0-4 range (Cruz et al. 2018), blue L dwarfs and L subdwarfs have weaker CO bands than their solar-metallicity counterparts (see e.g. Zhang et al. 2017). A strong CO band has been previously observed in the blue L1 dwarf 2MASS J17561080+2815238 .
Comparison of the spectral indices and equivalent widths presented in Table 4 as a function of the age and metallicity for these three systems leads to some preliminary considerations: • the "water-based" indices H • the H 2 OD and H 2 O-J indices seem sensitive to age (i.e. surface gravity) too; • the K i lines are sensitive to age (i.e. surface gravity) but also metallicity, becoming stronger (i.e. having larger equivalent width) as age increases, but weaker at higher metallicity. As a result, the young (≈100-370 Myr) L1: V478 Lyr C has K i lines of roughly equal strength as the older (3.0-5.9 Gyr) but metal rich L1 HD 164507 B (5.52, 8.20, 5.05, and 4.89Å for V478 Lyr C vs. 5.97, 6.87, 5.28, and 5.03Å for HD 164507 B).
Followup of a larger sample of benchmark L dwarfs is fundamental to better identify/quantify possible dependencies of the above spectral features on age and metallicity.

CONCLUSIONS
We have presented seven multiple systems discovered in Gaia DR2 data, identified as part of our GUCDS project. The systems presented here include an L1 companion to the G5 IV star HD 164507, an L1: companion to the RS CVn star V478 Lyr, three low-mass binaries consisting of late Ms and early Ls, an L2 companion to the metal-poor K5 V star CD-28 8692, and an M9 V companion to the young variable K0 V star LT UMa. The HD 164507 and CD-28 8692 systems are particularly important benchmarks, because the primaries are very well characterised and offer excellent constraints on the atmospheric parameters of the companion. While the HD 164507 AB system is slightly metal rich, the CD-28 8692 AB system is slightly metal poor, and therefore cover an exotic region of the parameter space, where observational constraints on theoretical models is currently scarce. The V478 Lyr ABC system is a nice addition to the sample of wide low-mass tertiary components to tight binaries, a population of crucial importance to validate formation theories for tight binaries.
We have also reported the discovery of the currently widest L+L binary known -the 2MASS J01390902+8110003 + 2MASS J01385969+8110084 system, with a projected separation of about 960 au. This system, together with the other two wide low-mass wide binaries presented here, pose an increasing challenge to models of formation and evolution of wide low-mass binaries.
A first, qualitative analysis of the sample reveals tentative correlations between spectral indices, equivalent widths, and age and metallicity for the ultra-cool dwarfs presented here. Analysis of a larger sample of benchmarks will provide stronger constraints on such correlations, and Gaia DR2 will play a cornerstone role in shaping our understanding of ultra-cool atmospheres. FM acknowledges support by the NASA Postdoctoral  Program            Features likely due to telluric absorption are labelled with the symbol ⊕. All spectra are smoothed down to a resolution of 3Å pix −1 to reduce the noise. The alkali lines in V478 Lyr C and HD 164507 B show remarkable similarity, while those in CD-28 8692 B are deeper and broader, confirming the known trend with metallicity (see e.g. Kirkpatrick et al. 2010). The CO band head at 2.293 µm appears deeper in CD-28 8692 B as well, while the CO band head at 2.322 µm is in a region of too low signal-to-noise-ratio.