A Possible Alignment Between the Orbits of Planetary Systems and their Visual Binary Companions

2.1.1. Identification with TESS

2.1.2. Sample of Visual Binaries with Planet Candidates

2.1.3. Follow-up Ground-based Time-series Photometry

2.2.1. Control Sample of Visual Binaries without Planet Candidates

We also identified a control sample of visual binary systems from the El-Badry & Rix (2018) catalog. Kepler has taught us that most of these stars likely host planetary systems of their own (e.g., Fressin et al. 2013; Deacon et al. 2016), but since they do not host any transiting planets, their inclinations will be unknown. Therefore, performing our analysis on a control sample and comparing the results against the sample of binaries with planet candidates helps give us confidence that any features we see in the resulting distribution of inclination angles are astrophysical and not due to selection effects. We specifically constructed our control sample to have nearly identical properties to the sample of binaries with planet candidates to make sure that our control sample incorporates any selection biases from the TESS planet-detection process. To achieve this goal, we defined a metric, ${ \mathcal M }$, to quantify the similarity between any two visual binary systems:

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal M } & = & {\left(\displaystyle \frac{{\rm{\Delta }}{G}_{1}}{4}\right)}^{2}+{\left(\displaystyle \frac{{\rm{\Delta }}{G}_{2}}{4}\right)}^{2}+{\left(\displaystyle \frac{{\rm{\Delta }}{{RP}}_{1}}{4}\right)}^{2}+{\left(\displaystyle \frac{{\rm{\Delta }}{{RP}}_{2}}{4}\right)}^{2}\\ & & +{\left(\displaystyle \frac{{\rm{\Delta }}{{BP}}_{1}}{4}\right)}^{2}+{\left(\displaystyle \frac{{\rm{\Delta }}{{BP}}_{2}}{4}\right)}^{2}+{({\rm{\Delta }}\varpi )}^{2}+{\left({\rm{\Delta }}s\right)}^{2},\end{array}\end{eqnarray} \tag{ 3 }$$

where s is the projected separation of the two stars in the binary system, ϖ the system parallax, and G, BP and RP are the stars' apparent magnitudes in the three Gaia passbands. Here, the Δ symbol represents the normalized fractional difference between the values for the two systems: a system with a transiting exoplanet and potential control sample system, and the subscripts 1 and 2 represent the primary and secondary star in each system. For instance, ${\rm{\Delta }}{G}_{1}=\tfrac{({G}_{1}-{G}_{c})}{{G}_{c}}$, where c represents the control sample. We arbitrarily divide all magnitude normalized differences by 4 so that not all weight is given to the magnitudes.

For each of the visual binary systems with non-false-positive exoplanets as of 2020 December 25, we identified the 12 systems from the El-Badry & Rix (2018) catalog with the lowest ${ \mathcal M }$ metric. Systems were not removed after each sampling procedure (i.e., they are allowed to be included twice); however, due to the large number of systems present in the El-Badry & Rix (2018) catalog, the resulting control sample has no repeated systems. A subset of our planet-candidate sample was also identified by El-Badry & Rix (2019) to have spectroscopic metallicity measurements from one of several large spectroscopic surveys (see Section 2.3.1). For these systems, we restricted our search for similar systems to those that also have an archival spectroscopic metallicity measurement from El-Badry & Rix (2019). In total, we identify a control sample of 960 systems with very similar distributions of parameters to the input sample of binaries containing planet candidates. Figure 3 shows various properties of the sample with exoplanets and control sample.

Zoom In Zoom Out Reset image size

2.2.2. Astrometric Parameters

2.2.3. Removing Incorrect Cross-matches, High-RUWE Solutions, and White Dwarfs

We apply a variety of cuts to both the control sample and sample with exoplanets in order to ensure that only high-quality astrometric parameters are preserved.

In the process of converting between Gaia DR2 and Gaia EDR3 IDs, a purely positional cross-match can contaminate the sample due to proper motion movement from the Gaia DR2 epoch (2015.5) to the Gaia EDR3 epoch (2016) and the addition of new sources in EDR3. To ensure that there are no incorrectly cross-matched stars in our sample, we exclude 17 binary systems in the control sample for which ∣G_EDR3−G_DR2∣ > 0.05.

The renormalized unit weight error (RUWE) can be used as an indicator of the quality of the Gaia astrometric solution for a star (Lindegren 2018). The RUWE is the square root of the reduced χ² divided by a correction function that eliminates dependence on G magnitude and BP − RP color. An RUWE of greater than 1.4 typically indicates a poor astrometric fit, so we eliminate any systems for which the RUWE for at least one of the stars is greater than 1.4. A high RUWE can indicate the presence of an unresolved companion (Belokurov et al. 2020).

We also remove any binaries where either the host star or companion star is a white dwarf; it is more difficult to estimate masses for white dwarfs than for main-sequence stars, and in these systems the binary orbit has been influenced by post-main-sequence mass loss. While these effects are very interesting in their own right, it is beyond the scope of this work to consider them.

After these cuts, there are 67 binary systems with exoplanets and 688 binary systems in the control sample. The distribution of the radii and the periods of the exoplanets in our sample are shown in Figure 4.

Zoom In Zoom Out Reset image size

2.2.4. Other Work Identifying Visual Binaries in Gaia

Fitting binary orbits using only instantaneous positions and proper motions from Gaia requires an estimate of the mass of each binary component, which in turn requires an estimate of each star's metallicity. To derive the metallicities of stars in our samples, we use both archival observations from large spectroscopic surveys and targeted follow-up observations of planet-candidate host stars made by the TESS Follow-up Observing Program (TFOP). Below, we describe the sources of our spectroscopic parameters and the procedure we used to determine the metallicities of the observed stars. Because the components of relatively wide-binary stars are known to have nearly identical elemental abundances in most cases (Hawkins et al. 2020), we assume the metallicity of both stars in the binary are the same when we only have metallicity measurements for one of the pair.

For stars that have more than one spectroscopic observation, we use an average of the metallicities derived from the separate spectroscopic observations and add the errors from each separate observation in quadrature.

2.3.1. Archival Spectroscopy

2.3.2. Las Cumbres Observatory/Network of Robotic Echelle Spectrographs

2.3.3. Tillinghast Reflector Echelle Spectrograph

2.3.4. FIbre-fed Echelle Spectrograph

2.3.5. CTIO High Resolution Spectrometer

We derive masses for the stars in our samples using the Isochrones Python package (Morton 2015). Given observable parameters like a star's apparent magnitudes in different bandpasses, its trigonometric parallax, and spectroscopic parameters, Isochrones determines the star's most likely physical parameters and their uncertainties by comparing measured parameters to those predicted by stellar evolution and atmosphere models. Isochrones supports several different suites of isochrone models; we use the MIST isochrones (Choi et al. 2016; Dotter 2016). We input each star's parallax, metallicity (when available from archival spectroscopy or NRES, TRES, FIES, or CHIRON), and apparent magnitudes in the G, BP, and RP bandpasses from Gaia and the J, H, and K bandpasses from the 2MASS survey (Skrutskie et al. 2006). For uncertainties in bandpasses, we use the error floors of Eastman (2017).

Not all of the stars in our samples have spectroscopic metallicity measurements. If a star does not have a metallicity value, the metallicity of its companion is used (see Hawkins et al. 2020). If neither star in a binary pair has a spectroscopic metallicity measurement, we use the isochrones default metallicity ⁹³ , a prior based on a two-Gaussian fit to the distribution of metallicities reported by Casagrande et al. (2011).

Model isochrones are not as well calibrated for M-dwarf stars (which constitute around half of the stars in the sample of visual binaries) as they are for Sun-like stars (Angus et al. 2019). For M dwarfs, we use the M_K − M_⋆ relationship of Mann et al. (2019) using the accompanying M_-M_K- Python package ⁹⁴ to estimate the stars' mass. M_k magnitudes, when available, are taken from the TIC (Stassun et al. 2019). Any stars in the range 4.5 < M_k < 10.5 (the more conservative option given by Mann et al. 2019) are treated as M dwarfs.

Rough mass estimates for these stars are also provided in the TIC (Stassun et al. 2019). We compare our mass estimates to the TIC estimates, as demonstrated for the control sample in Figure 5, to ensure that there are no systematic discrepancies in estimation. The median uncertainty of the masses in our sample is 0.016 M_⊙ and in the TIC 0.082 M_⊙. We also compare our mass estimations to those derived from spectral energy distribution modeling (Stassun et al. 2018), where stellar atmosphere model grids are interpolated in T_eff, $\mathrm{log}g$, and [Fe/H] and combined with spectroscopic $\mathrm{log}g$ measurements. There is general agreement between the masses derived from both techniques, an independent check of the masses assigned to the stars.

Zoom In Zoom Out Reset image size

Linear Orbits for The Impatient (LOFTI; Pearce et al. 2020) uses measurements of the relative positions of two stars (in both R.A., Δα, and decl., Δδ), relative proper motions, relative radial velocity (if available), and stellar masses to constrain orbits of visual binaries. In our analysis, we take all of these parameters from Gaia. LOFTI calculates relative proper motion and position as primary minus companion, while relative radial velocity is calculated as companion minus primary. LOFTI uses the Orbits for the Impatient (OFTI) sampling and rejection technique (Blunt et al. 2017) applied to Gaia astrometric parameters.

OFTI, on which LOFTI is based, optimizes orbit fitting via rejection sampling with the unique scale and rotate step. The algorithm randomly samples eccentricity (e), argument of periastron (ω), mean anomaly from derived epoch of periastron (t_o ), and cosine of inclination ($\cos (i)$) from uniform priors, then scales the semimajor axis (a) and rotates the longitude of the ascending node (Ω) to match the observed separation in binaries, rejecting orbits if the χ² probability (proportional to ${e}^{-{\chi }^{2}/2}$) is greater than a randomly chosen number in the range [0, 1] to form a posterior distribution of orbital parameters. OFTI is computationally advantageous over traditional Markov Chain Monte Carlo methods for visual binary systems, where orbital parameters are often poorly constrained (Blunt et al. 2017, 2020). Table 2 is a list of priors used in the LOFTI algorithm on sampled parameters, while Table 3 shows parameters calculated either from sampled parameters or via the scale and rotate method.

Table 2. List of Priors for Orbital Parameters that are Randomly Sampled in LOFTI

Parameter	Prior
Eccentricity (e)	Uniform(0, 1)
Inclination (i)	Sin(0, 180)
Argument of periastron (w)	Uniform(0, 360)
${t}_{\mathrm{const}}$	Uniform(0, 1)
Total mass (Mtot)	Normal(M0,Mstd)
Distance	Normal(D0,Dstd)

Note. ${t}_{\mathrm{const}}$ is a value used in the calculation of the epoch of periastron passage (t₀).

Download table as:  ASCIITypeset image

Table 3. List of Orbital Parameters that are Calculated Either from Sampled Parameters or Calculated Using the Scale and Rotate Method to Match Observed Data

Parameter	Formula
Period (T)	$\left(\tfrac{{a}^{3}}{{M}_{\mathrm{tot}}}\right)$
Epoch of periastron passage (t0)	${d}_{\mathrm{epoch}}-{t}_{\mathrm{const}}\ast T$
Semimajor axis a	Scale and rotate
Longirude of ascending node (Ω)	Scale and rotate

Note. d_epoch is the date of the epoch of the observations.

Download table as:  ASCIITypeset image

Orbital fitting using very short observational baselines (as done here) can lead to degeneracies in parameters. Particularly, for high eccentricity values, inclination is biased toward edge-on orbits (Ferrer-Chávez et al. 2021). The inclusion of a control sample in our study should allow us to determine if there is a preferential alignment regardless of the presence of this bias.

We ran LOFTI on the binary systems with exoplanets and those in the control sample. We continued sampling the posterior distributions until we reached 100,000 accepted orbits. We performed this computationally intensive task on the Maverick 2 and Lonestar 5 supercomputer clusters at the Texas Advanced Computing Center.

We show 100 orbits drawn from a sample LOFTI posterior in Figure 6 and a sample posterior distribution of parameters from LOFTI for one binary system in Appendix A. We also show archival ground-based images of two characteristic systems, one aligned and one misaligned, annotated with astrometric information from Gaia in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We assessed the statistical significance of the apparent alignment between visual binary orbits and their planetary systems using a Kolmogorov–Smirnov (K-S) test. A K-S test measures the maximum difference between the cumulative distributions of two empirical distributions and then calculates a p-value representing the probability that the two empirical distributions were drawn from the same underlying distribution. A small K-S value indicates that the distributions being compared are not from the same underlying distribution. We perform a K-S test between the sample with exoplanets and control sample's ∣90 − i∣ values, where i is the median value of the distribution of inclinations. The test rejects the null hypothesis that the two data sets are drawn from the same underlying distribution with p = 0.0037. To check that our significance value is not dependent on our somewhat arbitrary choice of RUWE cutoff, we adjust our RUWE cutoff to 1.3, 1.2 and 1.1, finding that the p-values are, respectively, p = 0.016, p = 0.051, and p = 0.045. Additionally, we use two other statistical tests specifically designed for comparing angular data: Wallraff's nonparametric test (Wallraff 1979) and Watson's nonparametric test (Watson 1983). The Wallraff test yields a p-value of 0.031, while the Watson test yields a p-value in the range [0.001, 0.01] (the test only gives a range of p-values).

None of these tests account for error in the individual measurements of inclination. The K-S test is shown visually as the maximum difference in cumulative distributions in Figure 11.

Zoom In Zoom Out Reset image size

Since most of the excess of aligned systems are found in binary systems with semimajor axes less than about 700 au, we also compare the distributions of ∣90 − i∣ for these close binaries. For systems with semimajor axes less than 700 au compared to the control sample, the K-S test shows stronger evidence that the inclinations in the exoplanet sample and control sample are drawn from different distributions (p = 0.000250). Finally, we investigate whether the distributions of inclination above and below 700 au are significantly different from each other. For systems with semimajor axes less than 700 au compared to those with semimajor axes greater than 700 au, the K-S test still shows evidence for a difference between the distributions, but with lesser significance: p = 0.0172.

The standard deviation in inclination is on average 14°, and a large proportion of samples have a skewed distribution of error. Additionally, error varies with inclination. The distribution of errors in our sample is shown in Figure 12(a), and the dependence of error on inclination in Figure 12(b). To explore the effect of this large and often asymmetrical error on the results of the K-S test, we perform a bootstrap procedure on the two distributions where we take a random sample from each binary system's LOFTI results and compare the overall distributions of those random samples. We repeat this process 10,000 times to derive a bootstrap distribution. This process calculates the distribution of p-values that we would expect to measure if we ran the same experiment many times, with slightly larger uncertainties on the inclinations (by about a factor of $\sqrt{2}$). We found that when we performed this analysis comparing the full samples, 61.91% of the bootstrap instances yielded a p-value less than 0.05, with the median p-value being 0.0300. We run the same bootstrap test on the subset of systems with a < 700 au. This test returns a median p-value of 0.0044, with 87.35% of the bootstrap instances being at least p = 0.05.

Zoom In Zoom Out Reset image size

We also simulated a K-S test of two $\sin (90-i)$ probability density functions with sample sizes similar to those of that in our study and various types of skewed error (e.g., skewed away from i = 90 or toward i = 90). All simulations returned uniform distributions of K-S p-values, indicating that any skewed error is likely not affecting the result of the K-S test.

For completeness, we tested whether the distribution of binary inclinations in the control sample is consistent with randomly oriented orbits. The expected probability density function of binary inclination angles for random orbital orientations is proportional to $\sin (i)$. When we compare the control sample to inclinations drawn from a $\sin (i)$ distribution with a K-S test, we find a significant difference between the two. We show a comparison of the control sample distribution and the expected isotropic distribution in Figures 8(b) and 11. This difference is likely caused by either subtle biases when fitting astrometric orbits in small-orbit arc fitting (Ferrer-Chávez et al. 2021) or selection effects in the El-Badry & Rix (2018) catalog. The control sample has slightly fewer face-on systems and more edge-on systems than we would expect from a perfectly random distribution of inclinations. Nevertheless, the difference between the distribution of inclinations in the control sample and the sample of binaries with exoplanets is still significant.

4.1.1. Inclination–Eccentricity Degeneracy

There is an apparent degeneracy between inclination and eccentricity, presumably related to the bias mentioned in Ferrer-Chávez et al. (2021). We plot inclination versus eccentricity for our systems in Figure 13, where it is clear that LOFTI does not find any systems with low eccentricity and inclination near 90°. This is not a feature of the OFTI method specifically: all orbits derived through short-arc fitting will have a similar degeneracy between median inclination and eccentricity. Such a degeneracy is only between summary statistics of eccentricity and inclination. The full parameter space is covered among all individual samples.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have no reason to suspect that the eccentricity of the binary orbit could somehow be preferentially modifying the orbit of the planet, so do not believe that this degeneracy can explain the overabundance of systems with inclination near 90°. In fact, the difference in the distribution of eccentricity in the control sample versus the sample with exoplanets is not statistically significant (p = 0.119).

We fit a mixture model to the orbital inclinations in our sample of binaries with planets in an attempt to determine what fraction of the binary systems with exoplanets can be explained by a random distribution versus what fraction is aligned. The probability density function of i can be expressed as $\sin (i)$ (although selection effects slightly bias the distribution, as shown in the previous section), while we represent the distribution of the aligned binaries' inclinations as a normal distribution. Our likelihood function is defined as

$$\begin{eqnarray}&&p(y| \lambda ,\sigma )=\prod _{n=1}^{N}(\lambda \cdot { \mathcal N }({y}_{n}| \mu ,\sigma )+(1-\lambda )\cdot \sin {y}_{n}),\end{eqnarray} \tag{ 4 }$$

where y are the binary orbital inclinations, ${ \mathcal N }$ is a normal distribution with mean μ and standard deviation σ, and λ is the mixture parameter or the fraction of systems in the aligned population. We impose the following priors: μ is fixed to be π/2 radians (perfectly aligned with the planet orbits), σ is constrained with a half-normal distribution with standard deviation $\tfrac{\pi }{2}$ radians, and λ, the mixture parameter, has a uniform prior from 0 to 1. Specifically,

$$\begin{eqnarray}&&\sigma \sim {\rm{Normal}}(0,\displaystyle \frac{\pi }{2}),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\sigma \gt 0,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\lambda \sim {\rm{Uniform}}(0,1),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\mu =\pi /2.\end{eqnarray} \tag{ 8 }$$

We run the mixture model in the probabilistic programming language Stan (Riddell et al. 2018). The resultant distribution of λ is peaked near 50%, with the 90% confidence interval being [0.33,0.72]. Though the 90% confidence interval is fairly wide, it does allow us to conclude that, assuming our Gaussian + $\sin (i)$ model for the distribution is a reasonable approximation, the process of alignment occurs in between 26% and 72% of the systems in our sample with 90% confidence.

While we were careful to construct a control sample of binaries with nearly identical properties to the exoplaet sample, there are a few subtle differences that could plausibly introduce a difference in the distribution of inclinations. Here, we show that these differences do not significantly affect the distribution of inclination angles within the control sample, and therefore are highly unlikely to significantly contribute to the apparent alignment between planetary systems and wide-binary orbits.

4.3.1. Metallicity Measurements

Some of the binary systems with exoplanets have metallicity measurements from our own TFOP follow-up spectroscopic observations (see Section 2.3.1). We do not have metallicity measurements for any binaries in the control sample aside from those with archival observations identified by El-Badry & Rix (2019), so it is plausible that including the TFOP spectroscopy could introduce a small difference between the samples. To test this, we compare the inclination distributions of binary systems with spectroscopic metallicity measurements to those without using a K-S test.

First, we tested whether there is a significant difference between the inclination distribution of binaries with and without TFOP spectroscopy in our exoplanet sample. A K-S test reveals no statistically significant difference between these two subsamples (p = 0.32). In fact, if anything, the addition of spectroscopy leads to less of an alignment, as can be seen in Figure 15(a).

Zoom In Zoom Out Reset image size

We also tested whether there is a significant difference between the inclination distribution of binaries with and without archival spectroscopic metallicity measurements in our control sample. Here, the K-S test also reveals no statistically significant difference between those samples in the control sample with and without archival spectroscopy (Section 2.3.1) as can be seen in Figure 15(b) (p = 0.699).

These tests show that, evidently, the presence or absence of a spectroscopic metallicity measurement has a negligible impact on the resulting distribution of binary orbital inclinations. This effect is therefore highly unlikely to explain the apparent alignment of wide-binary orbits with their planetary systems.

4.3.2. 2 Minute Cadence versus FFI Light Curves

As described in Section 2.1.1, TESS observes 20,000 of the best planet-search target stars at 2 minutes cadence (Barclay et al. 2018; Stassun et al. 2019) per sector, while the rest of the sky is observed at 30 minutes cadence (and more recently, at 10 minutes cadence in the extended mission). One of the selection criteria for 2 minute cadence observations is the extent to which nearby stars dilute the signal of the target star. Therefore, stars chosen for 2 minute cadence observations were less likely to be nearby other bright stars, potentially including visual binary companions. This could lead to an increase in edge-on binaries in 2 minute cadence observation since companions could be more likely to be very close to the primary star, and therefore less likely to be resolved in ground-based, seeing-limited imaging and consequently less likely to be rejected due to the companion's flux dilution. Because transit searches in 2 minute cadence observations are more sensitive to detecting planet candidates, this could plausibly introduce a bias in our results.

We tested whether the selection of 2 minute targets significantly affects the inclination distribution of wide-binary companions by comparing the distributions of stars in our control sample that were and were not observed in 2 minute cadence mode. As demonstrated in Figure 16, we find no significant difference in the distribution of inclinations for stars in the control sample observed in 2 minute cadence. A K-S test fails to reject the null hypothesis that binary inclinations for stars observed at 2 minute cadence and those observed in FFIs are drawn from the same population (p = 0.85). The TESS 2 minute cadence target selection is therefore highly unlikely to be a source of bias.

Zoom In Zoom Out Reset image size

4.3.3. Additional Biases

The Lidov–Kozai mechanism induces oscillations in a pair of objects (in this case the planet and star) from the secular influence of a massive, faraway companion. Bodies undergoing the Lidov–Kozai effect experience oscillations in their inclination and eccentricity, trading off between the two quantities on a characteristic timescale set by the system parameters. The Lidov–Kozai mechanism has been invoked to explain a wide variety of astrophysical phenomena, including the existence of hot Jupiters (Fabrycky & Tremaine 2007; Naoz et al. 2012; Dawson & Johnson 2018; Li et al. 2020).

Since the planet mass is much smaller than the inner stellar mass, the Lidov–Kozai effect will only drive substantial dynamical evolution when there is a significant mutual inclination (≳40°) between the orbit of the planet–star system and the orbit of the distant binary companion (Naoz et al. 2013). The fact that the Lidov–Kozai effect operates on systems with large mutual inclinations means that it could plausibly explain the alignment we observe between the orbits of wide-binary stars and their close-in planetary systems. If, for example, the inclinations of wide-binary orbits were randomly distributed with respect to the orbits of the planetary system, the Lidov–Kozai effect would only act to significantly alter the geometries of systems with large initial misalignments. If the Lidov–Kozai effect preferentially disrupts systems with large initial misalignments (either by driving the eccentricity high enough to cause a collision between the planet and the star, or by causing dynamical instabilities that lead to planetary collisions or ejections), the resultant observed distribution of relative inclinations could develop some nonuniformity, as seen in our sample.

To assess whether the Lidov–Kozai effect could explain the alignment between wide-binary orbits and planetary orbits we observe, we calculate the timescale for Lidov–Kozai oscillations to take place for each system, assuming a large enough primordial mutual inclination for the effect to take place. If the Lidov–Kozai timescale is significantly smaller than the total system age, it could have acted to sculpt the observed alignment. The oscillations induced by the Lidov–Kozai timescale occur on a characteristic timescale, τ_LK, which can be estimated as

$$\begin{eqnarray}&&{\tau }_{\mathrm{LK}}\approx \displaystyle \frac{8}{15\pi }\left(1+\displaystyle \frac{{m}_{h}}{{m}_{c}}\right)\left(\displaystyle \frac{{P}_{c}^{2}}{{P}_{p}}\right){(1-{e}_{p}^{2})}^{(3/2)},\end{eqnarray} \tag{ 9 }$$

where m_h is the mass of the planet host star (i.e., the primary star), m_c is the mass of the stellar companion, P_c is the period of the binary companion, P_p is the period of the planet, and e_p is the eccentricity of the planet (Holman et al. 1997; Antognini 2015). The period of the planet is very small compared to the mutual period of the stars. We calculate the Lidov–Kozai timescale for each system in our sample and show the results in Figure 17 as a function of the semimajor axis of the binary system.

Zoom In Zoom Out Reset image size

We find that in systems where the Lidov–Kozai timescale is significantly less than the total age of the system (the systems below the dashed line in Figure 17) the systems are preferentially aligned. This is a consequence of the fact that most systems with binary semimajor axes smaller than about 700 au tend to have orbits aligned with their planetary systems in our sample. As a result, it is plausible that the Lidov–Kozai effect would operate rapidly enough in systems below 700 au to contribute to the alignment we see.

However, we suspect that the Lidov–Kozai effect is likely not a dominant effect causing the observed alignment between binary and planetary orbital planes. One reason is that for planets orbiting very close to their stars, even small apsidal precession rates (due to either general relativity, tides, or the stellar quadrupole moment) can completely suppress the Lidov–Kozai effect (Sterne 1939; Murray & Dermott 1999; Fabrycky & Tremaine 2007). Also, any sculpting by the Lidov–Kozai effect would tend to leave systems with mutual inclinations with up to 40° misalignments, while the overdensity of aligned systems in our sample is strongest within 10° of alignment. Thus, it is likely that if the Lidov–Kozai effect is acting in systems, the effect is quite small and would need to be coupled with other effects to explain the observed alignment.

Another process that could explain the preferential alignment of wide-binary orbits with their close-in planetary systems is that the gravitational influence of the wide-binary star could torque the protoplanetary disk into alignment early in the system's history. Consider, for example, a wide-binary companion in an initially misaligned orbit around a star with a gas-rich protoplanetary disk. The mechanism of Batygin (2012; see also Bate et al. 2000) could be induced by the wide-binary companion, causing the protoplanetary disk to precess in the reference frame of the binary-star system, which would manifest as an oscillation of inclination in the reference frame of the host star's angular momentum vector. As the protoplanetary disk precesses, it can simultaneously dissipate energy and move toward its lowest-energy state, an alignment with the binary-star angular momentum and disk angular momentum vectors (Bate et al. 2000; Martin & Lubow 2017). Any planets that subsequently form from the disk would tend to be in well-aligned orbits with the wide-binary companion.

The warping of a disk can lead to dissipation of energy in the disk. The effect of this warping on the alignment of the disk has been studied in Bate et al. (2000) and Zanazzi & Lai (2018) in a gaseous disk. The timescale of alignment due to warping of the gaseous disk may be compatible with the timescales we observe assuming an outer disk radius of approximately 100 au according to Zanazzi & Lai (2018). However, it is well known that any gas disk, regardless of warping, will thermalize energy due to turbulence or viscosity (Nelson et al. 2000; D'Alessio et al. 2006). As a disk radiates away this energy during precession, it will move to its lowest-energy state, an alignment with the binary-star system.

To explore the possibility that the disk is dissipating energy and moving toward an alignment via a precession, we use the approximate closed-form expression of Batygin's (2012) supplementary materials Equation (6) to estimate the timescale of one precession, and thus one cycle of oscillation in inclination. The timescale of precession is not the same as the timescale of alignment, but given that sufficiently fast precession is a necessary condition for alignment, this timescale can be used to somewhat estimate if alignment could take place in the systems in our sample. Batygin (2012) gives

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{{\rm{\Omega }}}_{\mathrm{disk}}}{{\rm{d}}t}=-\displaystyle \frac{{\int }_{{a}_{\mathrm{in}}}^{{a}_{\mathrm{out}}}G{\rm{\Sigma }}{(r){M}_{c}\left(\tfrac{r}{{a}_{c}}\right)}^{2}{\tilde{b}}_{3/2}^{(1)}(0.95)\mathrm{dr}}{4{\int }_{{a}_{\mathrm{in}}}^{{a}_{\mathrm{out}}}{\rm{\Sigma }}(r)\sqrt{{{GM}}_{\star }{r}^{3}}\mathrm{dr}}\cos (i),\end{eqnarray} \tag{ 10 }$$

where a_c and M_c are the semimajor axis and mass of the companion star, respectively, M_⋆ is the mass of the host star, i is the inclination of the disk with respect to the binary orbit, Σ(r) is the disk density profile, a_out and a_in are the inner and outer boundaries of the disk, and ${\tilde{b}}_{3/2}^{(1)}$ is a Laplace coefficient with disk softening (see Batygin 2012, supplementary materials Equation (4) for details). Since binaries with smaller separations tend to have smaller outer disk radii (Manara et al. 2019), we calculate all disk timescales for a_out in the range $[{a}_{b}^{* }0.05,{a}_{b}^{* }0.25]$ if a_b < 400 and otherwise [20, 100], where a_b is the semimajor axis of the binary-star system. Manara et al. (2019) calculated the dust radii of disks, but the gas radii of disks (which is most responsive to torques) is expected to be larger than the dust radius, although the two tend to be proportional (Ansdell et al. 2018). In general, τ_Ω has a very weak dependence upon the value of a_out as long as a_out is within a reasonable range of values. We use a density profile of the form

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{c}{\left(\displaystyle \frac{r}{{R}_{c}}\right)}^{\gamma }\exp \left[-{\left(\displaystyle \frac{r}{{R}_{c}}\right)}^{2-\gamma }\right],\end{eqnarray} \tag{ 11 }$$

where γ is 0.9 (Andrews et al. 2010) and R_c is 110 au, the mean Andrews et al. (2010) value. Assuming the parameters of Equation (10) do not vary greatly over the precession cycle, the timescale of one complete precession cycle can be estimated as

$$\begin{eqnarray}&&{\tau }_{{\rm{\Omega }}}={\pi }^{2}\left(\displaystyle \frac{4{\int }_{{a}_{\mathrm{in}}}^{{a}_{\mathrm{out}}}{\rm{\Sigma }}(r)\sqrt{{{GM}}_{\star }{r}^{3}}\mathrm{dr}}{{\int }_{{a}_{\mathrm{in}}}^{{a}_{\mathrm{out}}}G{\rm{\Sigma }}(r){M}_{c}{\left(\tfrac{r}{{a}_{c}}\right)}^{2}{\tilde{b}}_{3/2}^{(1)}(0.95)\mathrm{dr}}\right).\end{eqnarray} \tag{ 12 }$$

Here, Equation (12) is averaged over possible inclination values (the mean of ${\cos }^{-1}(i)$ for isotropic inclinations is π/2). We calculated the disk precession timescales for each system in our sample and show the resulting timescales in Figure 18. We find that systems with disk precession timescales less than a few megayears tend to be preferentially aligned, while systems with longer precession timescales show more large misalignments. It is noteworthy that a timescale of a few megayears happens to be the typical lifetime for protoplanetary disks (Mamajek 2009), indicating that this mechanism could explain both the existence of an alignment and the greater level of alignment for binaries closer than about 700 au.

Zoom In Zoom Out Reset image size

One possibility is that the alignment between wide-binary and close-in planetary orbits is primordial. That is, the planets and binary each formed in nearly their current aligned configuration. Binary-star systems at relatively low separations can form through disk fragmentation and turbulent fragmentation. In disk fragmentation, the initial gravitationally unstable circumstellar disk fragments during the collapse of the system, forming a binary or higher-multiplicity system (Adams et al. 1989; Bonnell & Bate 1994; Sigalotti et al. 2018). Disk fragmentation is expected to form binary systems where the stellar angular momenta are aligned with the binary orbit and the plane of the protoplanetary disk. Evidence for disk fragmentation includes aligned bipolar outflows (Tobin et al. 2013) and the metallicity-dependent binary fraction (El-Badry & Rix 2018). Recent observations have shown that disk fragmentation primarily takes place for binary stars with semimajor axes less than about 200 au (Tobin et al. 2016; El-Badry & Rix 2019). More distant binaries likely formed via only turbulent fragmentation, which is not necessarily expected to produce binary stars with orbits in the same plane as any planetary system. Such systems can also migrate to smaller separations (i.e., below 200 au).

One way to explain our observation that there appears to be a preferential alignment between wide binaries and planets, and that the alignment is most prominent for binary stars with semimajor axes less than 700 au, is that close-in binaries form aligned with the protoplanetary disk, while more distant binaries form via a different mechanism that leaves their orbits misaligned. In this case, presumably this would imply that the close-in binaries in our sample formed mainly via disk fragmentation, while the more distant binaries formed via turbulent fragmentation and, more rarely, capture of the binary companion. If this were true, it would would necessitate that most binaries with semimajor axes less than 700 au formed from disk fragmentation, which would be in tension with the apparent ∼200 au cutoff for disk fragmentation (Tobin et al. 2016; El-Badry & Rix 2018), as well as the fact that some binaries below 200 au should have formed via turbulent fragmentation and later migrated inwards.

If the binaries in our system tend to all form aligned (which is not supported by the observed occurrence of turbulent fragmentation at separations similar to those seen in our study), then the binary orbits of the systems above the threshold ∼700 au could be sculpted by the dynamical influence of passing stars. Deacon & Kraus (2020) have shown that wide binaries are rare in open clusters, suggesting that open clusters are a highly dynamical environment. Before complete disruption of a binary system, the system can be dynamically perturbed, leading to the larger fraction of misaligned systems we see above ∼700 au. Such a dynamical sculpting need not take place in only clusters: Kaib et al. (2013) found that the orbits of binary systems can be altered by the galactic tide and passing stars even when systems are not in clusters.

However, as discussed in Section 6.1, observational evidence suggests that aligned systems form primarily below 200 au and thus the cutoff is probably not primarily due to a dynamical sculpting of the binary orbits.

The final possibility we consider is that some dynamical process has sculpted the population of planets in wide binaries into orbital alignment after formation. Wide-binary systems with separations similar to those probed in this work are expected to form primarily via turbulent fragmentation, where turbulence fractures the initial stellar core into multiple separate overdensities (Offner et al. 2010, 2016; Bate 2018). In turbulent fragmentation, the two stars tend to have a distribution of angular momentum vectors that is less aligned on average than disk fragmentation (Offner et al. 2016; Bate 2018). However, a larger portion of binaries that form via turbulent fragmentation are still expected to have aligned orbits than a completely isotropic distribution (Bate 2018). Indirect evidence for turbulent fragmentation in wide-binary systems includes outflow orientations (Lee et al. 2016) and direct imaging (Fernández-López et al. 2017).

The initial misalignment of some binary systems that form via turbulent fragmentation will cause the disk to precess by the mechanism described in Section 5.2, dissipate energy, and lead to an alignment between the disk and binary orbit. This scenario seems like it could explain both the observed alignment and the transition from predominantly aligned systems to more misaligned systems around ∼700 au, which corresponds to a disk precession timescale similar to the lifetime of protoplanetary disks.

We have presented evidence showing that the orbits of wide-binary stars are preferentially aligned with the orbits of close-in planets in those systems. Although the alignment appears to be statistically significant, it is based on a relatively small sample of planetary systems. Confirming the conclusions of this work with a larger sample of binary stars hosting transiting exoplanets would be highly valuable.

Fortunately, it should be straightforward to increase the size of our sample. Recently, El-Badry et al. (2021) released a sample of over a million wide-binary stars from Gaia EDR3, including 274 visual binaries found to have exoplanets by TESS as of 2021 February 1. With this new sample of binaries, it will be possible to triple the size of our current sample. In addition, with improved Gaia astrometry from future data releases, such as more systems with radial velocity and an astrometric acceleration term, the inclinations of systems can be better constrained. If our hypothesis that most of the binaries form via turbulent fragmentation and were dynamically sculpted into aligned systems is true, then the statistical significance of the difference in inclination distribution of binaries above and below 700 au should increase, and a clearer dividing timescale between aligned and randomly distributed systems should emerge.

A larger sample of binaries will also make it possible to subdivide the sample and look for correlations with planetary parameters. For example, initially misaligned wide-binary companions could help trigger the formation of hot Jupiters via Lidov–Kozai oscillations, so it would be interesting to see whether the inclinations of wide-binary companions in hot Jupiter systems are different from those in systems with small planets. Another important characteristic to investigate is the transiting planet system multiplicity. Compact multiplanet-systems' planets are susceptible to inclination oscillations from exterior companions that can change the planet's mutual inclinations enough to prevent them all from transiting (Becker & Adams 2017).

Another valuable avenue may be to make similar measurements for different samples of planetary systems. For example, ground-based, high-resolution imaging observations of binaries too close to be resolved by Gaia can help determine whether binaries with semimajor axes less than those probed in this work (a ≲ 100 au) also show preferential alignment. It may also be possible to increase the sample of particularly wide binaries (a ≳ 1000 au) by including planet candidates orbiting the more distant stars targeted by Kepler and K2 that reside in binary systems, which can be resolved by Gaia.

Finally, we note that as follow-up of TESS planets continues, the purity of the TESS planet-candidate sample will increase and systematics in TESS planet detection will become better understood. The ongoing TFOP observation campaign will continue identifying false positives among the TESS planet-candidate sample and increase the likelihood that any given candidate in the surviving sample is indeed a planet candidate. Identifying these false positives will remove noise from the distribution of wide-binary inclinations. As the systematics in TESS planet detection become better understood, it can be determined with a higher degree of confidence whether the observed alignment is affected by a bias in TESS's detection method.