A high binary fraction for the most massive close-in giant planets and brown dwarf desert members

TL;DR

This study reveals a high binary fraction among stars hosting massive close-in giant planets and brown dwarfs, indicating stellar companions significantly influence their formation and evolution, especially for higher-mass and tighter-orbit systems.

Contribution

It provides the first robust statistical evidence linking stellar multiplicity to the presence of massive close-in planetary and substellar companions.

Findings

Binary fraction of ~79% for these systems, higher than field stars.

Higher binary rate observed for systems with shorter orbital periods.

Peak in binary separation at around 250 AU, suggesting influence from wide binaries.

Abstract

Stellar multiplicity is believed to influence planetary formation and evolution, although the precise nature and extent of this role remain ambiguous. We present a study aimed at testing the role of stellar multiplicity in the formation and/or evolution of the most massive, close-in planetary and substellar companions. Using direct imaging observations, as well as the Gaia DR2 catalogue, we searched for wide binary companions to 38 stars hosting massive giant planets or brown dwarfs (M > 7 MJup) on orbits shorter than ~1 AU. We report the discovery of a new component in the WASP-14 system, and present an independent confirmation of a comoving companion to WASP-18. From a robust Bayesian statistical analysis, we derived a binary fraction of 79.0+13.2-14.7% between 20-10,000 AU for our sample, twice as high as for field stars with a 3-{\sigma} significance. This binary frequency was found to be larger than for lower-mass planets on similar orbits, and we observed a marginally higher binary rate for inner companions with periods shorter than 10 days. These results demonstrate that stellar companions greatly influence the formation and/or evolution of these systems, suggesting that the role played by binary companions becomes more important for higher-mass planets, and that this trend may be enhanced for systems with tighter orbits. Our analysis also revealed a peak in binary separation at 250 AU, highlighting a shortfall of close binaries among our sample. This indicates that the mechanisms affecting planet and brown dwarf formation or evolution in binaries must operate from wide separations, although we found that the Kozai-Lidov mechanism is unlikely to be the dominant underlying process. We conclude that binarity plays a crucial role in the existence of very massive short-period giant planets and brown dwarf desert inhabitants, which are almost exclusively observed in multiple systems.

Tables

Table 1. Table 1: Summary of our new observations.

Target	Observation Date	Telescope / Instrument	Filter	Field of View	Pixel Scale	Previous Observations
WASP-18	September 4, 2017	VLT / NACO	L’	$28″\times 28″$	0$\stackrel{\prime \prime }{.}$027	Ngo et al. (2015)
HD 162020	September 6, 2017	VLT / NACO	L’	$28″\times 28″$	0$\stackrel{\prime \prime }{.}$027	Eggenberger et al. (2007)
BD+24 4697	September 6, 2017	Gemini North / NIRI	Ks	$22″\times 22″$	0$\stackrel{\prime \prime }{.}$022	$-$
HD 77065	December 12, 2017	Gemini North / NIRI	Ks	$22″\times 22″$	0$\stackrel{\prime \prime }{.}$022	$-$
HD 134113	June 22, 2018	WIYN / NESSI	562nm, 832nm	$4.6″\times 4.6″$	0$\stackrel{\prime \prime }{.}$040	$-$
HD 160508	June 24, 2018	WIYN / NESSI	562nm, 832nm	$4.6″\times 4.6″$	0$\stackrel{\prime \prime }{.}$040	$-$

Table 2. Table 2: Relative differences in parallax and proper motion, with their associated errors, between the components of all Gaia binaries. Fractional differences are calculated relative to the first component listed for each system (single or binary), and our science targets always correspond to the first component given. The two new binaries identified in this work are marked in bold.

System	$\Delta \pi /{\pi }_0$	$\Delta {\mu }_{\alpha \ast }/{\mu }_{\alpha \ast ,0}$	$\Delta {\mu }_{\delta }/{\mu }_{\delta ,0}$
	(%)	(%)	(%)
11 ComA-B	10.74 $\pm$ 1.92	1.05 $\pm$ 0.30	1.67 $\pm$ 0.33
30 AriBC-A	1.03 $\pm$ 0.39	3.21 $\pm$ 0.11	42.23 $\pm$ 1.81
$\upsilon$ AndA-B	0.48 $\pm$ 0.38	0.10 $\pm$ 0.19	0.26 $\pm$ 0.14
AS 205A-BC	18.41 $\pm$ 1.55	27.25 $\pm$ 7.14	13.83 $\pm$ 1.66
HD 89744A-B	0.43 $\pm$ 2.24	0.53 $\pm$ 0.58	1.40 $\pm$ 0.75
HD 156846A-B	2.44 $\pm$ 1.70	5.46 $\pm$ 0.51	7.23 $\pm$ 0.41
HD 178911B-AC	17.02 $\pm$ 1.58	33.40 $\pm$ 1.20	5.73 $\pm$ 0.34
Kepler-13A-BC	8.61 $\pm$ 5.46	11.39 $\pm$ 6.55	4.85 $\pm$ 2.38
NLTT 41135B-A	0.55 $\pm$ 0.51	5.44 $\pm$ 0.18	0.26 $\pm$ 0.11
WASP-14AB-C	0.98 $\pm$ 1.17	4.34 $\pm$ 0.71	11.51 $\pm$ 2.69
WASP-18A-B	16.85 $\pm$ 14.60	6.30 $\pm$ 5.94	10.78 $\pm$ 9.03

Table 3. Table 3: Detection limits used.

Object ID	Facilities / Instrument	Filter	Phot. (mag)	Limits units	Curve flag	Data Set
30 Ari B	Keck / NIRC2	$J$	$6.080\pm 0.02$	$\Delta$mag	0	Kane et al. (2015)
70 Vir	AEOS	$I$	3.98	mag	2	Roberts et al. (2011)a
	2MASS	$J$	3.80	mag	2	Pinfield et al. (2006)
$\upsilon$ And	AEOS	$I$	3.35	mag	2	Roberts et al. (2011)a
	Lick / LIRC II	$K^{\prime }$	2.86	$\Delta$mag	0	Patience et al. (2002)
BD+24 4697	Gemini North / NIRI	$Ks$	$7.474\pm 0.023$	$\Delta$mag	0	This paper
CI Tau	Subaru / HiCIAO	$H$	$8.43\pm 0.04$	$\Delta$mag	0	Uyama et al. (2017)
EPIC 219388192	Keck / NIRC2	$K$	$10.666\pm 0.0216$	$\Delta$mag	0	Curtis et al. (private communication)
	Subaru / IRCS	$H$	$10.734\pm 0.021$	flux ratio	0	Nowak et al. (2017)
HAT-P-2	Calar Alto / AstraLux	$z^{\prime }$	$9.506\pm 0.001$	$\Delta$mag	1	Bergfors et al. (2013)
HAT-P-20	Keck / NIRC2	$K$	$8.601\pm 0.019$	$\Delta$mag	1	Ngo et al. (2015)
HD 5891	Calar Alto / AstraLux	$i^{\prime }$	$7.47\pm 0.01$	$\Delta$mag	0	Ginski et al. (2016)
HD 33564	Calar Alto / AstraLux	$i^{\prime }$	4.6	$\Delta$mag	1	Ginski et al. (2012)
HD 41004 B	Hipparcos	$Hp$	$8.785\pm {0.014}^{\ast }$	$\Delta$mag	2	Hipparcos and Tycho Catalogues
HD 77065	Gemini North / NIRI	$Ks$	$6.638\pm 0.020$	$\Delta$mag	0	This paper
HD 87646 A	Hipparcos	$Hp$	$8.203\pm 0.002$	$\Delta$mag	2	Hipparcos and Tycho Catalogues
HD 89744	AEOS	$I$	5.2	mag	2	Roberts et al. (2011)a
	UKIRT / UFTI	$H$	4.53	mag	0	Mugrauer et al. (2004)
HD 104985	Calar Alto / AstraLux	$i^{\prime }$	$8.302\pm 0.145$	$\Delta$mag	1	Ginski et al. (2012)
HD 114762	Keck / NIRC2	$K$	$5.888\pm 0.017$	$\Delta$mag	0	Patience et al. (2002)
HD 134113	WIYN / NESSI	$z^{\prime }$	7.6	$\Delta$mag	0	This paper
HD 156279	Calar Alto /AstraLux	$i^{\prime }$	$7.65\pm 0.03$	$\Delta$mag	0	Ginski et al. (2016)
HD 160508	WIYN / NESSI	$z^{\prime }$	7.4	$\Delta$mag	0	This paper
HD 162020	VLT / NACO	$L^{\prime }$	$6.539\pm 0.024$	$\Delta$mag	0	This paper
HD 168443	VLT / SPHERE	$H$	$5.325\pm 0.016$	$\Delta$mag	0	Moutou et al. (2017)
	VLT / SPHERE	$Ks$	$5.211\pm 0.015$	$\Delta$mag	0	Moutou et al. (2017)
HD 178911 B	Hipparcos	$Hp$	$6.835\pm {0.013}^{\ast }$	$\Delta$mag	2	Hipparcos and Tycho Catalogues
HD 180314	Calar Alto /AstraLux	$i^{\prime }$	$6.14\pm 0.05$	$\Delta$mag	0	Ginski et al. (2016)
HD 203949	VLT / SPHERE	$H$	$3.107\pm 0.200$	$\Delta$mag	0	Moutou et al. (2017)
	VLT / SPHERE	$Ks$	$2.994\pm 0.232$	$\Delta$mag	0	Moutou et al. (2017)
KELT-1	Keck / NIRC2	$K^{\prime }$	$9.437\pm 0.019$	$\Delta$mag	2	Siverd et al. (2012)b
Kepler-13 A	Palomar HALE / PHARO	$K^{\prime }$	9.958	$\Delta$mag	0	Adams et al. (2012)
WASP-14	Keck / NIRC2	$K$	$8.621\pm 0.019$	$\Delta$mag	1	Ngo et al. (2015)
WASP-18	VLT / NACO	$L^{\prime }$	$8.131\pm 0.027$	$\Delta$mag	0	This paper
XO-3	Calar Alto / AstraLux	$z^{\prime }$	9.798	$\Delta$mag	0	Bergfors et al. (2013)
Notes. Curve flags: 0 is contrast curve specific to the observations of the target; 1 is average limits of observed sample (or subset); 2 is typical sensitivity of instrument for the observational set up used (see text for more detail). a we used the typical performance curves given in Turner et al. (2006) for the same observational set up. b as no limits are provided in Siverd et al. (2012), we assumed a similar performance as in Ngo et al. (2015) which used a comparable observing strategy on the same instrument, and used the average detection level from that paper. ∗ magnitude of the primary (see text).

Table 4. Table 4: Summary of multiplicity properties from our study, with a comparison to field stars and hosts to lower-mass planets.

Sample	$N$	$M_2$	$\mu$	$\sigma$	$\gamma$	$f$	Binary sep.	Binary $q$	Reference
		(MJup)				(%)	(AU)
Full sample	38	7$-$70	${2.39}_{-0.15}^{+0.14}$	${0.68}_{-0.10}^{+0.12}$	$-{0.52}_{-0.32}^{+0.31}$	${79.0}_{-14.7}^{+13.2}$	20$-$10,000	0.05$-$1	This work
						$70\pm 10$	50$-$2000	0.05$-$1	This work
10 day subset	12	7$-$70	${2.36}_{-0.31}^{+0.28}$	${0.76}_{-0.16}^{+0.29}$	${0.03}_{-0.40}^{+0.38}$	${92.0}_{-19.0}^{+8.0}$	20$-$10,000	0.05$-$1	This work
$>$ 10 day subset	26	7$-$70	${2.40}_{-0.17}^{+0.18}$	${0.63}_{-0.12}^{+0.15}$	$-{0.89}_{-0.64}^{+0.55}$	${74.0}_{-15.9}^{+14.4}$	20$-$10,000	0.05$-$1	This work
FGK field stars	454	$-$	1.70	1.68	$\sim$0	$44\pm 2$	overall	overall	Raghavan et al. (2010)
						$36\pm 2$	20$-$10,000	overall	Scaled in this work
						$16\pm 1$	50$-$2000	overall	Scaled in Ngo et al. (2016)
Friends of HJs	77	4*	…	…	…	$47\pm 7$	50$-$2000	0.05$-$1	Ngo et al. (2016)
*5 objects from the Friends of Hot Jupiters survey have masses between 7$-$12 MJup, all of which are part of our studied sample (see text).

Equations

$$P(\Theta ,m)=1-e^{-\pi {\Theta }^2\rho (m)}.$$

$$P(\rho \mid \mu ,\sigma )\propto \exp [-{\left({\log }_{10}(\rho )-\mu \right)}^2/2{\sigma }^2]$$

$$P(q\mid\gamma)\propto\left\{\begin{array}[]{ll}q^{\gamma}\;\;\mathrm{for}\>\>\gamma\geqslant 0\\ (1-q)^{-\gamma}\;\;\mathrm{for}\>\>\gamma<0\end{array}\right.$$

Full text

††pubyear: 2019††pagerange: A high binary fraction for the most massive close-in giant planets and brown dwarf desert members–B.2

1 Introduction

In the search for analogues to the planets in our own Solar System, exoplanet studies originally firmly excluded known binary systems, despite the fact that about half of Solar-type stars are found in multiple systems (Raghavan et al., 2010). Serendipitous discoveries and subsequent dedicated surveys later revealed that a significant fraction of exoplanets are actually found in binary star systems (e.g. Patience et al., 2002; Desidera et al., 2004; Mugrauer et al., 2006; Mugrauer & Neuhäuser, 2009), mostly with binary separations of at least a few hundred AU. These findings led to numerous studies investigating how stellar binarity affects planet formation and the characteristics and demographics of planetary populations (e.g. Desidera & Barbieri, 2007; Eggenberger et al., 2007, 2011; Daemgen et al., 2009; Adams et al., 2012, 2013; Ginski et al., 2012). The dominant results that emerged from these surveys were a strong deficit of binary companions within $\sim$50$-$100 AU for planet hosts (Roell et al., 2012; Bergfors et al., 2013; Wang et al., 2014a, b; Kraus et al., 2016), and the idea that massive short-period planets appear to be preferentially found in multiple-star systems (Zucker & Mazeh, 2002; Eggenberger et al., 2004).

These studies, however, focused primarily on systems in which the planet had a mass less than $\sim$4 MJup. Theoretical calculations (Kratter et al., 2010; Forgan & Rice, 2011) and numerical simulations (Stamatellos & Whitworth, 2008; Stamatellos, 2013; Hall et al., 2017) both suggest that planets that form via disc fragmentation in gravitationally unstable discs (Boss, 1998) typically have masses above $\sim$4 MJup. Therefore the planets in these existing studies probably formed via the standard core accretion scenario (CA; Pollack et al., 1996), rather than via gravitational instability (GI).

When it comes to planets that formed via core accretion, binarity on close separations is generally considered to have a negative influence (see Thebault & Haghighipour, 2015 for a review of planet formation in binaries and the issues introduced by the presence of a close binary companion). Theoretical studies have concluded that close stellar companions can hinder planet formation by stirring up protoplanetary discs (e.g. Mayer et al., 2005), tidally truncating the discs (e.g. Pichardo et al., 2005; Kraus et al., 2012), or leading to the ejection of planets (Kaib et al., 2013; Zuckerman, 2014). More compact, truncated discs generally have just enough mass left to form a low-mass Jovian planet (Jang-Condell et al., 2008), and planet formation is then further complicated by the very short lifetime ($\lesssim$1 Myr) of these truncated discs (Kraus et al., 2012).

On the other hand, Batygin et al. (2011) and Rafikov (2013) predict that stellar companions should have little influence on planetesimal growth. It has also been proposed that the presence of an outer companion could raise spiral arms in protoplanetary discs, creating regions of high gas and particle densities, favourable for planetesimal formation (Youdin & Goodman, 2005) and pebble accretion (Johansen et al., 2007; Lambrechts & Johansen, 2014). For example, the spiral arm structures observed in the disc around HD 100453 (Wagner et al., 2015) may be due to perturbations from the M-dwarf companion (Dong et al., 2016), located at 120 AU from the primary and originally reported by Chen et al. (2006). Similarly, the asymmetric disc of HD 141569 is attributed to the stellar companions in this triple system (Augereau & Papaloizou, 2004). In the “Friends of hot Jupiters” series of papers, Knutson et al. (2014), Piskorz et al. (2015) and Ngo et al. (2015, 2016) found a binary fraction 3 times higher for hosts to hot Jupiters (mostly up to $M_{p}\sim 4$ MJup) than for field stars on separations of 50$-$2000 AU, and concluded that wide binary systems may either facilitate the formation of Jovian planets, or help the inward migration of planets formed at wider separations.

It has also been suggested that binary companions could induce the inward migration of planets through secular interactions such as the Kozai-Lidov mechanism (Kozai, 1962; Lidov, 1962). In this scenario, an outer companion with a large mutual inclination between the planetary and binary orbits can excite large periodic oscillations of the eccentricity and inclination of the planet. Tidal interactions between the planet and its host star can then drive the planet onto a final orbit with a very small orbital separation when compared to its initial location (e.g. Fabrycky & Tremaine, 2007; Naoz et al., 2012; Petrovich, 2015). In particular, the Kozai-Lidov migration process has been proposed to explain the high obliquities often observed in hot Jupiters, although recent studies indicate that this mechanism can only account for about 20$-$30% of the observed hot Jupiter population (Naoz et al., 2012; Ngo et al., 2016). Similarly, it has been suggested (Rice et al., 2015) that Kozai-Lidov oscillations could drive planetary-mass bodies that form on wide orbits via GI onto short-period orbits. Since disc fragmentation preferentially forms massive planets or brown dwarfs (Kratter et al., 2010; Stamatellos & Whitworth, 2008; Forgan & Rice, 2011), such a process would tend to be associated with more massive planets ($M_{p}\gtrsim 4$ MJup), or brown dwarfs.

Although the true influence of binarity on planet formation and evolution is still unclear, systems hosting planets with masses up to a few Jupiter masses have been extensively surveyed. In contrast, systems with more massive planets ($M_{p}>4$ MJup), and objects within the brown dwarf mass regime, have yet to be studied in detail in the context of stellar multiplicity. Zucker & Mazeh (2002) were the first to point out that massive ($M_{p}>4$ MJup) short-period planets tend to be predominantly found orbiting one component of a multiple star system and possess distinctive characteristics compared to planets orbiting single stars (Eggenberger et al., 2004).

Such massive planetary and substellar companions are very challenging to form at small separations. Giant planet formation, whether by CA or GI, is thought to occur preferentially in the relatively cool outer regions of protoplanetary discs, from a few AU for CA (Pollack et al., 1996), to several tens of AU for GI (Rafikov, 2005). Massive hot Jupiters are thus expected to have formed at wide orbital separations from their host stars, where the lower temperatures in the protoplanetary disc allow for planet formation to proceed (Bell et al., 1997), or be born under very different conditions than currently encompassed by most planet formation models. Recently, Schlaufman (2018) found evidence for two distinct populations of close-in giant planets, with a suggestion of a transition between CA and GI companions occurring at around $\sim$4$-$10 MJup. This is consistent with both semi-analytic (Kratter et al., 2010; Forgan & Rice, 2011) and numerical simulations (Stamatellos & Whitworth, 2008; Stamatellos, 2013; Hall et al., 2017) which suggest that objects that form via GI have masses above $\sim$3$-$5 MJup, and might suggest that these more massive close-in planets formed by GI rather than via CA.

In this work, we aim to constrain the multiplicity statistics of hosts to the most massive giant planets ($M_{p}>7$ MJup) and brown dwarfs found within $\sim$1 AU, in order to investigate the role of binarity in the formation or evolution of these systems. This will allow us to assess if a wide binary companion could be responsible for the observed orbital configurations of these objects, which are both scarce and challenging to explain with current formation theories. Our investigation will provide an indication of whether the Kozai-Lidov mechanism could play a role in the origin of the most massive short-period gas giant planets and brown dwarfs. This study will also help us determine if these massive companions are an extension of the population of lower-mass, CA giant planets, or if they belong to a separate population formed through a distinct mechanism (i.e. GI on wide orbits, followed by inward migration; Nayakshin, 2010; Rice et al., 2015). In particular, we will explore the binary properties of hosts to members of the “brown dwarf desert” (Marcy & Butler, 2000), depicting the significant deficit of brown dwarf companions found within a few AU around Sun-like stars (e.g. Grether & Lineweaver, 2006; Ma & Ge, 2014).

In Section 2 we present our selected sample of targets. Section 3 describes the direct imaging observations acquired for this project and the data reduction. In Section 4 we detail our search for wide companions, using past imaging surveys found in the literature to complement our new direct imaging data, as well as the Gaia Data Release 2 (GR2; Gaia Collaboration et al., 2016, 2018) catalogue. Section 6 describes our approach to the statistical analysis of our survey, and we present our results in Section 7. Finally, we discuss our interpretation of the obtained results in Section 8 and summarise the main results of our project in Section 9.

2 Sample selection

The aim of this project is to search for wide, substellar or stellar companions to stars hosting a massive planet or brown dwarf on a very short orbit. Recent findings suggest that GI forms planets with masses larger than $\sim$4 MJup (Stamatellos, 2013; Stamatellos & Herczeg, 2015; Hall et al., 2017), and the transition between CA and GI companions is thought to occur around $\sim$4$-$10 MJup (Schlaufman, 2018). Studies of core accretion populations found that CA rarely forms planets with masses larger than $\sim$5 MJup (Matsuo et al., 2007; Mordasini, 2018), and shows a steep drop and a strong metallicity dependence in the formation of higher-mass planets (Mordasini et al., 2012; Jenkins et al., 2017). In order to investigate the higher-mass planetary population, which likely formed by disc GI rather than through CA, we choose for this survey a lower limit on inner companion mass $M_{2}$ of 7 MJup, based on the studies mentioned above. This allows us to avoid the main region of overlap between CA and GI, while keeping a sufficiently large sample size for our study. We place an upper limit of 70 MJup (the hydrogen-burning limit) on the mass (or projected mass) of the inner companions, so as to limit our sample to likely substellar objects.

We place an upper limit of $P<400$ days (about 1 AU around a Sun-like star) on the orbital period of the close-in companions. It is now well-accepted that if planet formation via GI does occur, it typically takes place in the outer regions ($>$30 AU) of protostellar discs (Rafikov, 2005; Clarke & Lodato, 2009; Rice & Armitage, 2009). This thus ensures that all selected companions have undergone significant migration between their expected GI formation location and their current configurations, or that they had to be formed under considerably different natal environments than for standard planet formation in order to be born in-situ. We set an upper limit of 500 pc on the distances of our targets in order to be sensitive to wide companions from 50$-$100 AU around most stars in our sample. We use the distance estimates from Bailer-Jones et al. (2018) to infer distances for our targets. These distances are derived from the highly-precise parallax measurements provided by the Gaia DR2 catalogue, correcting for the nonlinearity of the transformation between parallax and distance. Finally, we only consider stellar primaries and place a limit on the host’s mass of $M_{*}>0.1$ M⊙.

Based on the arguments presented above, we selected all systems from the NASA Exoplanet Archive111https://exoplanetarchive.ipac.caltech.edu, the Exoplanet Data Explorer222http://exoplanets.org and the Extrasolar Planets Encyclopaedia333www.exoplanet.eu with confirmed transiting or radial velocity companions with well-constrained orbits that satisfy the following criteria:

$-$ inner companion mass $M_{2}$ (or $M_{2}\sin i$) between 7$-$70 MJup.

$-$ inner companion orbital period $P<400$ days.

$-$ distance within 500 pc based on Gaia DR2 parallax.

$-$ primary mass $M_{*}>0.1$ M⊙.

Our final sample consists of 38 objects, and includes some very short period ($P<10$ days) transiting systems, together with radial velocity objects extending to larger separations. Properties of the inner companions are presented in Table LABEL:t:planet_properties and the host stars are listed in Table LABEL:t:stellar_properties. We selected our sample without regard to the targets’ multiplicity, known or unknown. However, radial velocity and transit surveys are typically biased against binaries, excluding known multiple systems in target selection processes. As these biases are difficult to quantify and account for, our obtained results may somewhat underestimate the multiplicity rate of the population probed here, but we consider that our study is in no way biased towards the presence of wide companions.

About three quarters of the selected systems were discovered and characterised via radial velocity measurements. Mass estimates for companions discovered through this method only allow for the determination of a lower limit on the companion mass due to the unknown inclination $i$ of the system. Radial velocity systems are therefore expected to be more massive than the estimated $M_{2}\sin i$ as a result of the projection factor. Selected systems discovered via this method are thus likely to be more massive than the minimum masses reported in Table LABEL:t:planet_properties. Given the projected masses of these companions and assuming a uniform distribution of inclinations between 0 and 90 degrees, we can easily show with a Monte-Carlo approach that an average of 72% of the radial velocity systems considered here are statistically likely to be above the deuterium burning limit at 13 MJup. Combining this with our transiting systems, this means that more than $\sim$60% of our targets are likely in the brown dwarf mass regime, and close to 80% of our sample is expected to have a true mass $>$10 MJup. We therefore consider that our gathered sample of objects is representative of the population of the most massive planets found on tight orbits and provides a robust insight into short-period brown dwarf desert members.

We define in Table LABEL:t:planet_properties a tidal circularisation timescale $\tau_{\rm circ}$ for each planet and brown dwarf companion in our sample, given in $\log_{10}$[yr]. We estimate this parameter using the formalism presented in Rice et al. (2012), which is based on that developed by Eggleton et al. (1998) (see also Mardling & Lin, 2002 and Dobbs-Dixon et al., 2004). We assume the star has a tidal quality factor of $Q\textquoteright_{*}=5\times 10^{6}$ and the planet has a tidal quality factor of $Q\textquoteright_{p}=10^{5}$. We take the star mass, planet mass, orbital semi-major axis, and orbital eccentricity from Tables LABEL:t:planet_properties and LABEL:t:stellar_properties. We assume that the star has a rotation period of 20 days and that the planet is rotating synchronously. We estimate the circularisation timescale by simply evolving each system for a short period of time and determining the resulting change in eccentricity (i.e., $\tau_{\rm circ}=e/\dot{e}$).

Petrovich (2015) found that planets migrating via the Kozai-Lidov mechanism (Kozai, 1962; Lidov, 1962), under the influence of a distant companion, spend most of their lifetimes undergoing eccentric oscillations at separations $>$2 AU, or as Hot Jupiters at $<$0.1 AU. All the inner companions considered here have orbital separations smaller than 1 AU. If they migrated from wider separations to their current locations through the Kozai-Lidov scenario, they should be able to circularise onto Hot Jupiter orbits fairly rapidly. Inner companions with circularisation timescales longer than the age of the Universe are thus unlikely to be driven by secular perturbations such as the Kozai-Lidov mechanism. On the other hand, objects with timescales smaller than the age of the Universe (i.e., less than $\sim$10.2 in $\log_{10}$[yr]) could have migrated inwards via the Kozai-Lidov scenario. A total of 12 targets have tidal circularisation timescales shorter than that and may thus be consistent with a Kozai-Lidov migration process. The subset of Kozai-consistent objects corresponds to all the inner companions in our sample with an orbital period shorter than 10 days. This is in good agreement with the idea that planets migrating via the Kozai-Lidov mechanism spend most of their lifetime around their initial, wide separations, or on hot Jupiter orbits, as discussed in Petrovich (2015).

3 New Observations

3.1 Observations and data reduction

We used direct imaging facilities at the Very Large Telescope (VLT), Gemini North and the WIYN Observatory to acquire data for six objects in the sample presented in Section 2, four of which did not have any previously reported direct imaging observations. The observations are summarised in Table 1.

3.1.1 VLT / NACO observations

We obtained images in the $L^{\prime}$ filter ($3.8\,\mu$m) using the AO-assisted imager NACO at VLT (Lenzen et al., 2003; Rousset et al., 2003) for HD 162020 and WASP-18 (programme 099.C-0728, PI Fontanive). These new data were acquired with the aim to confirm or refute a candidate reported in Eggenberger et al. (2007) around the former target, and to achieve deeper detection limits than in currently available imaging data of the latter object (Ngo et al. 2015; see Appendix A). The observing setup included the L27 camera, and the data were taken in the pupil tracking mode, where the telescope pupil is held fixed, and the field rotates. Each target was observed using a three-point dither pattern, designed to avoid a bad quadrant of the NACO detector. We used short integration time (0.2 s) in order not to saturate the primaries, allowing photometric and astrometric calibrations.

Standard near-infrared data reduction techniques were applied using our custom IDL routines, including sky subtraction, flat-fielding and bad-pixel correction. Some of the frames were affected by the horizontal additive noise pattern, that sporadically appeared in the NACO data, and was variable in intensity and time. The pattern was removed following the procedure described in Hußmann et al. (2012). Individual frames were de-rotated according to the parallactic angle, and finally stacked together. We retrieved in our final images the unconfirmed candidate companion around HD 162020 reported by Eggenberger et al. (2007) and were able to refute the bound nature of this source based on our new data. The detailed analysis of the rejected candidate is presented in Section 3.3. No companion was detected around WASP-18 within the field of view of our images.

3.1.2 Gemini North / NIRI observations

We acquired images in $K_{s}$ band (1.95$-$2.30 $\mu$m) using the Gemini Near-Infrared Imager (NIRI; Hodapp et al., 2003) instrument at the Gemini North telescope for BD+24 4697 and HD 77065 (programme GN-2017B-Q-40, PI Fontanive). Targets were observed in the standard imaging mode, using the Gemini North adaptive optics (AO) system ALTAIR (Herriot et al., 2000) to obtain diffraction-limited images with the f/32 camera. Both our target were bright enough to be used as natural guide stars. The observing strategy adopted was similar to the one described in Lafrenière et al. (2008) and Daemgen et al. (2015). Each target was observed at five dither positions to allow for sky subtraction and bad pixel correction. At each dither position we acquired one non-saturated short exposure (divided into many coadds) in high read noise mode, followed by a longer exposure in low read noise mode. This prevents our observations from being limited by the high read out noise, resulting in a high observing efficiency and a large dynamic range, providing sensitivity at both small and large separations. Our targets were not saturated, even in the deeper exposures.

We followed standard procedures for near-infrared data reduction, using the Gemini NIRI IRAF package and our dedicated IDL routines. A sky frame was constructed by taking the median of the dithered images, masking the regions dominated by the target’s signal. The individual images were then sky-subtracted and divided by a normalised flat field, and bad pixels were replaced by a median calculated over their good neighbours. For all images, field distortion was corrected as described in Lafrenière et al. (2014). No candidate companion was identified around either target.

3.1.3 WIYN / NESSI observations

We acquired observations of HD 160508 and HD 134113 with the WIYN 3.5-m telescope at Kitt Peak National Observatory (KPNO). We used the NASA-NSF Exoplanet Observational Research (NN-Explore) Exoplanet and Stellar Speckle Imager (NESSI) in diffraction-limited speckle imaging mode. NESSI is based on an upgraded design of the Differential Speckle Survey Instrument (DSSI; Horch et al., 2009, 2012). Each target was observed simultaneously in two cameras, with a filter centered on 562 nm ($r$-narrow) on the blue channel and a bandpass at 832 nm ($z$-narrow) on the red channel. The standard NESSI observing strategy was followed, with typical integration times of 40 ms (see Scott et al., 2018). Data were reduced by the KPNO speckle reduction pipeline that generates reconstructed images and contrast limit curves for each observation (Scott et al., 2018). We did not identify any candidate companion in the obtained data around these two targets.

3.2 Achieved sensitivities

We estimated the limiting sensitivities reached around our observed targets in order to establish the full range of detectable companions covered by the obtained data. For the VLT / NACO and Gemini North / NIRI data, detection limits were determined from the final images described above. The 5-$\sigma$ noise curves were calculated as a function of radius by computing the standard deviation in circular annuli with 1-pixel widths, centred on the primary targets. Noise levels were then converted into magnitude contrasts by dividing by the peak pixel values of the targets (which were not saturated), and converting the obtained flux ratios into magnitude differences in the considered filters. The contrast curves generated from the custom KPNO pipeline in the 832 nm filter were considered for the WIYN data, as the redder filter is better suited for the detection of warm, low-mass companions. The achieved magnitude contrasts are presented in Figures 1, 2 and 3 for our NACO, NIRI and NESSI observations, respectively. The hydrogen-burning limits are shown for the first two data sets, showing that we are sensitive to substellar companions around these stars. We did not reach the stellar/substellar boundary in the WIYN observations and are only able to detect low-mass stellar companions in these data.

3.3 Refuted candidate around HD 162020

HD 162020 had previously been observed with NACO as part of the survey conducted in Eggenberger et al. (2007). Eggenberger et al. (2007) reported two point sources within 5″from the star, one of which was found by the authors to clearly be a background source. They found that the second candidate, at 4$\aas@@fstack{\prime\prime}$98 $\pm$ 0$\aas@@fstack{\prime\prime}$03, was more likely unbound than bound, although the low significance level of this result led them to report the companionship of this candidate as inconclusive based on their data alone. Both sources are retrieved in our new NACO images (Section 3). The positions of the detected sources were extracted using the StarFinder PSF-fitting algorithm (Diolaiti et al., 2000), employing an empirical PSF extracted from the primary.

To calibrate the pixel scale and the True North (TN) of the detector, we used the astrometric calibrator system $\theta_{1}$ Orionis C, observed on October 6 2017. Using the same procedure as described in Chauvin et al. (2012), we obtain the pixel scale of $27.10\pm 0.05\,$mas, and the TN position of $-0.45\pm 0.10\,$deg. However, as previously pointed out by Eggenberger et al. (2007) and Chauvin et al. (2012), additional systematic errors might be present in the determination of the TN of the NACO detector in a case where different sets of calibrator stars were used between the epochs. Since we do not know which calibrators were used to derive astrometry in the previous epochs by Eggenberger et al. (2007), we add $0.5\,$deg to the TN error budget (G. Chauvin, priv. comm.).

With a baseline larger than a decade between the observations used in Eggenberger et al. (2007) and ours, and given the proper motion of the primary from Gaia DR2, we were able to refute the bound nature of this companion. Figure 4 shows the relative positions of the primary and candidate at the various epochs available, together with the expected motion of a background object. The plot clearly demonstrates that the candidate does not share common proper motion with the primary. The fact that the relative positions at the old epochs of observation are not consistent with the primary, nor with the expected background positions, suggests that the source has a non-negligible proper motion of its own.

4 Search for wide companions

We searched for wide companions to the 38 systems from our core sample using our new data and published direct imaging observations, as well as the Gaia DR2 catalogue. A total of 16 objects were found to have at least one wide stellar or substellar companion confirmed to be comoving, listed in Table LABEL:t:confirmed_binaries. Another 7 candidate companions are reported in the literature around 4 of our targets and are presented in Table LABEL:t:candidate_binaries. One of the targets with a reported candidate is already a confirmed wide binary (HD 89744). Figure 5 displays the properties of the inner planets and brown dwarfs, showing the positions in the planet period-mass space of confirmed binaries (star symbols), targets with a candidate companion (triangles) and apparently single objects (circles). In Figure 6 we present the architecture of each identified hierarchical system, plotting the semi-major axes of the inner companions in blue and the projected separations of detected wide binary components in red, with symbol sizes proportional to the planetary masses and binary mass ratios, respectively.

4.1 Literature search and imaging surveys

We conducted an extensive literature search to compile available observations of all objects in our sample and gather existing knowledge about the multiplicity of our targets. We present our findings for each individual target in Appendix A, providing detailed information about every companion, candidate or confirmed, reported around our targets in imaging surveys or catalogues, as well as null-detections. A total of 30 targets are mentioned in the literature in the context of a search for wide companions (with or without detections), to which we add 4 of our 6 observed targets that had no previously reported observations (see Section 3).

Of these 34 objects, we found 21 targets with reported detections in the literature. Among those, 15 are confirmed bound systems (11 Com AB, 30 Ari ABC, $\tau$ Gem AB, $\upsilon$ And AB, AS 205 ABC, HAT-P-20 AB, HD 41004 AB, HD 87646 AB, HD 89744 AB, HD 114762 AB, HD 156846 AB, HD 178911 ABC, Kepler-13 AB, NLTT 41135 AB and WASP-14 AB), which we detail in Appendix A.1. 14 of these binaries or higher-order multiple were demonstrated to form physical pairs in the literature and we confirm the true companionship of the $\tau$ Gem AB system in this work (see Figure 20). We list all confirmed multiples in Table LABEL:t:confirmed_binaries. Three of these systems, HD 87646, HD 41004 and HD 178911, were identified as binaries in the Tycho-Hipparcos catalogues and for the latter two systems, our planet-host stars correspond to the fainter component of the binary system.

The remaining 6 targets are mentioned to have unconfirmed candidate companions and are discussed in Appendices A.2 and A.3. We discarded the 3 point sources reported by Moutou et al. (2017) around HD 168443, which are highly likely to be background contaminants given the crowded galactic latitude of the target (Moutou et al., 2017; see discussion of HD 168443 in Appendix A.3). In Section 3.3, we showed that the faint candidate reported around HD 162020 by Eggenberger et al. (2007) does not share common proper motion with the primary and thus rejected this candidate. We were also able to refute the candidate companion reported around XO-3 by Bergfors et al. (2013) and Wöllert & Brandner (2015) based on the inconsistent parallax and proper motion of this source and XO-3 in Gaia DR2. This leaves 3 targets with unconfirmed candidates, namely, 70 Vir (two candidates), EPIC 219388192 (three candidates) and KELT-1 (one candidate). A candidate companion is also reported by Roberts et al. (2011) around HD 89744, already known to be a confirmed wide binary (Mugrauer et al., 2004). These final 7 candidate companions retained for this study are presented in Table LABEL:t:candidate_binaries.

We are not able to make any clear statement on the physical association of these candidates based on the available data. We can however make a statistical argument on the chance of finding an unrelated background source at close angular separation from the primaries. For each source, we used the Trilegal galaxy models (Vanhollebeke et al., 2009) to calculate a probability of the observed candidates being true companions. This was done by estimating the surface density $\rho$ of background sources expected to be found within 30′from the primary targets, given the galactic latitude and longitude of the objects and the depth and wavelength of the obtained observations. From Brandner et al. (2000), the probability $P(\Theta,m)$ of detecting one or more background stars within an angular separation $\Theta$ (in arcsec) and down to a limiting magnitude $m$ is then given by:

[TABLE]

The probability of an observed candidate being physically associated to the primary is then given by the complement of the chance of alignment, that is, $1-P(\Theta,m)$. The resulting probabilities are listed in Table LABEL:t:candidate_binaries for each candidate. The two faint candidates identified beyond 7″from EPIC 219388192 by Nowak et al. (2017) were found to likely be background sources, with probabilities $<$15% of being physically associated. With the exception of the wider candidate around 70 Vir, most other candidates were found to have very high probabilities of being bonafide companions: the close candidates around 70 Vir, EPIC 219388192, HD 89744 and KELT-1 have $>$99% probabilities of being bound. While additional observations will be required to confirm their true companionship through common proper motion analyses, these objects are therefore highly likely to be true companions.

Finally, a total of 9 targets from our core sample are mentioned as single objects in the literature and have reported null detections from direct imaging surveys (CI Tau, HAT-P-2, HD 5891, HD 33564, HD 104985, HD 156279, HD 180314, HD 203949 and WASP-18; see Appendix A.4). We add to these objects 4 of our observed targets that had no previous observations (BD+24 4697, HD 77065, HD 134113 and HD 160508) and around which we did not find any companion. No published data were found in the literature for the targets 4 UMa, 59 Dra, HD 39392 and HD 112410, which we were not able to observe either.

4.2 Companions in Gaia DR2

Raghavan et al. (2010) found that the period distribution of binary companions to nearby FGK stars is approximately a Gaussian in the logarithm of the period, with a broad peak around 300 yrs ($\sim$50 AU), and a 1-$\sigma$ Gaussian interval spanning from 2 to 1500 AU, in reasonable agreement with previous studies by Duquennoy & Mayor (1991). Most of the imaging data considered here only allow for the detection of companions out to several hundred AU. Hence a significant number of wider companions could remain outside the field of view of these observations and be missed by direct imaging surveys. At these wide separations, outer companions are expected to be massive (i.e. a stellar binary) in order to be able to affect the formation or evolution of close-in planets or brown dwarfs. Such wide stellar companions are expected to be found in the Gaia Data Release 2 (DR2; Gaia Collaboration et al., 2016, 2018), which may thus be used to search for widely-separated comoving components to the objects in our sample. We therefore searched for Gaia DR2 sources with parallaxes and proper motions consistent with those of our targets, to complement our direct imaging search for wide companions in Section 4.1.

The recent release of the Gaia DR2 catalogue provides unprecedentedly-precise astrometric measurements on the parallaxes and proper motions of stars. However, these highly-precise measurements must be considered and handled with caution in the context of a search for common proper motion systems. Shaya & Olling (2011) accurately pointed out that astrometric missions spanning a few years only (e.g. Hipparcos, 3.5 year baseline) capture the reflex motions of multiple systems in their kinematics measurements. Indeed, the components of a binary wobble around the centre of mass of the system and a short-term proper motion measurement is highly likely to reflect this orbital motion. Longer time spans are required to ensure that the observed proper motions correspond to the true barycentric motion (e.g. the Tycho-2 catalogue which uses data from over a century timescale). The apparent changes in proper motions between short and long-term measurements can be as large as several tens of mas yr*-1*, based on the components masses, binary separation, orbital phase and parallax of the system (Shaya & Olling, 2011). These changes in proper motion may even be exploited as a way to search for hidden companions, as was done by Makarov & Kaplan (2005) with the Hipparcos and Tycho-2 catalogues. Despite the excellent precision of Gaia DR2, the catalogue is based entirely on data collected between July 2014 and May 2016, spanning a period of only 22 months, and the same problem as for Hipparcos is encountered.

While these effects are reduced at very wide separations, further complications can also arise from the presence of an unresolved binary. Shaya & Olling (2011) estimated that a tight system separated by a few AU could induce proper motion fluctuations of several mas yr*-1* on the primary, orders of magnitude larger than the errors on Gaia DR2 measurements. Close binaries not resolved in Gaia are treated as single objects in the second Data Release, which can lead to specious astrometric solutions (Arenou et al., 2018). A third component at a wide separation around an unresolved binary is therefore likely to show somewhat different astrometric parameters (proper motion and parallax) compared to its comoving, unresolved primary.

As a result, we adopted rather loose selection criteria to search for comoving companions to the objects in our sample using the Gaia DR2 catalogue. We considered the relative differences in parallax $\pi$, i.e. $\Delta\pi/\pi_{0}\equiv|(\pi_{0}-\pi_{i})/\pi_{0}|$, where the subscript 0 corresponds to our science target and $i$ to other Gaia sources. We then defined similar relative differences for the proper motion components, $\mu_{\alpha*}$ and $\mu_{\delta}$. To account for the uncertainties in the Gaia measurements, we generated, for each pair of objects, $10^{5}$ parallaxes and proper motions drawn from Gaussian distributions centred on the measured values, with a standard deviation set to the Gaia uncertainties. We then calculated $10^{5}$ corresponding fractional differences in $\pi$, $\mu_{\alpha*}$ and $\mu_{\delta}$ and set the final relative differences and associated uncertainties to the mean and standard deviation of the output distributions.

We selected sources that were consistent with relative differences of less than 20% in parallax and in at least one of the two proper motion components (including the correct direction), with a maximum relative discrepancy of 50% in the other proper motion component. We searched for such companions in the Gaia DR2 catalogue for all targets in our sample, out to angular separations corresponding to projected separations of $10^{4}$ AU. We found a total of 11 systems fulfilling the above selection criteria, 9 of which were previously known systems. These systems are listed in Table LABEL:t:confirmed_binaries, in which we give the Gaia DR2 parallaxes and proper motions for each binary component. The characterisation of the two newly-identified Gaia systems, WASP-14 AB-C and WASP-18 A-B, is detailed in Sections 4.2.1 and 4.2.2 below.

In Figure 7 we plot the relative differences in parallax and proper motion (RA and Dec coordinates) between the components of all identified Gaia binaries. The shaded area represents our arbitrary cut at 20% in the relative differences in parallax and proper motion. The 9 previously known systems are marked in blue and our the two new systems discovered here are shown in red. The obtained values and their associated uncertainties are given in Table 2 and are all consistent with our chosen constraints at the 1-$\sigma$ level.

In Appendix B we examined those systems more carefully, as well as other known binaries in our sample, in order to assess our selection criteria. We found that our selection method successfully identified all known binaries that were recoverable given the sensitivity and completeness of Gaia DR2, and consider that we unlikely missed additional binaries present in the Gaia DR2 catalogue. Based on the location of systems with a known unresolved component in Figure 7 (marked with stars), we conclude that most binaries should have relative discrepancies of $<$10% in all astrometric parameters ($\pi$, $\mu_{\alpha*}$ and $\mu_{\delta}$), while systems agreeing to within 20% in parallax and in one of the proper motion coordinates are likely to be hierarchical systems with an unresolved component (see Appendix B for details).

4.2.1 WASP-14

WASP-14 A is already known to have a 0.33 M⊙ bound companion at 300 AU (Wöllert & Brandner, 2015; Ngo et al., 2015) as discussed in Appendix A.1 (see Table LABEL:t:confirmed_binaries). This companion is not detected in Gaia due to the small angular separation (1$\aas@@fstack{\prime\prime}$45) and large magnitude difference ($\Delta$*J$=$*5.2 mag) of the WASP-14 A-B system (see discussion in Appendix B). We report a new companion to this system, WASP-14 C (Gaia DR2 1242084166679297920), at a separation of 11$\aas@@fstack{\prime\prime}$5397$\pm$0$\aas@@fstack{\prime\prime}$0001 and a position angle of 4.5827$\pm$0.0003 deg. The measured angular separation corresponds to a wide projected separation of 1900 AU at the distance of WASP-14 (see Table LABEL:t:confirmed_binaries). WASP-14 AB and C have measured Gaia DR2 parallaxes in excellent agreement, with a relative difference $<$1%. The relative discrepancies in proper motion are slightly larger but still in very good agreement: 4.34% in $\mu_{\alpha*}$ and 11.51% in $\mu_{\delta}$. Given the consistent parallax and small offsets in proper motion, we conclude that the two objects are comoving and form a physical pair. Comparing the placement of WASP-14 in Figure 7 to the other Gaia binaries in our sample also reinforces the idea that WASP-14 is a true binary and confirms our intuition that systems with an unresolved component tend to show larger disparities in their observed short-term proper motions.

WASP-14 C has a Gaia $G$-band magnitude of 17.32 mag, for a magnitude difference of $\Delta$G$=$7.67 mag with the unresolved WASP-14 AB primary. Photometry in the blue ($G_{\mathrm{BP}}$) and red ($G_{\mathrm{RP}}$) filters of Gaia indicate fairly red colours for this object, with $G_{\mathrm{BP}}-G_{\mathrm{RP}}=2.67$ mag. This suggests a mid-K to early-M main sequence star according to the Gaia DR2 HR diagram analysis in Gaia Collaboration et al. (2018). The new companion to WASP-14 is also found in the 2MASS catalogue, with magnitudes of $J=14.297\pm 0.054$, $H=13.801\pm 0.049$ and $K_{s}=13.592\pm 0.058$. According to Schmidt-Kaler (1982), the 2MASS colours correspond to a K5 V spectral type. This implies a bolometric correction of BCK of $\sim$2.3$\pm$0.1 mag (Masana et al., 2006). From these values, we calculated a bolometric luminosity and used the BT-Settl models (Allard et al., 2012) to infer a mass for WASP-14 C. Adopting a distance based on the Gaia DR2 parallax of the target and the age of the system given in Table LABEL:t:stellar_properties, we derived a mass $0.280\pm 0.016$ M⊙ for the newly discovered stellar component of the triple system WASP-14, making this companion the lowest-mass component of the system.

With its low mass and extremely wide separation, WASP-14 C is unlikely to have played a role in the formation or evolution of the 7.8 MJup planet on a 2.2 day orbit around WASP-14 A, as the closer and more massive WASP-14 B component would have had a much stronger influence (if any) on the inner substellar companion. In Figure 8 we show the relative positions of WASP-14 AB and C from the 2MASS and Gaia DR2 catalogues, confirming over a $\sim$20-year baseline that the companion is indeed comoving. The WASP-14 A-B pair is unresolved in both 2MASS and Gaia.

4.2.2 WASP-18

We report here the detection of a faint comoving companion at 26$\aas@@fstack{\prime\prime}$728$\pm$0$\aas@@fstack{\prime\prime}$001 ($\sim$3300 AU projected separation) and a position angle of 200.520$\pm$0.001 deg, found in the Gaia DR2 catalogue (Gaia DR2 4931352153572401152), outside the field of view of our VLT/NACO data (see Section 3). The same companion was independently identified by Csizmadia et al. (2019) at the same time. The source has a Gaia magnitude $G=20.92$ mag, and a $G-G_{\mathrm{RP}}$ colour of 1.79 mag. We do not consider the $G_{\mathrm{BP}}$ photometry of this source because of the known poor quality of fluxes in this bandpass for red sources with G $\sim$ 20$-$21 mag (Gaia Collaboration et al., 2018). At the distance of WASP-18, this places the companion on the M/L-dwarf sequence in the HR diagram for Gaia DR2 sources presented in Gaia Collaboration et al. (2018).

We searched for the same object in the 2MASS catalogue and retrieved a source at the same relative position (26$\aas@@fstack{\prime\prime}$71$\pm$0$\aas@@fstack{\prime\prime}$15 and 200.36$\pm$0.08 deg), in excellent agreement with the consistent proper motions of the two sources in Gaia DR2. Figure 9 shows the relative positions of WASP-18 A and B at the time of the 2MASS and Gaia DR2 observations, together with the expected motion of a background source between the two epochs. The two objects are clearly found to be comoving based on this analysis, confirming our findings from the Gaia DR2 catalogue. The 2MASS photometry of WASP-18 B is $J=16.289\pm 0.096$, $H=15.513\pm 0.083$ and $K_{s}=15.146\pm 0.121$. Using the relation between spectral type and MJ from Filippazzo et al. (2015), we infer a spectral type of M7.5 for the companion, in agreement with our rough estimate from the Gaia HR diagram. We used this spectral type to estimate a bolometric correction BCJ from the relations in Filippazzo et al. (2015) for field objects, and derived a corresponding bolometric luminosity. We finally interpolated the obtained luminosity into the BT-Settl evolutionary models (Allard et al., 2012) at an age of $0.90\pm 0.20$ Gyr (Bonfanti et al., 2016), to obtain a mass of $0.092\pm 0.003$ M⊙ for WASP-18 B.

WASP-18 B is our Gaia source with the largest uncertainties in its parallax and proper motion measurements, with significant errors of 1.52 mas in $\pi$, 1.98 mas in $\mu_{\alpha*}$ and 2.40 mas $\mu_{\delta}$. This is due to the fact that this source only has 169 observations in the Along-Scan direction and none in the Across-Scan direction. In comparison, WASP-14 A has 412 observations in both the Along-scan and Across-scan directions, allowing a much higher precision on its astrometric measurements (0.02$-$0.03 mas). Propagating the measurement errors we found fractional differences in parallax and proper motion between WASP-18 A and B of $16.85\pm 14.60$ %, $6.30\pm 5.94$ % and $10.78\pm 9.03$ %, respectively. While these errors are all very large and comparable to the obtained value for $\Delta\pi$, $\Delta\mu_{\alpha*}$ and $\Delta\mu_{\delta}$, they are still consistent with our selection criteria at the 1-$\sigma$ level. In particular, the 1-$\sigma$ intervals in the fractional difference of both proper motion components remain within 20%, which strongly indicates that both objects are travelling in the same direction. With the close angular proximity of the two sources (26″) and parallax measurements consistent with each other, we conclude that WASP-18 A and B form a physically associated pair.

This is further supported by the red colours of the companion, its placement on the Gaia HR diagram, and the consistent relative positions of the two components over 16 years between 2MASS and Gaia DR2 (Figure 9). The possible offset in the Gaia DR2 parallaxes could indicate that one of the components is an unresolved binary, as we found in our analysis that systems with an unresolved component tend to show higher inconsistencies in the short-term astrometric measurements available in Gaia DR2 (e.g. AS 205 A-BC and HD 178911 AC-B have relative differences in parallax close to 20%). Alternatively, this discrepancy could also reflect a much wider binary separation, along the line of sight, with one component located at a slightly larger distance. Further observations will be required to reduce the uncertainties in the parallax of WASP-18 B, which will be available in future data releases of the Gaia mission. Despite these large uncertainties, we are nevertheless able to confirm that WASP-18 A and B form of a common proper motion pair, in which the fainter component is a low-mass M dwarf close to the stellar/substellar boundary.

5 Detection limits

In order to take into account observational biases and survey sensitivities in our analysis, we gathered and generated detection limits for each target in our sample. We searched for existing contrast curves from the literature for the targets with previous observations (Section 5.1) and derived Gaia detection limits for all targets in our sample (Section 5.2). Combining those with our own sensitivity limits for our six observed targets (see Section 3.2), we are able to define a detection probability map which is presented and described in Section 5.3.

5.1 Imaging contrast curves

The data used to derive detection limits for all targets with existing imaging data (new or from the literature) are summarised in Table 3. Sensitivity limits were found to be available for a total of 29 objects, including our 6 observed targets. When multiple sets of observations were found, we chose the best limits available. For a few targets, the deepest contrast curves found only covered a limited range of separations. In those cases, we also consider the shallower detection limits and keep the best value available at any given projected separation.

We used the contrast curves presented in Section 3.2 for the 6 targets observed as part of this survey. Most targets with archival observations have contrast curves provided in the literature that are specific to the best set of observations for each target. We flag those with a 0 in the “curve flag” column in Table 3. A number of surveys only provide average detection limits for the observed sample (Ginski et al., 2012; Bergfors et al., 2013; Ngo et al., 2015), which we flag as 1. Finally, we considered the typical sensitivities achieved by specific facilities and instruments when no detection limits were available. These curves are flagged with a 2 in Table 3, and are detailed below.

Three of our targets have past observations with the AEOS telescope presented in the survey by Roberts et al. (2011). As the authors do not provide detection limits for their observations, we used the typical performance curves for the AEOS telescope given in Turner et al. (2006) for the same observation set up as that described in Roberts et al. (2011), and consider that they are representative of those data. Similarly, Siverd et al. (2012) acquired NIRC2 images of KELT-1 with Keck but do not present their achieved sensitivities. We therefore assumed a similar performance to that achieved in Ngo et al. (2015) for comparable NIRC2 observations and used the average detection limits from that work for this target.

Pinfield et al. (2006) searched for companions to 70 Vir out to 30″using 2MASS. As no detection limits are available in that paper, we generated a 2MASS contrast curve based on the typical resolving and completeness limits of the 2MASS survey. According to the 2MASS documentation (Skrutskie et al., 2006), close doubles with separations $<$ 5″are not reliably resolved by 2MASS and stellar PSFs can contaminate neighbour sources up to 10″. The $J$-band completeness limit is given at 16.0 mag. We therefore start our 2MASS contrast curve at 5″with $\Delta J=0$ (equal mass binary at the resolving limit). We then use the completeness limit $J=16.0$ from 10″out to 300″, the radial search limit from Pinfield et al. (2006), with a linear increase in $\Delta J$ between 5″$-$10″.

Finally, three targets in our sample were found to be Hipparcos-Tycho binaries, namely, HD 41004 B, HD 87646 A and HD 178911 B. For those three targets, we used the typical sensitivity to binaries in the Hipparcos catalogue based on the plot of separation against $\Delta H_{p}$ of all Hipparcos binaries found in the ESA documentation (ESA, 1997). We extend the separation range out to 30″, the given widest separation of identified Hipparcos-Tycho binaries. For two of these systems, our sample targets correspond to the fainter, lower-mass component of the binary system, which were detected as companions to the brighter primaries. We thus considered the magnitude of the primary of these systems to derive detection limits around the primary binary component.

Most of the obtained detection limits are in units of magnitude difference, $\Delta$mag, while a few are provided in magnitudes and one in flux ratio. These are indicated in the “Limits units” column of Table 3. All sensitivity limits are given as a function of angular separation. For all limits that were not in units of magnitudes we started by converting the contrast curves into apparent magnitudes using the photometry of our targets in the considered filters and given in Table 3. Using the distances from Table LABEL:t:stellar_properties, we then converted all magnitude limits into absolute magnitudes and the angular separations into physical projected separations, in AU. Adopting the ages from Table LABEL:t:stellar_properties for our targets, the obtained absolute magnitude curves were then interpolated into the BT-Settl evolutionary models by Allard et al. (2012) to derive corresponding minimum detectable companion masses. The BT-Settl models provide isochrones for numerous photometric systems. We were therefore able to use models corresponding to the specific facilities and filters considered and infer mass limits for each target. Finally, we used the stellar masses listed in Table LABEL:t:stellar_properties for our sample to convert the obtained mass limits into mass ratios $q$ (using the masses of the binary primaries for the two Hipparcos systems mentioned above). For the few targets with multiple entries in Table 3, we considered the lowest mass ratio value in the overlapping regions of separations in order to define a unique sensitivity curve for each object. The final mass and mass ratio curves for each target with direct imaging data are shown in Figure 10, together with the positions of all confirmed (black stars) and candidate (open circles) companions around these objects.

5.2 Gaia detection limits

Since all objects in our sample are found in the Gaia DR2 catalogue, we are able to derive Gaia detection limits for each of our targets. Gaia DR2 is found to be complete between $G=12$ and $G=17$, with a limiting magnitude of $G\sim 21$ and a bright limit of $G\sim 3$ (Gaia Collaboration et al., 2018). Ziegler et al. (2018) investigated the recoverability of close binaries in Gaia DR2 looking for known binaries from the Robo-AO Kepler survey (Law et al., 2014; Baranec et al., 2016; Ziegler et al., 2017) in the Gaia DR2 catalogue. They found that near equal-brightness binaries ($\Delta G<1$) were consistently retrieved from separations of 1″and that binaries down to $\Delta G=6$ were recovered at separations of $\sim$3″. Based on their results, we define our Gaia DR2 sensitivity limits to start at 1″and $\Delta G=1$, with a linear decrease to $\Delta G=6$ from 1″$-$3″. We then adopt a linear decrease out to 5″from $\Delta G=6$ to $G=21$, the Gaia faint limit, and use that limiting magnitude at wider separations, out to projected separations corresponding to $10^{4}$ AU.

Figure 11 shows the obtained sensitivity limits in terms of apparent $G$ magnitude and mass ratio for all objects in our sample. We plot on the left panels the limits for the targets without a known companion and on the right panels the limits for all confirmed multiple systems, with the positions the known companions. Magnitude limits were converted into corresponding mass ratio curves adopting the properties of our targets listed in Table LABEL:t:stellar_properties and following the approach described in the previous section with BT-Settl isochrones specific to the Gaia filter system.

While Gaia is essentially complete in the range *G$\sim$*12$-$17 mag, the catalogue has an ill-defined faint magnitude limit which depends on celestial position (Gaia Collaboration et al., 2018). In addition, the number of sources with a full 5-parameter astrometric solution (position, parallax and proper motion) decreases towards the faint end, where a larger fraction of sources only have positional measurements available, as discussed in the assessment of the Gaia DR2 astrometric performance by Lindegren et al. (2018). We must take into account the catalogue completeness in our detection limits to account for companions that are missed by Gaia, but also for those like $\tau$ Gem B and HAT-P-20 B that only have a 2-parameter solution and which we were not able to identify as Gaia companions in our analysis in Section 4.2. Arenou et al. (2018) report the completeness of the Gaia DR2 catalogue as a function of $G$-band magnitude in their catalogue validation work. The provided completeness level for sources with full astrometric solutions decreases from $\sim$99% at $G<17$ to $\sim$80% at $G=20$, before sharply dropping to 0% as $G$ approaches 21 mag (see figure A.1. in the appendix of Arenou et al., 2018). We thus use the completeness levels provided in that paper to account for these effects.

In Figure 12 we show an example of Gaia sensitivity curves for a primary of $G=8$ with a parallax of 20 mas, corresponding to a mass of 1 M⊙ at 3 Gyr, representative of the targets in our sample. In the top panel, we show the detection limits for such an object in terms of apparent magnitude, defined as described above. The horizontal dashed lines indicate the $G$ magnitudes associated with various completeness levels using the information from Arenou et al. (2018), down to the faint magnitude limit of Gaia DR2 at $G=21$ (grey line). The bottom panel shows the same contrast curve, converted into mass ratios using the BT-Settl models and adopting an age of 3 Gyr. Since we assumed a primary mass of 1 M⊙, the plot in the bottom panel is also representative of the corresponding mass limits in units of Solar masses. Figure 12 clearly demonstrates that when working in mass ratio space the range over which Gaia DR2 is not complete for sources with 5-parameter solutions ($<$99% completeness, below the red dashed lines) is significantly reduced relative to the span of the same incompleteness levels in magnitude space, going from *G$=$*17$-$21 to *q$=$*0.16$-$0.09 M⊙. This implies that in addition to the targets too faint for Gaia ($G>21$, below the 0% completeness grey dashed lines in Figure 12) only the lowest-mass companions have a high chance to of being missed due to survey incompleteness. We note however that the two known companions from our sample that only have a 2-parameter astrometric solution in Gaia DR2 ($\tau$ Gem B and HAT-P-20 B) have relatively bright $G$-band magnitudes of 9.42 mag and 12.80 mag, respectively. Their apparent magnitudes fall into a $>$99% completeness level of sources with full astrometric solutions according to Arenou et al. (2018). We conclude that it was statistically unlikely to have two bright companions in our sample present in the DR2 catalogue but lacking a 5-parameter solution, and that these two sources are not representative of the completeness of Gaia DR2.

For all targets in our sample, we converted the completeness curve obtained from Arenou et al. (2018) into mass ratios as we did for our example target in the bottom panel of Figure 12. We are thus able to assign a completeness factor to each mass ratio value for every detection limit presented in Figure 11. Instead of traditional sensitivities, where anything above the mass ratio curves is considered as detectable and anything below is not retrievable, we now associate every point in the separation-mass ratio space to a detection probability, given by the completeness level at any given mass ratio value. The part of the parameter space below the final limits remains at the zero detection probability level regardless of the associated completeness value. We will use these probabilities in the next section to define a 2-dimensional detection probability map for our sample.

5.3 Detection probability map

We combined all sensitivity curves obtained in Section 5.1 from imaging data and in Section 5.2 from Gaia DR2 to define a single detection probability map for our survey, as was done in Fontanive et al. (2018). For targets with both Gaia and imaging limits, we started by combining the two sets of contrast curves. Following the approach described in Section 5.1, we considered the best value available (lowest $q$ value) in the separation ranges where the Gaia and imaging limits overlapped, keeping track of the ranges over which the final curves corresponded to the Gaia limits. This allowed us to define a unique sensitivity curve for each object.

The mass ratio limits for each target in the sample were then placed on a grid of separations and mass ratios with a resolution of 0.002 in $q$ and steps of 0.01 in log($\rho$). For every cell in the grid, we then identified the number of targets around which a companion of given separation and mass ratio would have been retrieved given the data gathered for this survey. When the considered separation corresponded to the Gaia limit of a given target, we then scaled this detection by the Gaia completeness level at the associated mass ratio, which we previously computed in Section 5.2. The number obtained for each cell of the grid was then divided by the total number of objects in our sample, providing a value between 0 and 1 representing the average probability that a companion of projected separation $\rho$ and mass ratio $q$ would have been detected around our 38 targets given the data considered in this work.

Figure 13 shows the resulting detection probability map for our full sample of 38 objects, considering all available imaging data for the targets in our survey and the Gaia DR2 catalogue. Companions inside the 99% completeness region are essentially detectable around all targets in the sample. We are complete to companions with $q>0.2$ at separations $>$1000 AU around 90% of our sample, and down to $q\sim 0.1$ from separations of $\sim$100 AU around half of our targets (50% detection probability contour). Confirmed comoving companions were found to be relatively evenly distributed throughout the parameter space (both in separation and mass ratio). In contrast, most candidate companions are concentrated around $q\sim 0.1$, which we attribute to the fact that these fainter companions are not detected by Gaia and are more likely to lack a second direct imaging epoch, necessary to be astrometrically confirmed. Interestingly, no companion was found at projected separations $<$20 AU, despite reaching a completeness level up to 50%, while a number of companions (confirmed and candidates) were retrieved at the same detection probability level at wider separations and low mass ratios.

6 Statistical analysis

We used the statistical tool described in Fontanive et al. (2018) to constrain the multiplicity properties of our sample. Examining the binary statistics of these objects will allow us to investigate the possible role of binarity in the formation or evolution of massive, close-in brown dwarfs and planets. The code is based on a Markov Chain Monte Carlo (MCMC) sampling method, using the emcee Python package (Foreman-Mackey et al., 2013), and allows us to place robust Bayesian statistical constraints on the binary frequency and companion population distributions for the sample gathered in this study. We add a new capability to the tool in order to account for unconfirmed candidates, which we describe below.

The statistical approach uses the detection limits of the survey, in the form of a detection probability map (see Section 5.3, Figure 13), and the properties of detected companions (total number of detections, separations and mass ratios of identified systems) to derive posterior distributions of model parameters (binary fraction and parameters describing the shapes of companion distributions) most compatible with the gathered data. Based on previous studies of stellar multiplicity in the field (Duquennoy & Mayor, 1991; Raghavan et al., 2010), we adopt a lognormal distribution in companion separation $\rho$ (Equation 2) and a power-law in mass ratio $q$ (Equation 3):

[TABLE]

where $\mu$ and $\sigma$ are the mean and standard deviation of the underlying normal distribution in $\log(\rho)$. The mass ratio distribution ranges from 0 to 1 and is defined by the power-law index $\gamma$:

[TABLE]

so that negative and positive indices produce symmetric distributions about $q=0.5$ for the same absolute value of $\gamma$.

As was done in Fontanive et al. (2018), we truncated the model distributions at $\rho=$ 20$-$10,000 AU and $q=$ 0.05$-$1, in order to constrain the binary frequency on those separation and mass ratio ranges. We adopted flat priors for each model parameter, set to unity over the following ranges and to zero elsewhere: 0.5$-$4 for $\mu$, 0.1$-$3 for $\sigma$, -3$-$3 for $\gamma$ and 0$-$1 for the binary fraction $f$.

In order to also take into account the candidate companions identified around the targets in our sample, we used the probabilities of the candidates being physically bound derived in Section 4.1 and listed in Table LABEL:t:candidate_binaries. At each step in the MCMC chain, we drew for each target with a candidate companion a number between 0 and 1, and counted the candidate as a bonafide companion if the drawn value was below the companionship probability. This ensures that each candidate companion is selected in a fraction of MCMC steps that is representative of its probability of being physically associated to the primary target. For hierarchical systems (e.g. 30 Ari, WASP-14), we considered the properties of all detected components in the part of the code that constrains the shape of the separation and mass ratio distributions, accounting for candidate companions only when they were selected (e.g. 70 Vir, EPIC 219388192, HD 89744). For systems in which the binary companion is itself a tight binary (AS 205, HD 178911, Kepler-13), we used the combined mass of the binary component, since this total mass would be responsible for any dynamical effect on the close-in planet or brown dwarf. For the section of the tool constraining the binary fraction, we considered the number of multiple systems rather than the total number of companions in order to constrain the multiplicity rate of our sample.

7 Results

7.1 MCMC analysis for the full sample

We ran the MCMC sampler with 2000 walkers taking 5000 steps each on our full sample of 38 objects. We found that walkers were expanding from their initial positions to a reasonable sampling of the parameter space in less than 100 steps, and removed the first 100 steps of this “burn-in” phase. We found a mean fraction of steps accepted for each walker of 0.44, in good agreement with the rule of thumb acceptance fraction suggested by Foreman-Mackey et al. (2013) between 0.2 and 0.5, and trust the obtained value to be a reliable indication of convergence.

The full output from our MCMC analysis is presented in Figure 14. The best-fit values for the binary parameters of our core sample on separations in the range 20$-$10,000 AU are summarised in Table 4. Errors correspond to 68% confidence intervals, estimated using a highest posterior density approach to determine the boundaries of Bayesian credible intervals (see Fontanive et al., 2018). The highest density region method provides a set of the most probable values of a given parameter. All 4 model parameters were found to be well-constrained, converging to sharply-defined peaks in the posterior distributions, with the binary fraction $f$ showing the broadest posterior distribution.

The obtained posterior distribution for the binary frequency of our sample, $f=79.0^{+13.2}_{-14.7}$%, is found to be much higher than the overall multiplicity rate of FGK stars in the field, generally observed to be around 40$-$50% (Duquennoy & Mayor, 1991; Raghavan et al., 2010). The peak of the lognormal in separation corresponds to a value of $\sim$250 AU, also much wider than for the field population (Raghavan et al., 2010). We discuss these features further in Section 8 where we provide a detailed assessment of the possible sources and implications of these results.

7.2 Sample division at 10 days in inner companion period

We divided our sample into two subsets, with a cut at 10 days in the period of the inner planets and brown dwarfs, the commonly-accepted threshold for hot Jupiters (Wang et al., 2015; Dawson & Johnson, 2018). This allows us to investigate possible differences in the binary properties of the stars with companions on orbits comparable to hot Jupiters, and those with planets or brown dwarfs at slightly wider separations. The hot Jupiter-like subset includes 12 targets, 6 of which are confirmed binaries, with 2 additional targets having at least one high-probability candidate companion. The sample of wider inner companions contains 26 objects, including 10 confirmed multiples and 1 candidate binary.

Following the approach described above, we created detection probability maps for each subset, considering the available detection limits for all targets from each subsample. We then performed the same statistical analysis as that presented above to constrain the multiplicity rates and binary properties of our samples of objects with periods shorter and longer than 10 days, so as to assess whether statistically significant discrepancies are observed between the two populations. Figures 16 and 16 show the output of the MCMC sampler for the shorter and longer-period samples, respectively. As expected from the smaller sample sizes of the two subsets relative to the full sample, the walkers are slightly more widely spread throughout the parameter space than in Figure 14, and this effect is amplified for the smaller sample of $<$10 days companions (Figure 16). Nevertheless, all 4 model parameters are still well-constrained within the explored parameter space in both subsamples. The best-fit values and corresponding 1-$\sigma$ intervals are given in Table 4 for each subset.

The model parameters describing the companion separation distribution ($\mu$ and $\sigma$) peak at very similar values for the two subsets, indicating that no significant difference is found in the binary separation of these two populations. The power-law index $\gamma$ describing the mass ratio distribution appears to shift to slightly lower values for the sample of longer-period inner companions, which reflects the generally lower mass ratios of multiple systems found in that subset.

The binary fraction $f$, on the other hand, shows a larger discrepancy in the output posterior distributions. The obtained probability density function for the sample with inner companions on very short orbits (Figure 16) was found to peak at $92.0^{+8.0}_{-19.0}$% (68% confidence), consistent with a binary rate of 100% at the 1-$\sigma$ level. In contrast, the subset of wider inner companions (Figure 16) has a binary frequency of $74.0^{+14.4}_{-15.9}$%. We plot these two distributions in Figure 17, together with the obtained posterior distribution for $f$ for the full sample. The red line shows the corresponding multiplicity rate of field stars based on results from Raghavan et al. (2010), scaled to our probed separation range of 20$-$10,000 AU (see Section 8 for details). While all 3 distributions are consistent with one another, a clear shift is observed in the peak of the posteriors. In particular, the peak of the binary fraction distribution for the $<$10 day sample is located outside or at the edge of the 68% confidence interval of the other two posterior distributions. The resulting binary fractions are all much higher than the corresponding value expected for field stars.

8 Analysis and Discussion

8.1 Comparison with field stars

8.1.1 Multiplicity fraction

Raghavan et al. (2010) provided a comprehensive assessment of the multiplicity properties of Solar-type stars, searching for companions to 454 F6$-$K3 primaries in the field. Taking into account the completeness limits of their survey, the authors found that about $56\pm 2$% of stars are single, for an overall multiplicity fraction of $44\pm 2$%, in good agreement with previous results from Duquennoy & Mayor (1991). Our binary fractions derived in Section 7 were limited to separations in the range 20$-$10,000 AU. We must therefore restrict the overall binary rate from Raghavan et al. (2010) to this separation range in order to compare our findings to the general field population. Taking into account the shape of the distributions obtained by Raghavan et al. (2010) and excluding all companions from that study with separations outside our considered range, the fraction of stars found in binaries or higher-order multiples becomes $36\pm 2$%. This is more than twice as low as the binary rate obtained for our full sample of 38 objects, with a 3-$\sigma$ significance ($f=79.0^{+13.2}_{-14.7}$%; see Table 4). We also find the value for field stars to be lower than the binary rates derived for our two separate subsets, although these results have a lower significance ($\sim$2.5-$\sigma$ level) as a result of the smaller number statistics of the individual subsamples. In Figure 17 we compare the expected fraction of multiples on this separation range to the posterior distributions obtained for the full sample studied in this work and to the two subsets defined in Section 7.2. The plot clearly shows that our derived binary rates are statistically larger than the binary fraction from the overall FGK stellar population.

Zucker & Mazeh (2002) noted a substantial paucity of high-mass planets with short-period orbits around single stars. In contrast, this feature is not observed around binary star systems, which exhibit a prevalence of short-orbit massive planets (Zucker & Mazeh, 2002; Eggenberger et al., 2004; Desidera & Barbieri, 2007; Mugrauer et al., 2007). Our results are highly consistent with these observations, suggesting that almost all stars with a $>$7 MJup companion within $\sim$1 AU are part of multiple stellar systems. The statistically higher binary occurrence of hosts to massive planets relative to the general field population indicates that stellar companions may play an important role in the existence of the most massive giant planets and brown dwarfs observed on tight orbits, and that a binary companion may be required to explain their presence. The much higher binary fractions we find for our sample compared to field stars, despite the known biases from transit and radial velocity surveys against close binaries, reinforces the idea of a significant correlation between stellar binarity and the existence of the massive inner companions studied here. While the nature and magnitude of this role are not clear and cannot be established based on this study alone, a number of possibilities have been formulated and explored in the literature to explain the possible influence of binary companions on giant planet formation and evolution. We discuss these theories in Section 8.3.

A caveat of this analysis is that Raghavan et al. (2010) studied stars in a volume limited to 25 pc in distance. In order to compare our results to the overall field population, we extrapolated the measurements from Raghavan et al. (2010) out to distances of 500 pc. While the distributions found by Raghavan et al. (2010) are valid for 0.5$-$1.5 M⊙ stars, our sample contains two targets from young star-forming regions (one confirmed binary and one not known to have any companion), as well as a number of giants and subgiants. These populations may have different binary statistics than the main sequence solar-type stars probed by Raghavan et al. (2010) and the assumptions made in our analysis may not be entirely valid. The field study by Raghavan et al. (2010) was also heavily biased towards G stars, while our sample contains a number of more massive A and F stars, which are expected to have a higher binary fraction, as well as some M-dwarfs, expected to have a lower binary frequency. The mass dependence of stellar binarity could therefore be another factor affecting our results. The field stars sample may also be contaminated with planet hosts, and it must be pointed out that the results presented above are not a comparison between planet-free stars and planetary hosts, but rather an assessment of planet hosts multiplicity properties relative to the general stellar population. That being said, the extremely high binary fraction derived for our studied sample is still a robust and significant result by itself, even if the comparison to field stars may not be fully reliable.

8.1.2 Mass ratio distribution

Raghavan et al. (2010) found a roughly flat mass ratio distribution for binaries separated by more than 100 days. The value obtained for the power-law index in our full sample indicates a slight preference for lower-mass companions, but is fairly consistent with a flat distribution (i.e. $\gamma=0$, see Table 4). Our subset of $>$10 day inner companions indicates a moderately larger preference towards low mass ratio companions for these systems, and the subsample of short-period planets was found to exhibit a uniform distribution in wide companion mass ratio. Our probed samples are thus in reasonable agreement with the mass ratios observed around multiple stars in the field, and we find no evidence for distinct populations between our studied targets and the general field population.

8.1.3 Separation distribution

In contrast, we found larger and more significant disparities in binary companion separation between the distributions obtained in Section 7 and the expected distribution from FGK field stars, as shown in Table 4. Raghavan et al. (2010) reported a lognormal distribution in companion separation peaking at 1.70 in $\log_{10}(a)$, with a Gaussian width of 1.68, corresponding to a broad peak around 50 AU. This is significantly smaller than our derived value of $\mu=2.39^{+0.14}_{-0.15}$ for our full sample, with a mean located at $\sim$250 AU. We also found a much narrower separation distribution, with a Gaussian width of $\sigma=0.68^{+0.12}_{-0.10}$. The results obtained for our two subsets are in good agreement with each other and with the full sample (Table 4). Figure 18 shows the constraints obtained from our statistical analysis on the separation distribution of the multiple systems in our core sample and defined subsets. The red distribution represents the results obtained by Raghavan et al. (2010) for solar-type stars in the field, clearly demonstrating the preference for wider binaries among our targets and the strong deficit of closely-separated systems in our studied sample.

While we restricted ourselves to a 20$-$10,000 AU separation range to constrain the binary frequency $f$, the parameters describing the separation distribution ($\mu$ and $\sigma$) were explored over a broad parameter space. The MCMC walkers would have been able to converge to a distribution peaking near or even inside 20 AU had it been compatible with the observed data. As we noted in Section 5.3, the lack of tight binaries is unlikely to be the sole result of observational biases and limiting sensitivities. Indeed, we are sensitive to binary companions with separations of 5$-$20 AU around 20$-$60% of our targets depending on the system mass ratio, and a number of companions (confirmed and candidates) are retrieved at the same detection probability level at larger separations and low-mass ratios (see Figure 13). If the true underlying separation distribution in our sample was comparable to that of field stars, with a broad peak near 50 AU, we would have expected to detect numerous close binaries given the number of wide multiple systems already present in the sample. We thus consider that this observed feature is real and not due to observational limitations.

This deficiency of tight binaries could however be due to selection effects in exoplanet surveys, which are often biased against close binary systems. In addition, it may be harder to detect planets via the radial velocity method in the presence of a close, massive stellar companion. The spurious assumption that a planet host is single when it is in fact an unresolved binary will also lead to erroneous measurements of the planet’s physical properties. A population of massive planets and brown dwarfs in closely-separated binaries could hence exist and be underrepresented or misreported among detected exoplanets. If this is the case, the true multiplicity rate of systems hosting massive, close-in planets or brown dwarfs should then be even higher than what we found here.

Nevertheless, our observations are consistent with previous studies. The shortfall of close binaries among planet hosts has indeed been vastly reported in the literature (Roell et al., 2012; Bergfors et al., 2013; Wang et al., 2014a, b; Kraus et al., 2016) and is generally attributed to a hindrance of planet formation in very tight binaries. While observational constraints remain sparse, the current census is that binarity on scales comparable to the Solar system ($\lesssim$50$-$100 AU) has the potential to affect planet formation and evolution. However, different conclusions have been reached on the theoretical side, and it is not clear whether these processes are altered, inhibited or facilitated by the presence of a binary companion and on what separations these effects may take place. This is further discussed in Section 8.3.

Another possible explanation for the depletion with tight binary systems among our sample is that the inner planetary and brown dwarf companions to our targets may have formed at much wider separations than their current locations, at radii overlapping with the missing population of binary companions. Such massive planets are indeed expected to form outside at least a few AU for formation by CA (Mordasini et al., 2012), and more likely several tens of AU for GI, in regions of the circumstellar discs that are massive and cool enough to form massive giant planets (Rafikov, 2005). If inner companions have formed at such separations, additional, massive binary companions should not exist within a few tens to hundreds of AU around these systems, which would be reconcilable with our observations. However, as this trend is also observed around systems hosting low-mass planets, for which a formation at very wide separations is not required, this may not be the primary phenomenon responsible for this feature.

8.2 Binarity as a function of planet properties

8.2.1 Binary frequency versus inner companion period

In Section 7.2, we divided our sample into two subsets in order to investigate possible differences in the demographics of stars hosting planetary or brown dwarf companions within and beyond orbital periods of 10 days. While we found no evidence for distinct binary mass ratio or separation distributions between these two populations, our statistical analysis revealed a possibly larger binary frequency for the subset of shorter-period companions, with a peak at 92%, compared to 74% for the subsample of more widely-separated systems. These results are marginal due to the less stringent constraints we were able to place on the individual subsets, as shown by the broader posterior distributions in Figure 17 relative the one obtained for the full sample. Larger sample sizes will be required to confirm this tendency.

This theory is nonetheless supported by the similar trend seen for hosts to lower-mass planets. Surveys searching for wide companions to radial velocity planets from sub-Jupiter masses up to a few MJup and out to 5 AU found that less than $\sim$25% of these systems were part of binaries or multiple systems (e.g. Raghavan et al., 2006; Ginski et al., 2012), although these surveys may be biased or incomplete. In contrast, studies of slightly shorter-period transiting planets (typically $<$100 days) observed rather higher binary fractions, generally around $\sim$50%, for planets of comparable masses (Adams et al., 2012, 2013; Ngo et al., 2016). Furthermore, Tokovinin et al. (2006) found that 96% of spectroscopic binaries with periods $<$3 days have a third component, compared to only 34% of systems with periods longer than 12 days, albeit some selection biases may affect these results to a limited extent.

The marginal difference in binary occurrence observed in this work between very short-period transiting planets and brown dwarfs, and the somewhat wider population of radial velocity companions, thus indicates that this trend, if real, may also hold for the very massive inner companions studied here. This trend could suggest that binarity greatly helps the formation or migration of massive giant planets and brown dwarfs observed within $\sim$1 AU, and essentially becomes necessary for these companions to reach orbital periods shorter than 10 days.

8.2.2 Binary frequency versus inner companion mass

In the final paper of the Friends of hot Jupiters campaign (Knutson et al., 2014; Piskorz et al., 2015; Ngo et al., 2015), Ngo et al. (2016) found that 47$\pm$7% of hot Jupiter systems have a stellar companion between 50 and 2000 AU, a binary rate 3 times higher that for field stars in this separation range. The authors concluded that binary companions on these separations facilitate planet formation or help the inward migration of giant planets. Our study probed higher-mass planets than those considered in that survey, allowing us to examine trends in stellar multiplicity as a function of planet mass, including inside the brown dwarf regime. The Friends of hot Jupiters survey looked at systems with planet masses mostly limited to 4 MJup. Only five objects with more massive companions were studied in that work, with masses between 7 and 12 MJup, all of which are part of our selected targets (3 are confirmed multiples: HAT-P-20, WASP-14 and WASP-18; 2 are apparently single: HAT-P-2 and XO-3). The binary fraction derived here for more massive objects was estimated for separations between 20 and 10,000 AU. The corresponding binary rate restricted to the 50$-$2000 AU separation range becomes 70$\pm$10% for our core sample, 1.5 times larger than the value from Ngo et al. (2016) at the 2.3-$\sigma$ level. This fraction is 4 times larger for field stars on these separations (see Ngo et al., 2016), with a 5-$\sigma$ significance. These results are summarised in Table 4.

Our results for the shorter-period subset are less significant due to the looser constraints we were able to place on the smaller-sized subsample. We thus only compare previous studies with the binary fraction estimated for the full sample, and keep in mind that hosts to shorter-period planets may have an enhanced binary rate, as discussed previously. Our findings suggest that the trends characterised by Ngo et al. (2016) for hosts to hot Jupiters are also observed and even strengthened for the highest-mass close-in planets and brown dwarfs. These results are in excellent agreement with early observations by Zucker & Mazeh (2002) and Eggenberger et al. (2004), who determined that the most massive planets on orbits of a few days are consistently found in binary systems, suggesting that this planetary population does not exist around single stars.

We note that four targets in our sample have notably high mass estimates (40$-$60 MJup) relative to the rest of our sample, namely BD+24 4697, HD 77065, HD 134114 and HD 160508. These objects appear somewhat isolated in the period-mass parameter space in Figure 5. Given that their mass measurements are lower limits derived from radial velocity information, their true masses are most certainly even higher. Assuming a uniform distribution of inclinations, we may calculate the minimum value for the projected mass that corresponds to a true substellar mass $M_{2}$ $<$ 70 MJup at a given confidence level. This translates to $M_{2}\sin i<34$ MJup for a 68% confidence of a true mass $M_{2}$ below the hydrogen-burning limit. These four targets thus have a large chance of being stellar, and are therefore likely to have formed as tight stellar binary systems, rather than brown dwarf companions forming in a circumstellar disc around the host star. None of these systems were found to have a wide binary companion. Excluding these systems from our survey to focus on substellar companions only would hence have resulted in an even higher binary fraction for the rest of our sample. This further reinforces the idea that the most massive planets and brown dwarfs forming in discs and detected within $\sim$1 AU require a wide stellar companion to form or evolve to their observed orbital configurations.

8.3 Implications for formation and evolution processes

Our results demonstrate a very robust correlation between binary occurrence rate and the sporadic population of close-in massive giant planets and brown dwarf desert inhabitants. Whatever the underlying processes, this concurrence implies that wide binaries must have an influence on the observed population of short-period planetary and brown dwarf companions, which could occur at the stage of formation or during later evolution.

Zucker & Mazeh (2002) were the first to raise the possibility that planets in binaries may have a different mass-period distribution, a trend subsequently confirmed by Eggenberger et al. (2004), Desidera & Barbieri (2007), Mugrauer et al. (2007) and others. Our results are in very good agreement with these studies, suggesting that the most massive planets observed within $\sim$1 AU are almost exclusively found in binary systems, and that this feature is amplified as planets or brown dwarfs reach shorter periods. Desidera & Barbieri (2007) concluded that the presence of a stellar companion on separations $<$100$-$300 AU may be able to modify the formation or evolution of giant planets. Eggenberger et al. (2004) also found that massive planets in binary systems with periods shorter than 40 days have very low eccentricities, suggesting that these planets likely underwent some form of migration, possibly induced or driven by outer binary companions, to evolve to their current orbits. Duchêne (2010) investigated these observational trends and concluded that binarity does not affect the formation and growth of planetesimals (see also Batygin et al., 2011 and Rafikov, 2013). Duchêne (2010) proposed that planet formation in binaries tighter than $\sim$100 AU occurs at a similar rate but through different mechanisms than around wider binaries and single stars, possibly explaining the observed preponderance of very massive, close-in planets found in binaries but rarely seen around isolated stars.

Simulations by Kley (2001) showed that perturbations from a secondary star may alter the formation and evolution of a planet, in particular by enhancing the mass accretion and orbital migration rates. This could explain why the most massive short-period planets are found in multiple systems, the presence of stellar companions enabling massive planets to achieve smaller orbital separations than the corresponding limit for planets orbiting single stars (Eggenberger et al., 2004). Jensen & Akeson (2003) found that the distribution of disc mass in $>$200 AU binary systems among T Tauri stars is not always determined by the stellar masses and may be more asymmetric, with the primary retaining a much more substantial disc and the secondary being left with a very low-mass disc. Massive discs around primaries in wide binaries could thus provide larger reservoirs of material for planet formation, which is thought to be favourable to the formation of higher-mass planets, as discussed in Mordasini et al. (2012). The shorter lifetime of circumstellar discs in tight binaries (e.g. Kraus et al., 2012) argues for a formation via gravitational collapse of the circumstellar disc (thousand year timescale) rather than through core accretion, which requires 1$-$10 million years. A favoured formation by gravitational disc instability is further supported by the very high masses of the giant planets or brown dwarfs considered here. Furthermore, theoretical work by Boss (2006) suggested that a close stellar companion could rapidly induce gravitational perturbations and trigger the instabilities needed for gravitational fragmentation to proceed, even if the disc is not initially unstable to its own gravity. However, simulations by Forgan & Rice (2009) indicate that, rather than promoting fragmentation, perturbations from an outer companion are more likely to make the disc more stable.

The brown dwarf desert is thought to be a natural feature arising from formation around single stars, where massive objects with brown dwarf masses can only form at wide separations and can be challenging to bring inwards through disc migration alone. By modifying the circumstellar disc environment, allowing for different conditions facilitating in-situ formation, and/or by triggering migration processes, the presence of a binary companion could help populate the low-mass end of the brown dwarf desert and explain the puzzling existence of the scarce population of very close-in brown dwarf desert inhabitants.

It is worth noting that over half of detected binary companions in our study have projected separations larger than 200 AU. We therefore argue that wider binaries must also be able to impact, almost to the same degree, the formation and/or evolution of these systems. The processes described above must therefore also be possible from wider separations in order to account for the existence of the planets and brown dwarfs probed in this work. An easy way to facilitate this is to form the inner companions at significantly larger orbital distances (tens of AU), increasing the initial gravitational influence of the outer companion. As mentioned previously, this theory could tentatively explain the shortfall of binaries with separations $<$50$-$100 AU, which would then not be expected around such systems.

8.4 Scattering and migration via the Kozai-Lidov mechanism

One way binary companions could assist the migration of these systems is through the Kozai-Lidov mechanism (Kozai, 1962; Lidov, 1962). In this alternative scenario to produce hot Jupiters, an outer binary companion on an inclined orbit relative to the orbital plane of the planet triggers periodic oscillations of the planet’s eccentricity and inclination (Fabrycky & Tremaine, 2007; Wu et al., 2007; Dong et al., 2014; Petrovich, 2015). Combined with the effects of tidal dissipation, these secular interactions can result in a very short orbit for the inner companion (Rice et al., 2015), possibly with a high spin-orbit misalignment with its host star. The amplitude of these interactions mostly depends on the initial mutual inclination between the inner and outer companions (Fabrycky & Tremaine, 2007), allowing Kozai-Lidov cycles to be induced by very distant perturbers.

As noted in Section 2, the subset of targets with a planet or brown dwarf within a 10 day orbit corresponds to the systems in our sample with a tidal circularisation timescale shorter than $\sim$15 Gyr, which we consider to be fully consistent with the Kozai-Lidov scenario. In Section 7.2, we showed that this subsample may have a marginally higher binary fraction than the already very high multiplicity rate observed for our full sample. This is thus compatible with the idea that the inner companions from this subsets could have been driven to their current orbital configurations through Kozai-Lidov oscillations under the effect of wide companions. We also note that 4 of our 5 targets also studied in the Friends of hot Jupiters campaign have high obliquities (placed in the “misaligned” sample in Knutson et al., 2014), a feature often associated with the Kozai-Lidov mechanism.

Unfortunately, full orbital parametrisation, including inclination measurements, is not possible for wide, directly-imaged binaries. Nevertheless, we can determine if the observed wide companions could be responsible for a Kozai-Lidov scattering of the inner planets and brown dwarfs based on their masses and orbital distances. This is done by estimating the minimum companion mass required to excite Kozai-Lidov oscillations on a timescale shorter than general pericenter precession, as was done in Ngo et al. (2016). We adopted a primary stellar mass of 1 M⊙ and a mass of 15 MJup for the inner companion, close to the median of our Kozai-consistent sample. Equating equations (1) and (23) from Fabrycky & Tremaine (2007) for the Kozai-Lidov oscillation timescale and pericenter precession due to general relativity, respectively, we computed in Figure 19 the minimum masses and separations necessary to migrate a 15 MJup companion with initial orbital separations of 1 AU, 5 AU, 10 AU and 20 AU through this scenario. We assumed initially circular orbits for the inner companions and eccentricities of 0.5 for the outer companions, based on the roughly uniform eccentricity distribution between 0 and 1 found by Raghavan et al. (2010) for wide field binaries. We found that almost all detected binary companions could explain the presence of the inner companions in this subset via the Kozai-Lidov mechanism, assuming they formed at separations larger than at least 1$-$10 AU.

These hierarchical systems are hence compatible with a migration of the inner companions through the Kozai-Lidov scenario based on this simple analysis. However, the subset of objects inconsistent the Kozai-Lidov mechanism, based on our tidal circularisation timescale argument (see Section 2), was also found to have a particularly large binary frequency. This suggests that these systems do not primarily migrate via Kozai-Lidov oscillations. This is in good agreement with the theoretical study by Naoz et al. (2012) and observational constraints placed by Ngo et al. (2016), which concurred that only 20 to 30% of all hot Jupiters can be explained by the Kozai-Lidov migration process. Kozai-Lidov oscillations are therefore unlikely to be the dominant mechanism driving close-in massive planets to their current locations.

9 Summary and Conclusions

We have gathered a compilation of 38 planets or brown dwarfs with masses of at least 7 MJup and orbiting within $\sim$1 AU from their host stars, with the aim of examining the multiplicity statistics of these systems. We searched for wide binary companions to these objects using new direct imaging data, observations reported in the literature and the Gaia DR2 catalogue. A total of 16 confirmed multiple systems were found, and another 3 targets have at least one high-probability candidate companion. We report here the discovery in Gaia DR2 of a new mid-K tertiary component comoving with the WASP-14 binary system, and present an independent confirmation of the wide M7.5 companion to WASP-18. We used a robust MCMC statistical approach to constrain the binary properties of our sample, correcting for observational biases and incompleteness. Our main results are summarised below.

1. A very high binary fraction. Our analysis revealed a very large binary frequency of $f=79.0^{+13.2}_{-14.7}$% for these outer companions on separations between 20 and 10,000 AU, which is more than twice as high as for field stars on the same separation range, with a 3-$\sigma$ significance. These results demonstrate that wide binary companions greatly influence the formation or evolution of these close-in massive planets and brown dwarfs. The presence of a binary companion could allow for different natal environments in circumstellar discs, enabling in-situ formation at locations where giant planet formation is not normally possible. Stellar companions could also facilitate disc migration beyond the extent normally seen around single stars, or could trigger alternative migration processes through induced secular interactions.

2. A deficit of close binaries. The output of our statistical analysis also showed a strong preference for wide binaries, with a peak around 250 AU, compared to $\sim$50 AU for the overall field population. The apparent shortfall of $<$50$-$100 AU binaries is consistent with previous studies. It is not clear whether this deficiency indicates that planet formation is inhibited in tight binaries, that our probed planets formed near these separations, or if it is the result of selection biases in exoplanet surveys. Based on these observations, we argue that the mechanisms assisting planet formation or evolution in multiple star systems must be associated with widely-separated binaries, on distances larger than several hundreds of AU. However, we did find that the Kozai-Lidov mechanism is unlikely to be the dominant underlying process.

3. A higher binary rate for higher-mass planets. A comparison with the population of lower-mass planets suggests that binary occurrence increases with planet mass for these close-in objects. This is in good agreement with prior studies that found the most massive planets to be almost exclusively observed in binary systems, and indicates that the role played by binary companions in the existence of these systems becomes more critical for higher-mass planets.

4. A higher binary rate for shorter-period planets. Dividing our sample into two subsets, we found that stars hosting planets or brown dwarfs with orbital periods $<$10 days have a marginally larger binary rate than systems with longer-period inner companions, consistent with a 100% multiplicity fraction at the 1-$\sigma$ level. If confirmed, this trend could suggest that the influence of binarity on the formation/evolution of the most massive planetary companions is enhanced for shorter-period planets, and may even become a requirement for the very closest planets and brown dwarfs.

We conclude that wide binary companions have a crucial influence on planet formation and/or evolution and may be responsible for the sporadic population of high-mass planets and brown dwarf desert members observed on very tight orbital configurations, which seem to rarely exist around isolated stars.

Acknowledgements

We thank our anonymous referee for their insightful comments and suggestions. This paper comes from work undertaken during a visit funded by the Scottish Universties Physics Alliance (SUPA) Postgraduate, Postdoctoral and Early Career Researcher Short-Term Visits funding. This work benefited from the Exoplanet Summer Program in the Other Worlds Laboratory (OWL) at the University of California, Santa Cruz, a program funded by the Heising-Simons Foundation. KM acknowledges funding by the Science and Technology Foundation of Portugal (FCT), grants No. IF/00194/2015 and PTDC/FIS-AST/28731/2017. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). This study is based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), Ministério da Ciência, Tecnologia e Inovação (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). This work is based on observations collected at the European Southern Observatory under ESO programme 099.C-0728. Some of the observations in the paper made use of the NN-EXPLORE Exoplanet and Stellar Speckle Imager (NESSI). NESSI was funded by the NASA Exoplanet Exploration Program and the NASA Ames Research Center. NESSI was built at the Ames Research Center by Steve B. Howell, Nic Scott, Elliott P. Horch, and Emmett Quigley. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This research has made use of the Exoplanet Orbit Database, and the Exoplanet Data Explorer at exoplanets.org, as well as the Extrasolar Planets Encyclopaedia at exoplanet.eu. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the Uni- versity of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

Appendix A Notes on individual targets

A.1 Bound systems

11 Com (HD 107383, HIP 60202) is a common proper motion binary found in the Catalog of Components of Double and Multiple Stars (CCDM; Dommanget & Nys, 2000). The system has a magnitude difference $\Delta V=8.0$ and an angular separation of 9$\aas@@fstack{\prime\prime}$1, corresponding to a projected separation of 850 AU at the distance of 11 Com. From the reported magnitude difference, we infer a mass of 0.7 M⊙ for the secondary using the BT-Settl models (Allard et al., 2012) and stellar parameters given in Table LABEL:t:stellar_properties for the primary.

30 Ari B (HD 16232, HIP 12184) is part of a hierarchical system. Along with 30 Ari A (HD 16246, HIP 12189), it forms a physical pair with a projected separation of 38$\aas@@fstack{\prime\prime}$2 or 1520 AU (Shatsky, 2001). Both components of the F5V+F8V 30 Ari system are in turn close binaries. In addition to the 9.88 MJup planet orbiting 30 Ari B with a period of $335.1\pm 2.5$ days (Guenther et al., 2009), Riddle et al. (2015) found that 30 Ari B is also orbited by another companion, 30 Ari C, with a separation of 22 AU (0$\aas@@fstack{\prime\prime}$536). Roberts et al. (2015) subsequently demonstrated that the B-C pair is indeed comoving and inferred a mass of 0.5 M⊙ for the C component, classified as an M1V dwarf (see also Kane et al., 2015). Moreover, the primary component 30 Ari A is itself a spectroscopic binary with a 1.1 day period (Morbey & Brosterhus, 1974) and a total mass of 1.32 M⊙ (Guenther et al., 2009).

$\tau$** Gem** (HD 54719, HIP 34693) is reported in the CCDM and Washington Double Star (WDS; Mason et al., 2001) catalogues to have a candidate companion at a separation of 1$\aas@@fstack{\prime\prime}$9 and a magnitude $V=11$ mag. The system was determined to be most likely bound in (Mitchell et al., 2013), who estimated the companion to be a K0 dwarf with a mass of 0.8 M⊙ separated by 187 AU, if real. Roberts & Mason (2018) recently provided astrometry for this candidate using data obtained in 2004. They found a separation of 1$\aas@@fstack{\prime\prime}$76 at a position angle of 162.5 deg. This source is also found in Gaia DR2, although it only has a 2-parameter astrometric solution (position only) and therefore does not have parallax and proper motion measurements. From the relative positions of the primary target and candidate in Gaia we are able to confirm the bound nature of this system based on the 10 year baseline between the 2004 observations and Gaia measurements. Figure 20 shows the relative positions of of $\tau$ Gem and its companion at the two epochs, clearly demonstrating that the two objects share common proper motion and thus confirming that they form a physical pair. The CCDM also reports another candidate to $\tau$ Gem at 59″. However this latter source is found in the Naval Observatory Merged Astrometric Dataset (Nomad-I; Zacharias et al., 2004) to have a proper motion inconsistent with that of the primary (see Roell et al., 2012) and we therefore discard this candidate in our survey.

$\upsilon$** And** (HD 9826, HIP 7513) was found by Lowrance et al. (2002) to be a wide common proper motion pair on a 55″separation (750 AU). The secondary stellar component $\upsilon$ And B has a $J$-band magnitude of $9.39\pm 0.03$ mag and was estimated by Lowrance et al. (2002) to have an M4.5 V spectral type and a mass of 0.2 M⊙. The primary is the host to 4 close-in planets and substellar companions (Butler et al., 1997; Butler et al., 1999; McArthur et al., 2010; Curiel et al., 2011). This binary system is also mentioned as a physical pair in Raghavan et al. (2006); Raghavan et al. (2010). Patience et al. (2002) and Roberts et al. (2011) also observed the target but did not have a sufficiently wide field of view to detect the distant companion. Neither studies report any additional, more closely-separated candidates around $\upsilon$ And A.

AS 205 (V866 Sco, EPIC 205249328) is an extremely young ($\sim$0.5 Gyr) T Tauri star part of a hierarchical system in Upper sco (Reboussin et al., 2015). The K5 dwarf, and brightest component of the system, was found by Ghez et al. (1993) to form a common proper motion system separated by 1$\aas@@fstack{\prime\prime}$3 (corresponding to 166 AU at the distance of the system) with a low-mass spectroscopic binary (K7+M0; Eisner et al., 2005). Prato et al. (2003) estimated a mass of a mass ratio of *q$\sim$*0.2 between the A and BC components, suggesting a mass of $\sim$0.22 M⊙ for the binary secondary.

HAT-P-20 has a red M-dwarf companion at a separation of 6$\aas@@fstack{\prime\prime}$86 (500 AU) fainter by $\sim$2 mag (WDS catalogue), which was confirmed by Bakos et al. (2011) to form a common proper motion pair using Palomar sky survey archival data. The companion was successfully imaged in the Lucky Imaging survey by Wöllert & Brandner (2015) but was missed in observations from Ngo et al. (2015) due to the restricted field of view of their data. From the reported photometry and adopting the stellar parameters in Table LABEL:t:stellar_properties, we derive a mass of 0.57 M⊙ for this companion using the BT-Settl models at the age of the system.

HD 41004 B (HIP 28393) was identified in Santos et al. (2002) and Zucker et al. (2003) as the lowest-mass component of a K1V+M2V visual binary with a 0$\aas@@fstack{\prime\prime}$54 separation, corresponding to 23 AU. The system has a $V$-magnitude difference of 3.7 mag and is catalogued as a physical pair in the WDS, CCDM and Tycho-Hipparcos catalogues (see Roell et al., 2012). Both components are hosts to close-in substellar companions: the 0.7 M⊙ primary, HD 41004 A, is orbited by a giant planet at 1.33 AU with a projected mass of $2.54\pm 0.74$ MJup (Zucker et al., 2004), while HD 41004 B (0.4 M⊙) hosts a brown dwarf companion at 0.017 AU with a minimum mass of $18.37\pm 0.22$ MJup (Zucker et al., 2003).

HD 87646 A (HIP 49522) is flagged as a binary in the Tycho and Hipparcos catalogues with a magnitude difference in the Hipparcos $V$-band of $2.66\pm 0.97$ mag. Ma et al. (2016) acquired high-resolution images of the system and found a separation 0$\aas@@fstack{\prime\prime}$26 ($\sim$20 AU) between the G-dwarf primary and K-dwarf secondary. The authors derive a mass of $1.12\pm 0.09$ M⊙ for the A component and estimate a mass ratio of $q\sim 0.5$ for the system. In addition to the 12.4 MJup giant planet found at 0.117 AU around HD 87646 A, Ma et al. (2016) also report a 57 MJup brown dwarf candidate companion on an eccentric 1.6 AU orbit around the primary star.

HD 89744 (HIP 50786) is a wide binary on a 63″separation first identified spectroscopically by Wilson et al. (2001) and subsequently confirmed astrometrically to form a common proper motion pair by Mugrauer et al. (2004). The large angular separation of the binary corresponds to a projected separation of $\sim$2460 AU. Mugrauer et al. (2004) estimated a mass of 0.08 M⊙ for the secondary, near the hydrogen-burning limit. Raghavan et al. (2006) estimated an L0V spectral type for the companion. HD 89744 was also observed by Roberts et al. (2011) with Adaptive Optics on the AEOS telescope, who found a faint candidate companion at 5$\aas@@fstack{\prime\prime}$62 with a magnitude difference of $\Delta I=13\pm 2$ mag that is yet to be confirmed. Other sets of observations with PUEO-KIR at CFHT by Chauvin et al. (2006) or the UFTI data obtained by Mugrauer et al. (2004) do not go deep enough at that separation to retrieve this candidate. Given the observed magnitude difference, we infer a mass of 0.08 M⊙ for this candidate from the BT-Settl isochrones.

HD 114762 (HIP 64426) is a WDS 3$\aas@@fstack{\prime\prime}$2 (140 AU) binary pair confirmed to be comoving by Patience et al. (2002) using data from Keck/NIRC and Shane/IRCAL. Bowler et al. (2009) further characterised the system, estimating an M9 spectral type and inferring a mass of 0.09 M⊙ for the companion. The companion is also reported in the Adaptive Optics survey by Roberts et al. (2011).

HD 156846 (HIP 84856) is reported as a wide, bound binary in the WDS catalogue, with a separation of 5$\aas@@fstack{\prime\prime}$1 (250 AU). Tamuz et al. (2008) characterised the companion to the G0 planet host as an M4 dwarf of mass 0.59 M⊙.

HD 178911 B (HIP 94075) is the fainter component of a 16$\aas@@fstack{\prime\prime}$1 (790 AU) physical pair found in the Tycho-Hipparcos catalogue. The primary component HD 178911 AC is itself a 4.9 AU spectroscopic binary discovered by McAlister et al. (1987). The triple system was established to be comoving by Tokovinin et al. (2000) and subsequently confirmed by Raghavan et al. (2006). Tokovinin et al. (2000) estimated a combined mass of 1.9 M⊙ for the AC component, consistent with the value reported in Mugrauer et al. (2007), while the planet host HD 178911 B has a mass of 1 M⊙ (Mugrauer et al., 2007; Bonfanti et al., 2016).

Kepler-13 A (KOI-13) has been extensively targeted with direct imaging (Adams et al., 2012; Law et al., 2014; Shporer et al., 2014; Wöllert & Brandner, 2015; Kraus et al., 2016). Szabó et al. (2011) reported and confirmed Kepler-13 as a common proper motion system composed of two massive A stars, also found in the CCDM catalogue. Santerne et al. (2012) found the secondary component to be a spectroscopic binary. Johnson et al. (2014) later constrained the mass of Kepler-13 C to be between $0.4-0.75$ M⊙, for a total mass of $1.68\pm 0.10$ M⊙ for the BC component, and $1.72\pm 0.10$ M⊙ for Kepler-13 A, respectively (Shporer et al., 2014). The A-BC system has a projected angular separation of 1$\aas@@fstack{\prime\prime}$15, corresponding to a physical projected separation of 610 AU (Szabó et al., 2011; Adams et al., 2012; Law et al., 2014).

NLTT 41135 was identified by Irwin et al. (2010) as a physically associated companion to NLTT 41136 at 2$\aas@@fstack{\prime\prime}$4 separation (55 AU). From their characterisation of the system, the authors inferred masses of 0.16 M⊙ for NLTT 41135 and 0.21 M⊙ for NLTT 41136, respectively.

WASP-14 was found in Wöllert et al. (2015) to have a candidate companion at 1$\aas@@fstack{\prime\prime}$4, 5.4 magnitudes fainter than the primary in AstraLux Norte observations at the Calar Alto 2.2 m telescope. Ngo et al. (2015) independently identified the same candidate and were able to confirm the source to be a common proper motion companion to WASP-14 with a mass of $0.33\pm 0.04$ M⊙. We also found in this work a distant companion to the system at 1900 AU, identified in the Gaia DR2 catalogue. We characterise WASP-14 C as an 0.28 M⊙ K5 dwarf (see Section 4.2.1).

A.2 Unconfirmed Candidate Companions

70 Vir (HD 117176, HIP 65721) was observed by Roberts et al. (2011) using the Advance Electro-Optical System (AEOS) telescope. The authors report a candidate companion at a separation of 2$\aas@@fstack{\prime\prime}$86 (52 AU) around 70 Vir, which they classify as an M5 dwarf or later. With a magnitude difference of $\Delta I=11.4\pm 1.2$, we estimate a mass of 0.08 M⊙ for this candidate using the BT-Settl models. Pinfield et al. (2006) reported an L-dwarf candidate at 43″(848 AU) based on data from the 2MASS All Sky Catalogue. Common proper motion with the primary has yet to be determined for both candidates. We did not find the latter candidate as a Gaia DR2 source, most likely too faint for Gaia. 70 Vir had also previously been observed by Patience et al. (2002). Observations from this survey are not deep enough to retrieve the candidate found by Roberts et al. (2011) and do not have a large enough field of view to detect the wide source from Patience et al. (2002). Given the faint infrared 2MASS magnitude of the wide source ($J=15.84\pm 0.16$ mag), we estimate a mass of 0.07 M⊙ for the candidate adopting the age of the primary and the BT-Settl isochrones.

EPIC 219388192 is solar twin in the old Ruprecht 147 star cluster (Curtis et al., 2013; Nowak et al., 2017) which was found by Nowak et al. (2017) to host an eccentric transiting brown dwarf companion. The team acquired Subaru/IRCS+AO188 images of the target to search for nearby companions and found two point sources at 6″and 7$\aas@@fstack{\prime\prime}$5 with contrasts of $\Delta H=7.1$ mag $\Delta H=7.7$ mag, respectively. Nowak et al. (2017) estimated that the candidates, if found to be bound, would be late-type M dwarfs with masses less than $\sim$0.1 M⊙. Both sources are found in the Gaia DR2 catalogue with separations and position angles consistent with those reported by Nowak et al. (2017). However, given the relatively small proper motion of the target and the short time baseline between Gaia DR2 and the direct imaging data ($\sim$1 year), we are not able to confirm or refute either of those candidates. Curtis et al. independently studied the same object and found 4 km s*-1* offset between the center-of-mass radial velocity of the star and Ruprecht 147’s bulk velocity (announced at the Cool Stars 19 workshop 444https://doi.org/10.5281/zenodo.58758). As the star’s proper motion supports its cluster membership, Curtis et al. also obtained aperture-masking interferometry with Keck II and uncovered a 0.52 M⊙ companion at 82 mas (24 AU) with a magnitude contrast $\Delta K=2.24$, explaining the observed offset (Curtis et al. private communication).

KELT-1 was found by Siverd et al. (2012) to have a faint candidate companion at $588\pm 1$ mas ($154\pm 8$ AU) based on Keck/NIRC2 AO data. The relative brightness of the candidate was found to be $\Delta H=5.90\pm 0.10$ maf and $\Delta K^{\prime}=5.59\pm 0.12$ mag. The reported photometry suggests a mass of 0.2 M⊙ based on the BT-Settl models. Siverd et al. (2012) estimated an M4$-$5 spectral type and concluded that the companion is physically associated to the primary, with a $\sim$0.05% probability of being an unrelated background star based on Galatic models, in excellent agreement with our estimates (see Table LABEL:t:candidate_binaries). This target was more recently observed by Coker et al. (2018) with the WIYN 3.5 m telescope and by Wöllert & Brandner (2015) with the Calar Alto 2.2 m telescope, although neither of these sets of observations were deep enough the retrieve the candidate identified in Siverd et al. (2012).

A.3 Rejected Candidates

HD 162020 (HIP 87330) had previously been observed with NACO by Eggenberger et al. (2007), who found two point sources within 5″from the star. The first, closer candidate was found by Eggenberger et al. (2007) to be background, while the nature of the second source was inconclusive. With new NACO data for this target, we were able to refute the bound nature of this companion based on the Gaia DR2 astrometry of the primary and a decade-long baseline between the archival and new observations (see Figure 4). Our proper motion analysis of HD 162020 and this companion is presented in Section 3.3.

HD 168443 (HIP 89844) was observed with SPHERE at VLT in the survey conducted in Moutou et al. (2017). Three point sources are reported within 2$\aas@@fstack{\prime\prime}$5 of the primary in that paper. Moutou et al. (2017) found that given the galactic latitude of the target and the crowded field of view at wider separations around this object, the three identified sources are likely background contaminants due to the local environment of HD 168443. Using the Trilegal galactic models (Vanhollebeke et al., 2009) and following the approach described in Section 4.1, we infer probabilities $<$ 1% for any of these three sources to be physical associated to the primary and do not consider them as bonafide candidates for the purpose of our study.

XO-3 has a faint candidate companion ($i$ = 18.43 mag) first reported in Bergfors et al. (2013). The widely-separated candidate (6″or 1500 AU projected separation) was found by Bergfors et al. (2013) to likely be a physically unrelated background object if it is a main-sequence star based on a colour analysis, although the authors mention the possibility of a coeval white dwarf. This target was observed with Keck/NIRC2 in Ngo et al. (2015) but was not retrieved in the field of view of the images acquired for that survey. Wöllert & Brandner (2015) also imaged XO-3 in a search for wide companions and detected the same source in AstraLux Norte data. A faint source is found in Gaia DR2 ($G$ = 18.45 mag) at the same angular separation and position angle as the detected candidate. Given the comparable photometry and the short timespan between the Gaia DR2 epoch (2015.5) and the observations from Wöllert & Brandner (2015) (March 2015), we conclude that this is indeed the same source. The Gaia DR2 source (GDR2 470650457698311296) has a full 5-parameter astrometic solution and has parallax and proper motion measurements highly inconsistent with those of XO-3 in Gaia DR2. We thus conclude that it is an unrelated background object and rule out this candidate.

A.4 Null-detections

BD+24 4697 (HIP 113698) was observed with Gemini North/NIRI as part of this survey. Our data did not reveal the presence of any candidate companion within the field of view and detection limits of our observations.

CI Tau (EPIC 247584113) is a $\sim$2 Myr T Tauri star located in the Taurus star-forming region with an infrared excess in its SED and a dics resolved by Andrews & Williams (2007). It was observed by Uyama et al. (2017) with the Subaru Telescope, using the NIR camera HiCIAO together with the AO188 adaptive optics system, in quad PDI and standard ADI modes. The authors did not find any candidate companion within the 20″$\times$20″field of view of their observations. This targets is also reported to be single in Kraus et al. (2012) based on analyses of 2MASS images and in the HST young binary survey by White & Ghez (2001).

HAT-P-2 (HD 147506, HIP 80076) was found by Lewis et al. (2013) to have a long-term radial velocity trend, suggesting the presence of an outer companion in addition to its known 9 $M_{\mathrm{Jup}}$ planet on an eccentric 5 days orbit. Bonomo et al. (2017) subsequently placed lower limits on the period and mass of this possible outer companion of $\geqslant$49.2 yrs ($\sim$13 AU) and $\geqslant$39.5 MJup based on radial velocity data. This is consistent with results from Knutson et al. (2014) who constrained the companion properties to $M_{2}\sin{i}=8-200$ MJup and $a=4-31$ AU. Observations with NIRC2 on Keck II (Lewis et al., 2013; Ngo et al., 2015) and with AstraLux Norte (Bergfors et al., 2013) did not reveal any companion but only excluded the presence of an $\sim$equal-mass binary from $\sim$10 AU and companions near the hydrogen-burning limit from $\sim$50$-$100 AU. A companion responsible for the observed RV trend could therefore still remain undetected.

HD 5891 (HIP 4715) was observed by Ginski et al. (2016) with the Lucky Imaging instrument AstraLux at the Calar Alto 2.2 m telescope and did not find any companion, achieving contrasts of 4 mag at 1″and 9.5 mag at 5″.

HD 33564 (HIP 25220) is listed in the CCDM as a 25″binary although the 2 components display inconsistent proper motions and do not form a physical pair (Roell et al., 2012). Ginski et al. (2012) acquired observations of HD 33564 and excluded the presence of companions down to the substellar limit on separations of 20$-$100 AU. This star is also reported as a single object in Eggleton & Tokovinin (2008).

HD 77065 (HIP 44259) is one of the two targets we observed with NIRI on Gemini North. We did not find any candidates around this target in our images, ruling out companions at the hydrogen-burning limit from separations of 5 AU and substellar companions with masses $>40$ MJup from 70 AU.

HD 104985 (HIP 58952) was observed with the lucky imaging camera AstraLux on the Calar Alto 2.2-m telescope by Ginski et al. (2012). The team did not find any candidate around this target.

HD 134113 (HIP 74033) is part of the Arcturus moving group. We observed this target with the WIYN telescope and did not find any companions within our detection limits. HD 134113 has no previous direct imaging observations reported in the literature.

HD 156279 (HIP 84171) was observed by Ginski et al. (2016) with the AstraLux instrument on the Calar Alto 2.2-m telescope. No companion was detected in the obtained lucky imaging data.

HD 160508 (HIP 86394) was observed as part in this work using the WIYN imaging facilities. We did not detect any companions around this object within the field of view of our images.

HD 180314 (HIP 94576) was targeted by Ginski et al. (2016) with lucky imaging at Calar Alto. No source was uncovered in the obtained data within 12″, down to low-mass stellar companions.

HD 203949 (HIP 105854) was observed with VLT/SPHERE in Moutou et al. (2017). That survey does not report the detection of any candidates around this target.

WASP-18 (HD 10069, HIP 7562) was part of our observed sample and no source was detected in the field of view of our images. This object had already been observed with Keck II/NIRC2 in Ngo et al. (2015). No candidate was reported around WASP-18 in that survey. We achieved a better contrast than that reported in Ngo et al. (2015) both at diffraction and background-limited separations and our observations allowed us to rule out the presence of lower-mass companions around WASP-18. A comoving object was however found in this work in GDR2 at 3300 AU, outside the field of view of the direct imaging data (26$\aas@@fstack{\prime\prime}$7), for which we estimated a spectral type of M7.5 and a mass of 0.092 M⊙ (see Section 4.2.2).

Appendix B Gaia DR2 Analysis

In Section 4.2 we searched for sources in the Gaia DR2 catalogue with fractional differences of less than 20% in parallax and at least one proper motion component relative to the Gaia astrometry of our targets. Using these selection constraints, we identified a total of 11 binaries in Gaia DR2 among the targets in our sample, 9 of which were previously known. We now examine those systems more carefully as well as the remaining systems from Table LABEL:t:confirmed_binaries in order to evaluate and refine our selection criteria, if needed.

B.1 Binary Completeness

For completeness, we first searched for other known binaries in our sample that may have been missed by our chosen constraints. A total of 7 known, comoving systems are missing from our identified Gaia binaries, corresponding to the companions with no parallax or proper motion listed in Table LABEL:t:confirmed_binaries. From those, 30 Ari BC, HD 41004 AB and HD 87646 AB have angular separations $<$1″, the resolving limit of Gaia DR2, and were therefore missed because of angular resolution limitations.

While near-equal brightness binaries ($\Delta G<1$ mag) are typically resolved with Gaia from separations of $\sim$1″(e.g. Kepler-13, 1$\aas@@fstack{\prime\prime}$15 separation, $\Delta$mag = 0.2 mag; AS 205 AB, 1$\aas@@fstack{\prime\prime}$3, $\Delta G=1$ mag), larger separations are required to resolve lower mass ratio binaries. Ziegler et al. (2018) estimated that companions with $\Delta G$ down to $\sim$6 mag are consistently recovered at separations of 3″$-$5″, with a roughly linear decrease in the recoverable magnitude difference between 1″$-$3″. Based on these results, it is not surprising that systems such as WASP-14 AB (1$\aas@@fstack{\prime\prime}$45, $\Delta$*J$=$*5.2 mag) and HD 114762 AB (3$\aas@@fstack{\prime\prime}$2, $\Delta$*J$=$*7.6 mag) are not retrieved in Gaia DR2. We thus conclude that these companions are missing from our Gaia binary sample because they are fainter than the completeness level of Gaia DR2.

Finally, the last missing binaries are $\tau$ Gem AB and HAT-P-20 AB. In both cases, the two binary components were found to be Gaia DR2 sources, but the fainter component only had a two-parameter astrometric solution (position only) rather than the full 5-parameter solution (position, parallax and proper motion). With no parallax and proper motion measurements, we were not able to select these systems in our analysis and we attribute the fact that we missed them to the remaining incompleteness of Gaia DR2 and not to our selection criteria. We thus conclude that our selection method successfully identified all known binaries that were recoverable.

B.2 Binaries with excessive astrometric disparities

Table 2 reports the relative differences in parallax and proper motion, together with their associated uncertainties, obtained for all identified Gaia binaries (see also Figure 7). While the majority of the errors are comparable in size to the calculated values themselves, all systems remain fully within our arbitrary cuts at 20% at the 1-$\sigma$ level (with the exception of the newly discovered WASP-18 AB system which is discussed in Section 4.2.2).

A number of binaries in Table LABEL:t:confirmed_binaries are part of hierarchical systems and we find that 4 of the 9 previously known Gaia systems have an unresolved component in Gaia DR2 (30 Ari BC, AS 205 BC, HD 178911 AC and Kepler-13 BC), which correspond to the blue stars in Figure 7. Looking at the positions of these specific systems in the parameter-space from Figure 7, we find that they have the largest relative offsets in parallax and/or proper motion, and are the only systems for which the relative difference in proper motion was larger than our 20% threshold in one of the coordinates (outside the shaded area in Figure 7).

This is consistent with the idea that unresolved components can have a significant effect on the measured astrometry of binary pairs, reinforcing the argument for loose constraints in order to ensure that such hierarchical systems are not missed. With the exception of AS 205 and HD 178911, all known binaries detectable by Gaia would also have made a more stringent cut at $\sim$10% in the relative difference in parallax and in one of the proper motion components. Furthermore, the 5 known binaries that are not known to have an unresolved component (blue circles in Figure 7) also make that 10% cut in both proper motion components. We thus conclude that most binaries should have relative discrepancies of $<$10% in all astrometric parameters ($\pi$, $\mu_{\alpha*}$ and $\mu_{\delta}$), while systems agreeing to within 20% in parallax and in one of the proper motion coordinates are likely to be hierarchical systems with an unresolved component.

We note that such wide companions are not necessarily presently bound systems. Formerly physically associated components of a binary system may continue to travel along a nearly identical trajectory. However, we are seeking in this study companions that may have affected the formation or early evolution of inner companions and therefore also consider as bonafide any pair that previously constituted a bound system. We also point out that such a configuration would likely result in small discrepancies in the observed astrometric parameters of the individual components, an additional argument for the loose constraints considered above. In conclusion we trust that systems passing the selection criteria described above have consistent astrometric parameters and kinematics, and may be treated as binaries for the purpose of this work.