systems (Schwarz et al. 2016). 1 It is therefore evident that planets in multi-star systems are not unusual. What is less well-known is the prevalence of planetary systems in wide binaries, and, moreover, whether the evolution of such planetary systems may lead to alterations in stellar surface chemistry to a degree that would inflate observed chemical differences between the stellar components of wide binaries.

From a spectroscopic follow-up of the 131 radial velocitydetected candidate exoplanet systems known as of July 2005, Raghavan et al. (2006) estimate that 23% of exoplanet-hosting stars have one or more wide companions. Larger multiplicity studies, which aim to identify new companions to exoplanet host stars from the complete sample of >4000 known exoplanet hosts (e.g., Mugrauer & Neuhäuser 2009;Roell et al. 2012;Abushattal et al. 2019), find that the multiplicity rate of exoplanet hosts is ∼15%. However, the separation of the wide stellar binary significantly impacts the formation of circumstellar (S-type) planets. Moe & Kratter (2019) find that binary stars with separations of 10 AU host close planet companions at ∼15% of the occurrence rate of single stars, though WBs with separations >200 AU have a negligible impact on the formation of close planets. Additionally, Deacon et al. (2016) find that a wide stellar companion, with a separation greater than 3000 AU, does not impact the occurrence rate of short-period (<300 days) planets orbiting solar-like stars.

Current formation scenarios for stellar binaries imply that stars in such systems are born at the same time from the same protostellar cloud, and the stars in a binary ought to exhibit identical chemical compositions. Checking that hypothesis is a central means by which to test models of star formation and evolution. However, it is critical to employ binaries for which no later alteration in surface chemistry is expected. Thus, binaries composed of little-evolved (e.g., main sequence) components are preferable. Moreover, the stars should not be close to one another, both to avoid any possibility of previous or current mass transfer or other interaction, but also for the practical reason that the stars must be separated sufficiently for accurate spectroscopic analysis. Given the above preferences, wide main sequence binaries provide an ideal opportunity for such experiments.

The overall consistency in the metallicity of the stellar components of wide binaries had been explored in detail for a number of systems (e.g., Gratton et al. 2001), while a few studies ventured beyond overall metallicity to study individual elements, including Fe (Desidera et al. 2004(Desidera et al. , 2006)), V (Desidera et al. 2006), and Li (Martín et al. 2002). However, recently, the picture has dramatically changed through two more comprehensive analyses (Andrews et al. 2019;Hawkins et al. 2020) performed on larger wide binary samples (containing dozens of pairs) and testing a large number of chemical elements in these systems. While bearing exquisitely on the question of the "identical twin" nature of wide binary stars, both studies were motivated by an alternative goal, which was to assess the viability of chemical tagging as a means to identify Galactic field stars with their siblings from a common origin. Wide binaries have been suggested, as a complement to open cluster studies (Bovy 2016;Ness et al. 2018), to be a potentially unique sample with which to calibrate and test chemical tagging (De Silva et al. 2007;Andrews et al. 2018bAndrews et al. , 2019)). The results of these wide binary tests are promising for the prospects of chemical tagging as a useful tool for stellar population studies.

For example, Andrews et al. (2019) used high resolution, high signal-to-noise infrared spectra from the Apache Point Observatory 1 https://www.univie.ac.at/adg/schwarz/multiple.html Galactic Evolution Experiment (APOGEE; Majewski et al. 2017) survey to test the chemical consistency of 62 stars in 31 wide binaries. A key (and necessary) advantage of this study is that both components were observed and analyzed in a completely consistent manner. Andrews et al. (2019) found that the differences in abundances of individual elements are consistent with being due to measurement uncertainties alone. The authors conclude that the consistency of individual elements between components is too similar to imply anything other than a common origin for the wide binary pairs in this list.

Meanwhile, the Hawkins et al. (2020) study of 25 wide binaries derived the abundances of 24 chemical species spanning most nucleosynthetic pathways (including light, α, odd-Z, ironpeak, and neutron capture elements) from high resolution optical spectra. Their analysis found remarkable chemical homogeneity between the stars in the binaries, for the most part consistent to the level of their very precise measurement uncertainties (e.g., to less than 0.02 dex for [Fe/H]), but with differences for some elements reaching as high as 0.10 dex. However, for elements such as e.g., Al, Ca, Ti, these differences are greater than the reported measurement uncertainties.

The observed chemical differences that are slightly above the errors may either be due to underestimated errors or due to some alteration in stellar chemistry occurring even in wide binaries. Some processes have been suggested that might lead to abundance differences in wide binaries; in particular, differences in the history of planet formation or accretion in two otherwise identical stars may leave signatures in chemical abundances (e.g., Li in Carlberg et al. 2012;Oh et al. 2018, refractory elements). Elemental abundances of wide binaries that host planets have been studied in detail (e.g., Mack et al. 2014;Teske et al. 2015Teske et al. , 2016;;Liu et al. 2018). These studies find significant differences (typically on the order of 0.05 dex, but up to ∼0.2 dex, in metallicity, Oh et al. 2018) in the abundance patterns of components of wide binaries that host (or may have hosted) planets, which has been attributed to (1) the fact that forming more gas giants or rocky planets leads to an overall depletion of metals in the gas that will eventually accrete onto the host star (Biazzo et al. 2015;Ramírez et al. 2015) and (2) the accretion of planetesimals onto host stars (Mack et al. 2014).

Though wide binaries have been suggested as potential candidates for chemical tagging efforts (De Silva et al. 2007;Andrews et al. 2018bAndrews et al. , 2019)), if the processes described above do in fact impact a large number of wide binary pairs, including pairs which have undergone such processes in chemical tagging training sets might decrease the precision to which those studies can identify stars with similar abundances. Therefore, by identifying and removing from the training sets those wide binary pairs which have formed companions or accreted planetestimals will improve chemical tagging efforts.

In the present work, we aim to identify triple and higherorder multiple candidates and to explore the possibility for the detection of the chemical imprints of exoplanet formation, utilizing precise RV measurements from the Sloan Digital Sky Survey (SDSS-IV; Blanton et al. 2017) APOGEE spectroscopic survey for 37 wide binaries selected from the wide binary sample presented in El-Badry & Rix (2018). In Section 2, we outline the process we use to identify our sample of wide binaries in the APOGEE-Gaia cross-matched catalog. We detail the two methods utilized for detecting candidate close companions in Section 3. Section 4 gives detailed notes about individual triple-and higher-order systems, as well as a brief discussion of the abundance differences of all 37 wide binaries, versus just those wide binaries found to host stellar-

Companions in Stellar Wide Binaries 3 and substellar-mass companions. Finally, in Section 5 we provide some conclusions and ideas for future directions.

The APOGEE-Gaia Wide Binary Catalog is a subset of the wide binary sample described in El-Badry & Rix (2018), which was constructed by searching Gaia DR2 (Gaia Collaboration et al. 2016, 2018) for pairs of stars within 200 pc of the Sun, with positions, proper motions, and parallaxes consistent with being gravitationally bound. The El-Badry & Rix (2018) catalog contains >50,000 wide binaries, consisting of either main-sequence (MS) or white dwarf (WD) components, with separations of 50 s < 50, 000 AU.

The catalog provides Gaia measurements of proper motion, parallax, magnitude, RV (where measured), and their associated errors, for each star in the pair, as well as the physical separation of the pair and the binary class (i.e., MS/MS, MS/WD, or WD/WDwhere "MS" and "WD" mean main sequence and white dwarf, respectively). The APOGEE survey (Majewski et al. 2017) is a highresolution (R ∼ 22, 500), high signal-to-noise ratio (S/N > 100), infrared (1.51 µm to 1.69 µm), spectroscopic survey using twin spectrographs (Wilson et al. 2019) attached to the 2.5-m Sloan telescope at APO in the Northern Hemisphere (Gunn et al. 2006) and the du Pont telescope at Las Campanas Observatory (LCO) in the Southern Hemisphere. The latest APOGEE-2 data release 2 (SDSS DR16; Ahumada et al. 2020;Jönsson et al. 2020) comprises spectra for ∼437,000 unique stars, along with stellar parameters (effective temperature T eff , surface gravity log g, metallicity [M/H], αelement abundance [α/M], etc.) and abundances for over 20 individual chemical elements (for some stars), derived from the APOGEE Stellar Parameter and Chemical Abundance Pipeline (ASPCAP; García Pérez et al. 2016;Shetrone et al. 2015;Smith et al. 2021). In addition to elemental abundances, APOGEE provides multi-epoch RV measurements, which have an internal precision better than ∼0.1 km s -1 (Nidever et al. 2015) with accuracies ∼0.35 km s -1 (over all T eff and log g; Anguiano et al. 2018).

The two catalogs (El-Badry & Rix 2018 and APOGEE-2) were cross-matched adopting a 1 arcsec positional tolerance. The APOGEE sample is made up of red giants stars (∼70%) and MS stars (∼30%); stellar remnants (e.g., WDs) are too faint to be observed in the ∼1 hr APOGEE visits (in the H-band) and are therefore excluded from the APOGEE input catalog. For this reason, the cross-match of the two catalogs contains only MS/MS wide binary pairs. Additionally, El-Badry & Rix (2018) required both stars to have a measured Gaia G BP -G RP color, and to have well-resolved photometry (i.e., uncontaminated by nearby sources, where the degree of contamination is assessed using the Gaia phot_bp_rp_excess_factor). As a consequence of these requirements, the El-Badry & Rix (2018) catalog has an effective resolution limit of ∼2 arcsec. Note, however, that APOGEE fibers for sources observed on a single SDSS fiber plugplate have a much larger collision limit-71.5 arcsec for APOGEE-2N and 56 arcsec for APOGEE-2S (Wilson et al. 2012(Wilson et al. , 2019))-so that the closest separation of any wide binary pair in this catalog is ∼56 arcsec. In summary, the catalog used in the present work contains only those wide binaries for which we have APOGEE spectra (including stellar parameters and abundances from ASP-2 Available at: https://www.sdss.org/dr16/irspec/ CAP) for both components of the binary, as well as high-quality astrometry and photometry from Gaia.

RVs were not used in the selection of candidate binaries in the El-Badry & Rix (2018) catalog, but prove a useful check that binaries are gravitationally bound. In genuine binaries, the RVs of the two stars should agree within ∼a few km s -1 (within the RV error; El-Badry & Rix 2018). For this reason, to be considered wide binaries we require candidate wide binaries to have APOGEE RVs consistent to within 5 km s -1 (i.e., RV primary -RV secondary ≤ 5 km s -1 , where RV primary and RV secondary are the S/N-weighted average radial velocities for the primary and secondary stars in the WB, respectively). In general, this cut may remove stellar triples in which a close inner sub-companion induces large RV variations on one of the wide binary components; however, this criterion is also highly efficient at removing random alignments (see e.g., Andrews et al. 2018a,b), and is imperative for producing a catalog with minimal contamination. In this particular sample, the cut only removes one WB pair; the stars in this pair had only 1-2 visits, and one star had an APOGEE warning flag set (BRIGHT_NEIGHBOR), indicating that the average RVs for both stars in the pair may not have been reliably measured. We also select only stars with low APOGEE median visit RV error, σ v,med ≤ 0.2 km s -1 , and high S/N (≥ 25) to remove stars for which low-quality spectra may cause errors in the derived parameters and RVs. These criteria ensure minimum contamination of the wide binary sample by false-positive, random alignment contaminants.

In addition, we remove stars flagged with the ASPCAPFLAGS STAR_BAD or that have no result from ASPCAP, and pairs where one or both stars do not have calibrated [M/H] measurements, to ensure that we are considering only stars with reliably determined stellar parameters and abundances.

The resulting APOGEE-Gaia Wide Binary Catalog contains 37 genuine MS/MS wide binaries (shown in T eff -log g in Figure 1 and listed in Table 5) with separations between 300 s < 50, 000 AU, all within 200 pc of the Sun. All of the stars in the catalog have been observed exclusively by the APOGEE-2N instrument at APO; this is by and large due to the relatively few stars that had been observed by the APOGEE-2S instrument at LCO at the time the DR16 database was assembled. The sample spans the temperature range from <3500 K up to 6500 K, and confirms that all binaries in the catalog are composed of MS/MS pairs. We compare the Gaia proper motions and parallaxes for the stars in the 37 wide binaries in Figure 2. In most cases, measurement uncertainties are smaller than the data points. The dotted line in each panel of the figure shows the moving median (in windows of 2500 stars) proper motion ∆µ and distance ∆d differences for all pairs in the El-Badry & Rix (2018) catalog. Using the Gaia astrometric solution for each star, El-Badry & Rix (2018) requires that pairs must have a distance difference ∆d < 2 × s and that the difference in proper motion ∆µ of the two stars must be consistent with a bound Keplerian orbit (i.e, ∆µ ∝ s -1/2 ). Based on the construction of the catalog, the correlation between the proper motions and parallaxes (or distances) of the wide binary components with increasing separation is to be expected. As most pairs show better consistency in the astrometry and RVs than the stars in the full El-Badry & Rix (2018) catalog, our sample is free from contamination by random alignments, and contains only gravitationally bound pairs.

The left panel of Figure 3 the components of the 37 wide binaries in our catalog. The scatter in the APOGEE RV measurements (V scatter in the APOGEE catalog) are shown as the error bars, but are typically small compared to the size of the data points. Note that stars with only one visit in the APOGEE catalog have a V scatter equal to zero; however, we propose a 0.072 km s -1 single-visit error floor, which we assume as the V scatter for systems with only one visit (see Section 3.2 for additional discussion of the derived error floor). In the center panel, for the 12 cases where both components of the wide binaries have RVs measured by Gaia, a comparison of the Gaia RVs is shown. That the RV differences between most primary and secondary components of the 37 systems are consistent with RV primary -RV secondary 1 km s -1 , which is the median |∆RV| for all wide binaries in the El-Badry & Rix (2018) catalog, supports the veracity of the classification of these systems as genuine wide binaries. Finally, the right panel of Figure 3 shows the difference between the RVs measured by the two surveys-for the 12 pairs where both components of the wide binaries have RVs measured by both surveys-as a function of Gaia G magnitude. The RV measurements are in good agreement between the surveys.

In Figure 4 we show the comparison between metallicities of both stars in the wide binary pairs (blue). (Note, though some recent WB catalogs use metallicity as a selection criteria for genuine WB pairs, e.g., Godoy-Rivera & Chanamé 2018, the El-Badry & Rix 2018 catalog does not.) The Spearman-ρ correlation coefficient for the metallicities of all pairs in the catalog is 0.80 (99% CI = 0.58 to 0.91, derived using the jackknife technique), where a value of unity corresponds to a perfect correlation (zero is equivalent to no correlation). The relatively consistent measured [M/H] between components of the wide binary sample provides further evidence that our wide-binary sample is likely free of contamination. We report the root mean square error (rms∼0.13) about the expected ∆[M/H] = 0 line in this plot as well; this value will be important for comparing with the rms of other sub-samples of this dataset.

It has been shown (e.g., Souto et al. 2020) that the derived APOGEE DR16 metallicities for cool stars (particularly, late-K/M dwarfs) are too metal-poor. The right panel of Figure 4 also shows a comparison of the metallicity differences in the pairs, as a function of temperature; it does appear that this known effect causes larger scatter in the abundance differences for the coolest stars/pairs in our sample. Specifically, the two pairs with the largest abundance differences (i.e., those with |∆[M/H]| ∼ 0.4) contain the coolest dwarfs. Additionally, ASPCAP uncertainties do not account for any systematic trends in abundances as a function of stellar parameters; therefore, comparing the abundances of stars with large differences in surface gravity or effective temperature may over estimate the real metallicity differences of these stars. Andrews et al. (2019) show that, for stars with repeat observations as part of multiple APOGEE fields that are passed through the ASPCAP pipeline as two distinct objects (but with necessarily identical T eff and log g), typical metallicity differences are ∼0.03 dex (their Table 1). They also find that the difference in metallicity between WB components approaches ∼0.03 dex for pairs with similar T eff . Though the expected abundance differences certainly increase for stars with significantly different effective temperatures and surface gravities, we adopt the typical metallicity difference observed in the APOGEE repeat observations as a more appropriate measure of the metallicity uncertainties in this sample.

In this catalog, because all stars are dwarfs, differences in log g are relatively small ( 0.5 dex) and do not appear to influence the differences in the derived abundances. On the other hand, the differences in effective temperature between WB components do appear to have a significant impact on the observed abundance differences for these dwarf stars. If we look at just the subsample of pairs for which both stars in the pair have T eff > 4000 K and the difference in temperature between those stars is |∆T eff | < 200 K (red circled points in Figure 4) the rms is reduced to rms∼0.06, and is more consistent with the metallicity errors (0.03 dex).

In this section we outline the multiple methods utilized to detect candidate close companions in our catalog of wide binary systems.

An unresolved binary system comprising two identical stars (i.e., equal in mass) has the same color but twice the luminosity of an equivalent single star; such a system, comprising two equal-mass MS stars, will appear in the color-magnitude diagram 0.753 mag brighter, irrespective of the wavelength bands used (e.g., Hurley & Tout 1998). We show in Figure 5 the Gaia G BP -G RP color against the absolute Gaia G magnitude. The grey Hess diagram in the background shows the density of stars from the cross-match of APOGEE with the full Gaia DR2 catalog, including MS and giant stars. We overplot a 5 Gyr, solar-metallicity Dartmouth isochrone (solid line; Dotter et al. 2008;Evans et al. 2018), as well as the expected location of a corresponding binary sequence, represented as the same isochrone shifted 0.753 mag in absolute magnitude.

The figure confirms that all of the stars in our sample are on the MS, but also identifies four stellar mass companions to components of four separate wide binaries that lie near the equal-mass binary sequence. These four systems are discussed briefly in Section 4.1.

In addition to stellar parameters, the APOGEE reduction pipeline provides RV measurements with derived visit-level uncertainties typically better than ∼0.1 km s -1 for the APOGEE-2N instrument (Nidever et al. 2015), where an APOGEE "visit" is defined as the spectrum of a source from a single fiber plugplate's observation on a night (∼1 hr of exposure). The total survey number of APOGEE visits scheduled for a star depends on its H-band brightness, with fainter targets requiring more visits to obtain a survey-required minimum total S/N of 100; at the completion of the survey, most stars are expected to have a minimum of three visits to accomplish this goal. Though the DR16 catalog does not represent the final APOGEE database, all of the stars in the APOGEE-Gaia Wide Binary Catalog are brighter than H ∼ 12.5, and happen to already have combined combined spectra with high S/N, even though they may be composed of fewer than three visit spectra. Uncertainties derived by the APOGEE RV pipeline for DR16 (Nidever et al. 2015) are known to be systematically underestimated, as noted by Cottaar et al. (2014) and Badenes et al. (2018); if not properly accounted for in our analysis, underestimated errors will cause the spurious detection of binaries, whereas overestimated errors may hide truly variable systems. For the sample of stars in APOGEE DR16 (1) which are not telluric standard stars, (2) with five or more visits, and (3) with median visit-level RV uncertainty σ v,med < 0.5 km s -1 from the pipeline, Figure 6 shows the observed RV scatter, V scatter , versus σ v,med . The sample analyzed here contains stars with a broad range of temperatures, surface gravities, and H-band magnitudes, and so should be broadly applicable to any population of stars in the APOGEE survey. For any given single star (i.e., one that is not RV-variable), we expect V scatter to be comparable to the RV uncertainties; however, for the bulk of the stars the pipeline-derived RV uncertainties appear to be underestimated by a factor of 2-3×, as shown by the red one-to-one line (i.e., V scatter = σ v,med ) in the figure . Using the trend observed in Figure 6, we characterize the true uncertainties σ v,true more accurately by rescaling the visit RV uncertainties σ v via the best-fit to the RV scatter as a function of median RV uncertainty according to the parameterization

where σ v and σ v,true have units of km s -1 . These coefficients given in the equation represent the best fit to the median V scatter values over 25 equal count bins in σ v,med .foot_1

The true visit velocity errors may also depend on stellar parameters, log g and T eff . These errors are found to be constant with effective temperature (see Appendix A), but depend more strongly on surface gravity, likely because stellar surface oscillations (i.e., "jitter")-a significant contributor to the RV scatter-are inversely dependent on log g (Hekker et al. 2008). See Appendix A for further discussion on the trends remaining in the RV errors following

the application of Equation 1. Because we are only dealing with MS stars (log g > 4) in the present analysis (yielding a sample spanning only ∼1 dex in surface gravity), the impact of the surface gravity on the velocity errors is minimized, though systems with low-metallicity ([M/H] -0.5, Badenes et al. 2018) will tend to have larger errors. For many of the systems presented in this work, Equation 1 yields σ v,true ∼ 0.1 km s -1 , which approaches the lower-limit (0.072 km s -1 ) set by the equation. This σ v,true is smaller than the typical RV accuracies derived by e.g., Anguiano et al. (2018) or Mazzola et al. (2020); however, the systems presented in this work are all MS stars, whereas at least some fraction of the samples presented in Anguiano et al. (2018) and Mazzola et al. (2020) are composed of giant stars. Because surface oscillation/stellar jitter (which can be observed as changes in RV) are inversely correlated with surface gravity (Hekker et al. 2008), the expected APOGEE uncertainties are also expected to be inversely correlated with log g. Therefore, for the sample of stars in this work (all of which have log g>4) APOGEE RV uncertainties should be small compared to the typical RV accuracy for a star in the APOGEE dataset, which is composed of ∼2/3 giant stars. The systems presented in this work are also moderately brighter (median H-band magnitude ∼10.5 mag) than the average APOGEE target (median H-band ∼12.5 mag), so the RVs are expected to be more precise than the typical APOGEE star.

To ensure that the true errors supplied by Equation 1 are neither underestimated nor overestimated, we calculate the reduced chisquared statistic,

where n is the number of APOGEE visits and ν = n-1 is the number of degrees of freedom in the data for the star, v i and σ v,true,i are a single visit RV and corresponding error calculated from Equation 1, and v is the weighted average velocity of the star. The probability of exceeding a given χ 2 ν is P( χ 2 , ν), under the assumption that all stars are not RV variable binaries (i.e., the null hypothesis that all stars in the sample do not have close binary companions) and have Gaussian errors. The histogram of these probabilities for all stars in this WB catalog with two or more high-quality RV measurements in the APOGEE survey is shown in Figure 7.

If no close binaries are present in the catalog and the errors are properly inflated via Equation 1, the probability distribution should be uniform over all probabilities; if the errors are overestimated or underestimated the histogram should be strongly biased towards higher or lower probabilities, respectively. If the errors are properly estimated and close binaries do exist in the sample, we should see a sharp spike in the lowest probability bin, 0 < P( χ 2 , ν) < 0.01. 4The latter is precisely what is observed in the included figure; given these results, we can conclude that for this sample of stars the errors resulting from Equation 1 are not under-or overestimated. It should be noted that analyses of samples of APOGEE stars that span a larger range of log g or that include giant stars should include an Furthermore, for the widest systems (s 3000 AU), the primaries and secondaries should have RVs within ∼0.2 km s -1 (unless they have a close companion, most of which should not) to be gravitationally bound. For the widest systems in the sample, we do find that a majority of the pairs (12 out of 17 pairs with s 3000 AU) have RV differences (plus the quadrature sum of the inflated errors) less than 0.2 km s -1 . Of the five systems with RV differences >0.2 km s -1 , one is identified in this paper as a stellar triple (Section 4.1.1) and one is a system with a reported substellar-mass companion (Section 4.1.2).

With the properly inflated RV uncertainties, the data are sufficient for the detection of RV oscillations expected from relatively short-period companions down to a few Jupiter-masses (∼2 M Jup ; e.g., Troup et al. 2016;Price-Whelan et al. 2020). For stars with ≥ 8 visits, over a sufficient temporal baseline (i.e., longer than the period of the orbit), these high-quality RV measurements are suitable to detect RV variability and to derive Keplerian orbital fits.

For stars with fewer than 8 RV measurements, there are many orbit models that are consistent with the data, which leads to a challenging degeneracy of solutions spanning a large range of possible orbital parameters. The Joker (Price-Whelan et al. 2017, 2020) provides a Monte Carlo sampler, intended for systems with sparse RV measurements, that produces samplings in orbital parameters for data that include as few as 3 visits; however, for 5 or fewer visits, the returned samples are highly multi-modal, with the samples tending to form a harmonic series. For this reason, we limit our orbit fitting using The Joker to those stars in the APOGEE-Gaia Wide Binary Catalog with 6 or more RV measurements with APOGEE. In the APOGEE-Gaia Wide Binary Catalog 16 stars in 11 wide binary pairs have 6 or more visits, and so are eligible for Keplerian orbit fitting with The Joker.

To search for substellar mass companions, we generate 2 26 prior samples for the nonlinear parameters described in (Price-Whelan et al. 2020), evaluate the marginal likelihood of each sample, and rejection sample to produce a minimum of 256 posterior samplings in the nonlinear parameters. Additionally, we reject any samples with e > 0.970, which corresponds to the largest eccentricity of any known exoplanet, HD 20782b.foot_3 If fewer than the requisite 256 posterior samples are returned (indicating that the data constrain the orbit very well), we initialize a Markov chain Monte Carlo (MCMC) run, which generally returns a unimodal set of samples. Lewis et al. We report the stellar and orbital parameters from APOGEE and The Joker, respectively, as well as the minimum companion mass m sin i for each candidate. The minimum companion mass is derived using the returned posterior samplings from The Joker and primary stellar masses derived via the Torres et al. (2010) relations for MS stars. We define a planetary-mass companion as one with a minimum mass below the deuterium-burning limit (0.013 M ∼ 13.6 M Jup ), and a stellar-mass companion as any candidate above the hydrogenburning limit (0.080 M ∼ 83.8 M Jup ) (e.g., Burrows et al. 2001;Auddy et al. 2016). A brown dwarf (BD) companion is any candidate that falls between these two limits.

Orbit fitting with The Joker was performed for all stars with 6 or more APOGEE RV measurements. From this, we detect and report the presence of four substellar-mass companion candidates orbiting their WB-component hosts in S-type orbits (i.e., in a circumstellar orbit, where the planet orbits one of the two stars, in contrast to P-type orbits where the planet orbits the entire binary in circumbinary motion): one hot-Jupiter (i.e., planet mass) candidate and three candidates-with either an unconstrained orbital period or eccentricity-that have minimum masses that range between the planet and BD regimes. These four substellar-mass candidates orbit stars contained in only two wide binary systems, indicating that we have detected candidate quadruple systems with companions on each star. These systems are described in detail in Section 4.1.

Note that we do not report systems for which The Joker returns highly multi-modal posterior samples, where a "mode" is defined as a period at which more than 10 samples are returned. We only report those systems with relatively few modes (<5 modes) in the returned samples.

The following are notes on individual systems for which we have detected stellar-(Section 4.1.1) and substellar-mass companion candidates (Section 4.1.2) using the methods previously outlined. The systems are identified by the 2MASS ID of the primary component of the wide binary, as defined by El-Badry & Rix (2018); additional information about the wide binary systems-including observational details about the systems (e.g., number of visits passing our quality cuts, average radial velocity, etc.)-can be found in Table 5. For the candidate systems presented in Section 4.1.2, we have also examined the cross-correlation functions (CCFs) of the individual visits with the best-matching RV template, and have checked the individual visit spectra to confirm that the RVs derived from the APOGEE pipeline are accurate. We find no evidence that the pipeline has produced flawed RVs for any of the visits included in the orbital analysis of these stars. Please see additional discussion on this analysis at the end of the Section 4.1.2.

In each of the following cases where the presence of a stellar companion orbiting one of the components of a wide binary is assumed from the Gaia color-magnitude diagram (CMD, see Section 3.1), the systems do not have a sufficient number of APOGEE visits for Keplerian orbit fitting, so the orbital periods, eccentricities, etc. of these systems remain unknown. Therefore, for all of the triple systems described in this section, the companion is assumed to be of nearly equal-mass to its host. J03023780+0019430: triple system with a stellar companion in a circumstellar orbit about the secondary -In the Gaia CMD, the secondary of this wide binary pair (G BP -G RP = 1.21, M G = 5.37) clearly lies on the equal-mass binary MS, while the primary component of the wide binary (G BP -G RP = 1.21, M G = 6.15) falls on the MS. Therefore, we infer the presence of a stellar-mass companion around the secondary of this wide binary system, but not around the primary. The mass of the stellar companion orbiting the secondary can be assumed to be of equal-mass to the secondary, m ∼ 0.8 M .

J03093344-0053352: triple system with a stellar companion in a circumstellar orbit about the primary -The primary of this wide binary pair (G BP -G RP = 0.97, M G = 4.56) lies on the equal-mass binary sequence. Its secondary (G BP -G RP = 0.84, M G = 4.61) lies on the MS. Again, we infer the presence of a stellar-mass companion around the primary component of this wide binary from the Gaia photometry, but do not infer the presence of any companion around the secondary. The likely mass of the companion orbiting the primary is m ∼ 1.0 M . J16342492+4146352: triple system with a stellar companion in a circumstellar orbit about the primary -In the Gaia CMD, the primary of this system (G BP -G RP = 1.79, M G = 7.09) lies on the binary MS, but the secondary (G BP -G RP = 0.96, M G = 5.26) does not. As above, the presence of the companion, with mass m ∼ 0.6 M orbiting the primary is assumed from the Gaia photometry, but we do not infer the presence of any companion around the secondary. J19434227+4926329: triple system with a stellar companion in a circumstellar orbit about the secondary -The secondary of this wide binary pair (G BP -G RP = 0.91, M G = 4.38) lies on the equal-mass binary MS. The primary component of the wide binary (G BP -G RP = 1.11, M G = 5.91) does not lie on the binary MS. The mass of the stellar companion orbiting the secondary is assumed to be m ∼ 1.0 M .

In the following four cases of the detection of a substellar companion, the derived m sin i from all samples returned by The Joker is 2 M Jup . Again, we do not report all systems with the requisite 6 APOGEE RV visits here; we only report those systems for which The Joker returns samples with fewer than 5 period modes. J12530085+2734189: quadruple system with two substellar companions, one orbiting each the primary and secondary in Stype orbits -The primary component of this wide binary has eleven visits in the APOGEE survey meeting our quality cuts. From the maximum a posteriori (MAP) sample returned following MCMC analysis by The Joker, we find the RV variations observed are indicative of a companion with a minimum-mass of 8.3 M Jup in a 24 day orbit (i.e., a hot-Jupiter). However, depending on the true eccentricity of the system (shown by the colorbar in Figure 8), the minimum companion mass may fall anywhere in a range from 5 M Jup (for the lowest eccentricity orbit) to 40 M Jup (for the highest eccentricity orbit), potentially placing it in the BD mass regime. For this reason, the companion is classified only as a substellar-mass candidate (i.e., we do not specify planet-vs. BD-mass); additional RV measurements will better constrain the eccentricity and minimum mass of the system. Additional discussion of the likelihood of this candidate being a planet vs. a BD is included at the end of this section.

The secondary, J12525132+2735131, has twelve high-quality visits in the APOGEE survey. From the MAP sample returned by a standard analysis by The Joker (i.e., no MCMC), we find the RV

variations indicate a companion with a minimum-mass of 6.3 M Jup in a 35 day orbit. The surviving samples are multimodal, with a second, longer period mode at ∼386 days, with samples at the longer period tending to have slightly higher minimum masses (∼15 M Jup ).

As is the case for the companion candidate orbiting the primary of the system, additional RV visits are required to constrain the eccentricity and period of the system, and therefore accurately constrain the minimum mass of the companion. For the shortest period orbits, the minimum companion mass ranges between 5 and 10 M Jup (i.e., a planetary-mass companion), but for the longer period orbital modes, the companion mass may exceed ∼60 M Jup , placing it within the BD-regime. Again, the candidate is classified only as a substellar-mass object.

The orbital parameters for the close companions to the stars in this wide binary are given in Tables 1 and2, respectively. For the unimodal system (i.e., with samples all falling at the same period), the orbital parameters given are those of the orbit. The returned samples are approximately normally distributed in period P, eccentricity e, and barycentric velocity v 0 , so the reported errors are given by the standard deviation of the returned samples in each of those parameters. For parameters where the returned samples do not follow a Gaussian distribution, velocity semi-amplitude K and minimum companion mass m sin i, we report asymmetric errors, where the upper and lower errors are given as the difference between the MAP sample and the third and first quartile values, respectively, for those parameters. For the system with multiple period modes, we provide one sub-table for each period mode-beginning with the mode containing the MAP sample-where the parameters and errors are calculated as described for the unimodal system, for only the samples in that period mode. The phased-folded orbits are shown in Figure 8 and9. Stellar parameters from the APOGEE survey for the wide binary components are given in Table 5.

The importance of this quadruple system in the context of previously discovered hierarchical systems will be discussed at the end of this section. M Jup J12582326+2630091: quadruple system with two substellar companions, one orbiting each the primary and secondary in Stype orbits -The primary component of this wide binary has eleven visits in the APOGEE survey. From the MAP sample returned following an MCMC analysis, we find the RV variations observed are indicative of a companion with a minimum-mass of 3.4 M Jup in a 12 day orbit. Even with the large spread in possible orbital eccentricities given by the returned samples, the minimum companion mass for all samples fall within the planet-mass regime. However, if the companion is in a highly-inclined orbit (i.e., the system is viewed close to face-on) about the primary, the true mass m of the companion may exceed the 13.6 M Jup -limit, above which we would define the companion to be a BD; in this case, if i 15 • (where i = 0 • is face-on and i = 90 • is edge-on), m exceeds 13.6 M Jup . Because we cannot constrain the inclination of the system given these data, we classify the companion only as a planet candidate.

The secondary, J12582222+2630162, has thirteen good visits in the survey. From the MAP sample returned by a standard analysis by The Joker, we find the RV variations indicate a companion with a minimum-mass of 4.8 M Jup in a 123 day orbit. However, the surviving samples are multimodal in period, with modes as short as ∼20 days and as long as ∼200 days. Additional RV visits are required to constrain the period and eccentricity of the system, and therefore accurately constrain the mass of the companion. Depending on the true period and eccentricity of the system, the minimum companion mass may fall in a range from 2 M Jup to 50 M Jup , making it a possible BD candidate. Again, the companion can be classified only as a substellar-mass candidate.

The orbital parameters for the companions are given in Tables 3 and4, respectively, and the stellar parameters from the APOGEE survey for the companion hosts are reported in Table 5. The phasedfolded orbits are shown in Figures 10 and11.

These two quadruple systems are highly extraordinary among the stellar multiples systems with companions that have so far been discovered. Of the ∼75 binary systems known to have at least one companion in an S-type orbit, only five of those systems have companions orbiting both stars -HD 133131 Abc Bb (all planet mass companions), HD 20782 b / HD 20781 bc (all planet mass), HD 41004 Ab Bb (1 planet, 1 BD mass), WASP-94 Ab Bb (both planet mass), XO-2 Sbc Nb (all planet mass). Though additional RV followup is necessary to constrain the masses of the candidate companions reported here, if confirmed, the two candidate quadruple systems presented would make up more than 25% of all known quadruple systems containing two substellar companions in S-type orbits about each wide binary component (Schwarz et al. 2016).

Only one of these previously discovered systems has companions in both the planetary and BD-mass regimes (as defined here), and none have two BD-mass companions. Only HD 41004 Ab Bb is known to have a confirmed planetary-mass companion (HD 41004 Ab, 2.5 M Jup ; Zucker et al. 2004) orbiting the other. Furthermore, of the ∼160 confirmed substellar companions in binary-and higher-order systems with reported masses, only five of those companions have minimum masses in the BD regime (Schwarz et al. 2016).

For the systems presented here, with orbital periods 100 days and minimum masses that span the planet and BD mass regimes, the candidate companions are several times more likely to be planets than BDs, given that planets are more abundant (e.g., under the notion of the existence of a BD desert - Grether & Lineweaver 2006;Csizmadia et al. 2015-however, cf. Troup et al. 2016). In the end, if any of the three candidates that we currently classify as substellar-mass objects are confirmed to be BD-mass companions by additional RV follow-up, it would represent an increase in the known number of BDs in stellar multiples by as much as 50%.

In the context of previous studies, the prodigious identification of hierarchical systems here is surprising. We note that Fernandes et al. (2019) 2020, Fontanive et al. 2019, and references therein) find that wide binaries with separations >200 AU do not have a significant influence on the formation or occurrence rates of close planets, and that while substellar-companion hosts do show an enhancement in their WB fraction compared to stars that are not close companion hosts, it is at most a factor of ∼2× enhancement for the most massive, closest-separation substellar companions. In contrast, we find an enhancement in the occurrence rate of close substellar companions in our WB stars of more than an order of magnitude compared to the above studies-i.e., four stars with substellar companions out of 16 stars for which orbit fitting is performed in Section 4.1.2, or an occurrence rate ∼25%, instead of the nominal ∼2%. The probability of finding four close binary systems in a sample of 16 stars when we might expect to find 0.2-0.3 close substellar-mass binaries is P ∼ 0.0001, according to Poisson statistics; thus, either our yield of such systems is vastly overestimated or the expected occurrence rate is vastly underestimated.

We discern no part of our survey selection function that should act to enhance our occurrence rate. The most likely source of systemic bias in our analysis would seem to lie in the interpretation of the RV time series data. In particular, as is well known, orbit fitting of poorly sampled RV variations is prone to error. While there is nothing we can do about the existing phasing or number of RV visits within the APOGEE survey for our stars, it is possible to investigate whether our solutions to these data may have led to improper or unlikely combinations of companion mass, period, and/or orbital eccentricity. Of particular concern are the effects that RV outliers may have in driving a particular orbital solution. To this end, we have investigated the individual RV measurements for each of the four candidate companion hosts by checking carefully the RV CCFs and by comparing the visit-level spectra corrected by their calculated visit RVs to the rest-frame velocity, and we find no reason to doubt the veracity of any individual RV measurement. We have also performed analyses of orbital samples returned by The Joker where we have removed either the most extreme outlier RV (e.g., the observation at phase ∼0.95 with RV ∼ -47.3 km s -1 in Figure 8) or single random RVs, and we find that the solutionswhile multimodal-by and large, have relatively large velocity semiamplitudes (K 0.5 km s -1 ) and indicate that the observed RV variations may be induced by companions with minimum masses >2 M Jup . Clearly the best way to resolve this dilemma is to obtain more high-quality data utilizing the APOGEE spectrograph; additional data for the systems explored here, as well as multi-epoch data for WB pairs in the southern hemisphere, may be available in the forthcoming final APOGEE database.

These candidate quadruple systems are remarkably uncommon examples of multiple systems, a surprising success of a relatively small sample. However, Hwang et al. (2020) propose that the occurrence rate of multiple close companion hosts in wide binaries (like the candidate quadruple systems we report here) is higher than the occurrence rate of close companion hosts in the field (e.g., Mayor et al. 2011;Zhou et al. 2019). This is explained by the fact that systems hosting close companions-specifically those containing hot Jupiters-tend to have higher metallicities (e.g., Fischer & Valenti 2005); if that hot Jupiter host has a wide stellar companion, assuming that the wide companion has a very similar metallicity (shown to be true here, and in Andrews et al. 2019 andHawkins et al. 2020), we can assume that the wide companion has a higher chance of also hosting a close, hot Jupiter companion. The sample of known quadruple systems with close, substellar-mass companions is small though, so additional work on larger samples of such systems are required to confirm this finding.

For the remaining wide binary components-those without a close companion assumed from Gaia photometry, nor a companion detected from Keplerian orbit fitting (or too few visits for orbit fitting to be applied)-given the RV data from APOGEE or the photometry from Gaia, it is not currently possible to ascertain whether close-companions might exist without additional information. In future works, we intend to use the APOGEE ∆RV max (e.g., Badenes et al. 2018;Mazzola et al. 2020) and/or the Gaia renormalized unit weight error (RUWE), which indicates high astrometric noise (e.g., Belokurov et al. 2020), to detect unresolved companions in these systems; however, quality cuts used in making the El-Badry & Rix (2018) catalog may have removed such systems from our sample. It is unlikely that these systems will receive substantially more APOGEE visits during the remainder of the survey (i.e., some may be observed 1-2× more, but few will receive enough additional visits to reach the minimum 6 visits required here), but the APOGEE RVs may provide the basis for a longer monitoring campaign to search for massive close companions orbiting WB components.

We compare the differences in metallicity and abundances for eleven individual APOGEE-measured elements for the two components of the 37 wide binaries in the sample. We first focus on the consistency of the metallicity between components of the wide binaries, as shown in Figure 12, for those pairs that contain close companions compared to those pairs that do not. As shown previously, in Figure 4, the metallicities of the population of pairs in the catalog show a very strong, positive correlation, with Spearman-ρ = 0.80 (99% CI = 0.58 to 0.91, rms∼0.13, and rms∼0.06 when limited to pairs with |∆T eff | < 200 K). For the sample of stars without any companions detected (red points in Figure 12), the strength of the correlation between wide binary metallicities is Spearman-ρ = 0.82; while this is larger than the Spearman-ρ coefficient for the entire population, it is within the CI of the population, so is not statistically significant. The sample of systems with detected stellar-or substellar-mass companions show a weak positive correlation, with Spearman-ρ = 0.37 (blue points in Figure 12); however, because the sample contains only 6 data points (4 wide binaries with detected stellar companions and 2 wide binaries with detected substellar companions), the correlation coefficient is not well constrained, and should not be taken to mean that the metallicities between wide binary components with companions are less correlated than metallicities between wide binary components without companions. Additionally, the rms of the sub-sample of pairs without detected close companions (red in Figure 12) is neither significantly different from, nor improved, as compared to the full sample.

Figure 13 expands this comparison of elements to individual elemental abundances of C, O, Mg, Al, Si, K, Ca, Ti, Mn, Fe, and Ni. Because we've already shown that the metallicities of each pair are correlated in Figure 12, we therefore expect to see similar correlations between the pairs for [X/H] (where where X can be substituted for the element of interest). Here, instead of showing elemental abundances with respect to H, we show abundances with respect to metallicity, such that

(3) ASPCAP-measured abundances that have been excluded from this analysis either have too few pairs where both components have measured values (e.g., Na), or have large measurement uncertainties that are not useful for the comparison (e.g., Co). Based on the entire population of wide binaries in the APOGEE-Gaia Wide Binary Catalog-including those with companions-most of the elements shown here have moderate, positive correlations (aside from Ti, which shows a negative correlation, though this is one of the less reliably elements derived by ASPCAP). Like the trends observed for the metallicity of the pairs, for a majority of the elements shown here, after removing systems with close companions from the sample, the correlation between individual abundances does not change statistically. Similarly, systems with detected companions tend to show moderate to weak, statistically insignificant correlations, which should be expected to be positive; the fact that some elements show negative correlation coefficients (indicating that they are weakly anticorrelated) is likely due to the small sample size.

The chemical abundance correlations for the sample without detected companions and the sample with detected companions do not differ significantly from one another, and therefore no conclusions can be drawn about the impact of massive companions on the chemistry of wide binaries. The sample explored here is relatively small; a large sample of wide binaries with companions may reveal significant chemical differences between wide binary components. In addition to the small sample sizes in the metallicity and abundance analysis presented here, the large scatter in the metallicity differences are likely dominated by systematic trends in abundance with stellar parameters (as discussed in Section 2.1), making it impossible to differentiate the chemistries of stars with and without massive, close companions.

From a relatively small catalog-the APOGEE-Gaia Wide Binary Catalog, which provides Gaia photometry and astrometry, as well as APOGEE stellar parameters, visit-level radial velocities, and chemical abundances-we identify four triple systems, containing close stellar companions, and two candidate quadruple systems, each containing a total of two close substellar-mass companions. All systems in the catalog show very similar abundances, though the removal of systems containing close companions from the full sample does not improve the abundance consistencies observed for systems without close companions, as predicted.

Of the four substellar candidates, (ignoring inclination) one certainly falls within the planetary mass regime, and all four orbit in quadruple systems-a wide binary with each wide component having a closer, massive companion in an S-type orbit. None of these candidates are previously reported in the literature, and require further follow-up with RV observations (as part of the APOGEE survey or with another instrument) to be confirmed, since many of the systems have unconstrained periods or eccentricities, which impact the assumed minimum mass. Though relatively few quadruple systems with two substellar companions, each in close orbits about a component of a wide binary, have been confirmed (<10 systems Schwarz et al. 2016), this small sample presents two additional candidates. The orbital solutions and companion masses for these candidates must be carefully confirmed, as they indicate an order of magnitude enhancement in the expected occurrence rate of >1 M Jup companions with periods <300 days (occurrence rate of ∼25%, rather than the expected ∼2%).

We do not observe any statistically significant differences in the chemistries of wide binaries which have detected companions, compared to the chemistries of those without detected companions.

Given the small dataset and the magnitude of the abundance errors achieved with APOGEE, the influence or correlation of additional close companions is not discernible from the assumed identical chemistries of these pairs. A larger catalog of wide binaries hosting close companions, perhaps assumed statistically from the Gaia RUWE rather than derived from RV variation, with precise metallicities and abundances will be needed to better quantify any such effects. (which shows all stars in APOGEE DR16 that are not telluric standards, having five or more visits and σ v,med < 0.5 km s -1 ), the pipeline-derived, median RV errors are severely underestimated, with V scatter σ v,med for stars over all temperatures, surface gravities, metallicities, and apparent magnitudes. These errors are most significantly underestimated for the brightest stars in the sample (i.e., those typically with the highest SNR). For H 11, the errors are only underestimated by a factor of 2-3× (a similar result is also presented in Cottaar et al. 2014), but for brighter stars (say H ∼ 8), σ v,med is a factor of ∼20-30× smaller than the RV scatter. For future APOGEE data releases (DR17 in late 2021; see Holtzman et al. in prep.), a new RV pipeline has been implemented, and this significant underestimation of the RV errors of bright stars has been resolved.

By applying Equation 1 to the pipeline-derived RV errors, the true errors σ v,true are inflated substantially, and have a minimum RV error of 0.072 km s -1 . The ratio of V scatter to the true RV errors as a function of the explored stellar parameters and observational properties is shown in Figure A2. This expression for the error inflation has corrected the global 2-3× underestimation of σ v,med relative to the RV scatter that was observed across all parameter space, as well as the observed trend with decreasing H magnitude.

Even with the correction applied from Equation 1, there remains a significant trend in the ratio of V scatter to σ v,true for stars with decreasing surface gravity (i.e., giant stars). This is likely a real astrophysical effect caused by the fact that the RV scatter for these low-log g stars is inflated by stellar surface oscillations manifesting as velocity "jitter". These oscillations are expected to be inversely correlated with surface gravity (Hekker et al. 2008), as shown in Figure A2 for stars with log g 1.5. Though we do not correct for the effects of stellar "jitter" in this work, if the APOGEE RVs are to be utilized in any form (including the RV scatter, visit-level RVs, or RV errors), an additional "jitter" term should be added in quadrature to the true velocity errors derived in Equation 1(see e.g., Price-Whelan et al. 2020).