Stellar Multiplicity Meets Stellar Evolution and Metallicity: The APOGEE View

Among the targets in our sample, most (62%) have three RVs that pass quality cuts, 22% have only two RVs, and 16% have four or more. A simple figure of merit to characterize these sparsely sampled RV curves is the difference between the highest and lowest measured RVs for each star, ΔRV_max = max(RV_n) − min(RV_n) (see Badenes & Maoz 2012 and Maoz et al. 2012, for discussions). In Figure 5 we show the values of ΔRV_max as a function of log g for the targets in our main sample. The maximum value of ΔRV_max measured by APOGEE is a strong function of log g, with MS stars showing values as high as ∼300 $\mathrm{km}\,{{\rm{s}}}^{-1}$, and stars near the tip of the RGB only going to ∼30 $\mathrm{km}\,{{\rm{s}}}^{-1}$.

Zoom In Zoom Out Reset image size

This trend in ΔRV_max versus log g can be best appreciated in the cumulative histograms shown in Figure 6, which highlight the tail of the ΔRV_max distribution above a few $\mathrm{km}\,{{\rm{s}}}^{-1}$. Here, we have divided our main sample of non-RC stars into eight log g bins, one for subgiant and MS stars (log g > 3.25), and seven for RGB stars down to log g ∼ 0. These bins have roughly the same number of stars (∼8000), except for the MS/subgiant bin, which has ∼19,000 (see Table 1). The MS/subgiant bin could in principle be resolved into smaller bins, but given the known limitations of the ASPCAP pipeline for high log g stars, we felt this was not justified. The histograms in Figure 6 show that there is a clear attrition of high ΔRV_max values as stars climb the RGB, and shows that the distribution of ΔRV_max in the RC closely resembles that near the tip of the RGB.

Zoom In Zoom Out Reset image size

Table 1.  Systems with ΔRV_max > 10 $\mathrm{km}\,{{\rm{s}}}^{-1}$

	log g	Median
Sample	Range	log g	N	${N}_{{\rm{\Delta }}{\mathrm{RV}}_{\max \gt 10\mathrm{km}{{\rm{s}}}^{-1}}}$	${N}_{{\rm{\Delta }}{\mathrm{RV}}_{\max \gt 10\mathrm{km}{{\rm{s}}}^{-1}}}/N$
MS/Subg.	3.250–4.932	4.203	19045	1051	0.0552
RGB 6	3.009–3.250	3.133	6211	200	0.0322
RGB 5	2.791–3.009	2.876	8388	183	0.0218
RGB 4	2.649–2.791	2.722	8385	114	0.0136
RGB 3	2.339–2.649	2.524	8412	158	0.0188
RGB 2	1.931–2.339	2.141	8393	135	0.0161
RGB 1	1.455–1.931	1.715	8378	85	0.0102
RGB 0	0.069–1.455	1.176	8366	62	0.0074
RC	2.250–3.250	2.766	15667	100	0.0064

Download table as:  ASCIITypeset image

The highest value of ΔRV_max in a large sample of stars will correspond to the edge-on (i = 90°) binary systems with the shortest possible orbital periods. This should be the critical period for RLOF,

$$\begin{eqnarray}&&{P}_{\mathrm{crit}}=2\pi { \mathcal R }(q)\sqrt{\displaystyle \frac{{R}^{3}}{{GM}}},\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal R }$(q) is the ratio between the radius of the Roche Lobe and the orbital separation (which can be approximated by 0.38 for q = 1, Eggleton 1983), and M and R are the mass and radius of the photometric primary (i.e., the star that contributes most of the flux in the APOGEE bands). For any period P, the peak-to-peak amplitude of the RV curve of an edge-on orbit, ΔRV_pp, is twice the semi-amplitude K,

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{RV}}_{\mathrm{pp}}=2K=\displaystyle \frac{2}{\sqrt{1-{e}^{2}}}{\left(\displaystyle \frac{\pi {GM}}{2P}\right)}^{1/3},\end{eqnarray} \tag{ 2 }$$

where e is the eccentricity. It follows from Equations (1) and (2) that the value of ΔRV_pp at P_crit for a circular orbit with q = 1 is uniquely determined by the mass and radius of the primary. For stars whose effective gravity g is measured directly, which applies to our APOGEE targets, these equations can be combined into a simple expression,

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{RV}}_{\mathrm{pp},{P}_{\mathrm{crit}}}=0.87{({GMg})}^{1/4}.\end{eqnarray} \tag{ 3 }$$

Therefore, a theoretical maximum for ΔRV_pp can be estimated for any measured log g just by assuming a mass, as shown by the solid lines overlaid on Figure 5. For M = 1 M_⊙, P_crit is 0.35 days in the MS (log g ∼ 4.5) and 2.1 years at the tip of the RGB (log g ∼ 0), and ΔRV_pp should vary between 400 and 30 $\mathrm{km}\,{{\rm{s}}}^{-1}$, as shown schematically in Figure 7. These theoretical maxima compare well with the ΔRV_max measurements in Figures 5 and 6, allowing for the fact that the observed systems have a distribution of primary masses, mass ratios, and inclinations, and their RV curves are sparsely sampled by APOGEE.

Zoom In Zoom Out Reset image size

Orbital eccentricities will also affect ΔRV_max, even though tidal interactions will circularize the orbits with the shortest periods, i.e., those contributing the largest ΔRV_max at any given log g, on timescales that are comparable to the evolutionary timescale on the RGB (Meibom & Mathieu 2005). For Sun-like MS stars, Raghavan et al. (2010) measured a circularization period of P_circ ∼ 12 days (log P_circ ∼ 1.1). Unfortunately, there is no comprehensive observational study of P_circ in RGB stars with varying log g. From the theory of tides, Verbunt & Phinney (1995) derived the value of P_circ as a function of the depth of the convective envelope in the RGB. Taking the depth of the convective envelope from the 1 M_⊙ model of solar metallicity by Choi et al. (2016), the expressions in Verbunt & Phinney (1995) predict a steady increase of P_circ with decreasing log g, with log P_circ ∼ 4.2 at the tip of the RGB (Figure 8). This provides a baseline theoretical scenario for orbital circularization in the RGB. The maximum eccentricity allowed at any period above P_circ is limited by angular momentum conservation, and can be approximated as

$$\begin{eqnarray}&&{e}_{\max }=\sqrt{1.0-{\left(\displaystyle \frac{{P}_{\mathrm{circ}}}{P}\right)}^{2/3}}\end{eqnarray} \tag{ 4 }$$

(Mazeh 2008), which has been shown to be in good agreement with Kepler observations of "heartbeat stars" (Shporer et al. 2016). Note that the period exponents in Equations (2) and (4) cancel out, therefore the maximum ΔRV_pp remains constant for periods above P_circ (horizontal blue and red lines in Figure 7). In practice, however, ΔRV_max values close to this theoretical upper limit for highly eccentric systems are hard to observe at periods much longer than P_circ because (1) there is a limit to the temporal baseline that can be probed in any RV survey, (2) systems with e ∼ e_max are rare, and (3) at high eccentricities it becomes increasingly difficult to capture the full dynamic range of RV with a sparse sampling, since most of the variation happens in a brief time interval close to periastron. In any case, orbital eccentricities will break down the simple relationship between P and ΔRV_pp shown in Equation (2) for periods longer than P_circ.

Zoom In Zoom Out Reset image size

To help visualize the effects of the APOGEE sampling on the distribution of ΔRV_max, we have added an inset to Figure 7 that shows the fraction of stars that could in principle probe the full value of ΔRV_pp, because they have temporal baselines longer than half a given period. The temporal baselines of individual APOGEE targets range between 0.8 and 1043 days, with a median of 40 days. The fraction of targets that can fully sample the maximum RV range as a function of period falls rapidly above log P ∼ 1.7, but it remains above 10% at log P ∼ 2.8, which is the value of P_crit at the tip of the RGB for systems with 1 M_⊙ primaries. There are no large systematic variations of the distribution of temporal baselines or the number of visits with either log g or metallicity in our APOGEE targets.

Even though the APOGEE data reduction pipeline reports RV errors below 0.1 $\mathrm{km}\,{{\rm{s}}}^{-1}$ (Nidever et al. 2015), Cottaar et al. (2014) found evidence that these errors might be underestimated by as much as a factor ∼3. The ΔRV_max distributions measured by APOGEE also indicate that either the average RV errors are larger than reported by the pipeline or there is an additional source of scatter in the individual RV measurements.

In a survey with a large number of sparsely sampled RV curves, the ΔRV_max distribution usually has a core of low values dominated by measurement errors and a tail of high values from bona fide RV variables (Maoz et al. 2012). Even without a detailed understanding of the RV errors, the transition between tail and core can often be determined empirically, and objects above the transition can be identified as bona fide multiple systems (see Badenes & Maoz 2012; Maoz & Hallakoun 2016). In the left panel of Figure 9 we show the ΔRV_max distribution in our APOGEE MS/Subgiant sample, as well as that of the stars in the highest and lowest [Fe/H] terciles within this sample. The core/tail transition can be clearly seen at ∼0.7 $\mathrm{km}\,{{\rm{s}}}^{-1}$, with the high and low-metallicity terciles having somewhat lower and larger transition values, respectively. Also shown in Figure 9 are models where a Gaussian distribution of RV errors is sampled with the APOGEE visits. A comparison to these models indicates a value of σ_RVerr, which is close to 0.2 $\mathrm{km}\,{{\rm{s}}}^{-1}$ for the MS/Subgiant sample. These ΔRV_max distributions, and the inferred RV errors, are representative of the majority of our log g samples. The one exception is shown in the right panel of Figure 9. The APOGEE targets with the lowest log g values (0.07 > log g > 1.46) have noticeably broader ΔRV_max cores, with a core/tail transition as high as ∼5.0 $\mathrm{km}\,{{\rm{s}}}^{-1}$ for the stars in the lowest [Fe/H] tercile. Some of these stars might have RV errors as high as σ_RVerr ∼ 1 $\mathrm{km}\,{{\rm{s}}}^{-1}$.

Zoom In Zoom Out Reset image size

This behavior of the RV errors is qualitatively consistent with the properties of the cross-correlation procedure used by the APOGEE data reduction pipeline to measure RVs, which is sensitive to parameters that affect the number, depth, and broadness of spectral lines, like log g and [Fe/H]. Quantitatively, the extent of our ΔRV_max cores seems to require larger RV errors than those reported by the pipeline (see Section 8.3 in Nidever et al. 2015). This might be due to a systematic underestimation of the RV errors by the pipeline or to the presence of a non-Gaussian tail of RV errors, but a more likely explanation is RV jitter, which can introduce a scatter of up to a few $\mathrm{km}\,{{\rm{s}}}^{-1}$ in the spectra of RGB stars, especially at low log g (Carney et al. 2003; Hekker et al. 2008). To ensure that RV errors, whatever their origin, will not affect our analysis, and that we will be able to compare stars of different log g and [Fe/H] without biases, we define a conservative threshold of ΔRV_max = 10 $\mathrm{km}\,{{\rm{s}}}^{-1}$ to identify a system as a bona fide multiple. This value is well within the tail of the ΔRV_max distribution even in the most unfavorable cases shown in Figure 9, but is low enough to yield a total of 2088 bona fide multiple systems, with enough systems in each log g sample for statistical analysis (see Table 1). From Equation (2), the minimum period required to get ΔRV_max = 10 $\mathrm{km}\,{{\rm{s}}}^{-1}$ is log P = 4.3 (2 × 10⁴ days, or 54 years), which is close to the expected circularization period for Sun-like stars at the tip of the RGB (Figure 8). However, the temporal baselines in APOGEE cannot probe such long periods (see Figure 7). The longest period we are sensitive to with this ΔRV_max threshold is that which can produce an RV shift of at least 10 $\mathrm{km}\,{{\rm{s}}}^{-1}$ in the longest baselines available in APOGEE (∼10³ days), which is about log P ∼ 3.3, or 5.5 yr.

The relationship between ΔRV_max and log g shown in Figures 5 and 6 can be understood qualitatively through the equations introduced in Section 3.2 and the interplay between stellar multiplicity and stellar evolution. After ∼1 M_⊙ primaries exhaust H in their cores and leave the MS, they climb the RGB and their log g drops from ∼4.5 to ∼0 as their radii increase from ∼1 to ∼170 R_⊙. For those in multiple systems, the maximum value of ΔRV_max allowed for P = P_crit drops from ∼400 to ∼30 $\mathrm{km}\,{{\rm{s}}}^{-1}$ (Equation (3)). Because we cannot find any multiple systems in APOGEE with ΔRV_max values above these limits (i.e., to the right of the solid lines in Figure 5), all systems with P < P_crit must have been removed from the sample by some efficient process. This removal of short P systems also results in a lower number of stars observed at all values of ΔRV_max, as seen in Figure 6, due to the projection effect of random orbital inclinations (multiply Equation (3) by a factor sin i that is randomly distributed between 0 and 1). After core He ignition, the stars settle on the RC and their radii decrease again to ∼10 R_⊙, but their ΔRV_max distribution remains similar to that of the larger stars at the tip of the RGB, because their short-period companions have already been removed during shell H burning.

A detailed quantitative evaluation of this scenario, including constraints on multiplicity fractions and period distributions, would require forward modeling of the APOGEE ΔRV_max measurements within a hierarchical Bayesian scheme, taking into account all the relevant stellar properties and the details of the mass distributions and tidal interactions for the entire sample (e.g., see Maoz et al. 2012; Walker et al. 2015). We leave that analysis for future work, and here we examine the main physical implications of our observations using a simpler method.

We generate artificial populations of stars that can be sampled with the APOGEE epochs using a Monte Carlo code. Our code assumes that all photometric primaries are 1 M_⊙ (see Figure 4 and accompanying discussion), that the distribution of mass ratios is flat (a good first approximation for short-period companions to Sun-like stars, Moe & Di Stefano 2017), and that orbital inclinations are random (i.e., the distribution of cos i is uniform). For each run, we choose a MS multiplicity fraction f_m and²⁴ an effective gravity log g. We adopt the period distribution of Raghavan et al. (2010); see Figure 1 and accompanying discussion), which we truncate at the value of P_crit that corresponds to the chosen log g (Equation (1)). We assume that all systems with P < P_crit have lost their companions and can be considered single. We calculate P_circ from the theory of Verbunt & Phinney (1995) and the 1 M_⊙, [Fe/H] = 0 model of Choi et al. (2016) (Figure 8), and assume that all orbits with P_crit < P ≤ P_circ are circular. For longer periods, the eccentricity is drawn from a uniform distribution (Moe & Di Stefano 2017) between 0 and e_max (Equation (4)). We generate N systems with these parameters, each of which is sampled with the epochs (number of visits and time lags between visits) from a random APOGEE target, with the orbital phase of the first visit drawn from a uniform distribution between 0 and 2π. Thus, the Monte Carlo code captures the main physics affecting the values of ΔRV_max with only three free parameters: f_m, log g, and N.

In Figure 10 we compare our Monte Carlo simulations with the fraction of targets with ΔRV_max ≥ 10 $\mathrm{km}\,{{\rm{s}}}^{-1}$ observed by APOGEE in each of the log g samples and in the RC catalog (Table 1). The attrition of high ΔRV_max (short P) systems as stars climb the RGB is clearly seen in the APOGEE data, but the trend seems to break at log g ∼ 2.7 (sample RGB 4 in Table 1). This is the region of the HR diagram where the density peaks of RGB and RC stars overlap (see Figure 2), which suggests that the lower fraction of high ΔRV_max systems in this sample might be due to a contamination from unidentified RC stars with lower ΔRV_max values. This is not surprising, since the RC catalog of Bovy et al. (2014) is designed to be pure, rather than complete (see Price-Jones & Bovy 2018).

Zoom In Zoom Out Reset image size

Setting aside the issue of the RC contamination, our Monte Carlo simulations reproduce the observed behavior of ${N}_{{\rm{\Delta }}{\mathrm{RV}}_{\max }\geqslant 10\mathrm{km}{{\rm{s}}}^{-1}}/N$ as a function of log g quite well for RGB stars. The value observed in the RC is close to that predicted by our simulations at the tip of the RGB, as expected in the framework of our physical model. In the RGB, we obtain the best match to the observations for f_m = 0.35. The RGB sample with the highest log g (∼3.1) and the MS/subgiant sample seem to require higher values of f_m, but these deviations should be interpreted with caution. The APSCAP pipeline was not designed for stars with high log g, particularly at low T_eff, so the fitted stellar parameters in these samples might be subject to systematic uncertainties (García Pérez et al. 2016; note the mismatch between data and models in this region of Figure 2). Even taking the fitted stellar parameters at face value, it is clear from Figure 2 that these samples must have a primary mass distribution that is significantly different from the bulk of RGB and RC stars in APOGEE. Our Monte Carlo code assumes a single value of f_m for all stars, which is not a valid approximation if there is a broad enough range of primary masses (Lada 2006; Moe & Di Stefano 2017). Furthermore, the mass ratio distribution, MS period distribution, and tidal circularization model in our code are only appropriate for 1 M_⊙ primaries. In view of this, the Monte Carlo simulations for the two highest log g samples might not reflect the underlying properties of the APOGEE targets. We will discuss our measurement of f_m in more detail in Section 4.

The relationship between stellar multiplicity and metallicity is a controversial topic. Some numerical simulations of star formation show that metallicity should have a strong effect on multiplicity, with lower-metallicity clouds showing higher fragmentation rates and smaller initial separations for binary systems (Machida 2008; Machida et al. 2009), but others predict no significant impact of metallicity for values of [Fe/H] between −1.0 and 0.5 (Bate 2014). Observationally, Gao et al. (2014, 2017) found an inverse correlation between [Fe/H] and multiplicity fraction in field F-, G-, and K-type dwarfs with multiple RVs from the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST, Cui et al. 2012). Hettinger et al. (2015) found the opposite trend in a smaller sample of F-type dwarfs from SEGUE. However, these studies are based on optical spectra with low resolutions (R ∼ 1800), which makes it hard to confidently measure small RV shifts.

Here, we address this topic with the high-resolution IR spectra from APOGEE. In Figure 11 we break down the ΔRV_max distribution in each of the samples from Table 1 into [Fe/H] terciles. The boundaries and median values of these terciles are shown in Figure 12. There is some variation of the tercile boundaries with log g, but for the purpose of this discussion we can consider [Fe/H] ≲ −0.5 as the low tercile, and [Fe/H] ≳ 0.0 as the high tercile.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In all non-RC samples from the MS/subgiants to the tip of the RGB, the ΔRV_max distribution in the lowest [Fe/H] tercile is clearly above that in the highest [Fe/H] tercile. This is true at all values of ΔRV_max above our 10 $\mathrm{km}\,{{\rm{s}}}^{-1}$ threshold, provided there are enough objects to have small Poisson noise (i.e., ${N}_{{\rm{\Delta }}{\mathrm{RV}}_{\max }\geqslant 10\mathrm{km}{{\rm{s}}}^{-1}}/N$ above a few times 10⁻³ for most of our samples). The ratio between the fraction of systems with ΔRV_max ≥ 10 $\mathrm{km}\,{{\rm{s}}}^{-1}$ in the high and low [Fe/H] terciles is shown in Figure 13, together with the 1σ probability interval obtained from comparing random terciles in our Monte Carlo simulations. The ratio is roughly between 2 and 3 for all non-RC samples except the sample at log g ∼ 2.1, where it reaches 4.5. The only sample that shows no significant metallicity effect is the RC. This could be related to the fact that the spread of metallicities in the RC catalog is smaller than that in the other APOGEE samples (see Figure 12). Stellar structure in the RC is itself sensitive to metallicity (Girardi 2016), so very-low-metallicity He-burning stars ([Fe/H] ∼ −1) will actually be on the horizontal branch instead of the RC. This introduces a selection effect that leads to a narrower [Fe/H] distribution, where a multiplicity trend with metallicity can be easily washed out by errors in the APSCAP stellar parameters.

Zoom In Zoom Out Reset image size

The most salient features of Figures 11 and 13 can be interpreted in simple terms. It is clear that low-metallicity and high-metallicity stars in APOGEE must have different multiplicity characteristics. Either the MS multiplicity fraction or the MS period distribution (or both) must be metallicity-dependent, though the clear segregation at all values of ΔRV_max seen in Figure 11 favors a metallicity-dependent multiplicity fraction (see discussion in Section 3 of Maoz et al. 2012). Provided that the MS period distribution remains roughly lognormal, with a sharp truncation at P_crit for RGB stars of all metallicities, the multiplicity fraction of APOGEE targets in the lower-metallicity tercile should be a factor 2–3 higher than that in the higher -metallicity tercile, to reproduce the ratios shown in Figure 13. If f_m ∼ 0.35 for the mean ([Fe/H] ∼ −0.3), the multiplicity fraction for metal-poor ([Fe/H] ≲ −0.5) stars should be close to 0.46, and that for metal-rich ([Fe/H] ≳ 0.0) stars should be close to 0.23. The rising trend in the ratio of high ΔRV_max systems for log g ≲ 3 might be related to the shift toward lower metallicity in the [Fe/H] distributions in that log g range seen in Figure 12.

Although a metallicity-dependent multiplicity fraction does provide the simplest explanation for the APOGEE observations, we emphasize that our simple analysis can only offer a partial glimpse into what is likely a complex situation. For example, it is well known that T_eff, which we have ignored, is strongly correlated with [Fe/H] in the RGB (see Figure 12 in Holtzman et al. 2015), with lower [Fe/H] stars having higher T_eff at constant log g. Since [Fe/H] is an intrinsic property of the star that is set at birth, we have chosen to present our results as a function of [Fe/H], rather than T_eff; but see Gao et al. (2017) for a discussion of the interplay between the two parameters. Another factor that could in principle influence the ΔRV_max distributions is the mass of the primary, but we have found no large systematic differences in the masses determined by Ness et al. (2016) for the high- and low-metallicity terciles in our samples, so this effect, if present, must be small.

The best-fit value of f_m ∼ 0.35 obtained by comparing our Monte Carlo simulations to the APOGEE observations of RGB stars is ∼60% lower than the canonical multiplicity fraction for Sun-like MS stars (0.50 ± 0.04 for log P < 8, which is equivalent to 0.56 over the entire period range; Moe & Di Stefano 2017). A direct comparison between these two numbers, however, is not trivial. Our code assumes the Raghavan et al. (2010) period distribution (lognormal, with a peak at log P ∼ 5.0, Figure 1), but the APOGEE baselines are only sensitive to systems in the tail of this distribution, below log P ∼ 3.3. Stellar multiplicity surveys usually have a small number of systems with such short periods, so the value of f_m is more uncertain in this period range than the integral over all periods. After correcting the Raghavan et al. (2010) sample for completeness, Moe & Di Stefano (2017) listed 32 systems with log P < 3, which implies a $\sqrt{N}/N=18 \%$ uncertainty in the value of f_m from Poisson statistics alone, much larger than the 8% given for the entire sample. The tidal circularization model in our Monte Carlo code, which is derived from theory alone, might be inaccurate, and this could impact the number of high ΔRV_max systems produced by our simulations. In Monte Carlo runs where the value of P_circ shown in Figure 8 is divided by 2 and 5, we find that the relative increases in ${N}_{{\rm{\Delta }}{\mathrm{RV}}_{\max }\geqslant 10\mathrm{km}{{\rm{s}}}^{-1}}/N$ are ∼3% and ∼5%, respectively. Finally, Moe & Di Stefano (2017) discussed the probability that the photometric primary in a SB1 binary might in fact be the mass secondary in a post-common envelope system with a white dwarf companion (see their Section 8.3 and Figure 29). The fraction of short-period (log P < 3.3) SB1 systems with white dwarf companions depends on the age of the stellar population, and is ∼15% at 1 Gyr and ∼30% at 10 Gyr. Because the RGB stars in APOGEE have a median age of 4.6 Gyr (Ness et al. 2016), this fraction should be around ∼20% for the bulk of our sample, but the impact on high-metallicity and low-metallicity stars will be different (see below). The combination of the effects we have discussed here (∼18% from small sample sizes at short periods, ∼5% from uncertainties in the tidal circularization model, and ∼20% from white dwarf companions) is comparable to the ∼60% difference between our measurement and that of Moe & Di Stefano (2017). Therefore, we cannot claim that our measurements of the multiplicity fraction are mutually inconsistent.

The metallicity trend shown in Figures 11 and 13 will also be affected by the presence of white dwarf companions to short-period SB1 binaries. However, the increase of the fraction of white dwarf companions from 15% to 30% between high-metallicity and low-metallicity populations estimated by Moe & Di Stefano (2017) cannot explain the factor 2–3 increase in the fraction of high ΔRV_max systems seen in APOGEE. In principle, a systematic difference between the masses of stars in the high- and low-metallicity terciles of each log g sample could explain our results, but this difference would have to be very large. If the close (log P ≲ 4) binary fraction scales as M^0.7 (Moe & Di Stefano 2017), the median mass of the stars in the low-metallicity terciles would have to be a factor ∼2.6 higher than the median mass of the stars in the high-metallicity terciles to increase the fraction of systems with ΔRV_max > 10 $\mathrm{km}\,{{\rm{s}}}^{-1}$ by a factor of 2. In our RGB samples, the median masses of stars in the high- and low-metallicity terciles measured by Ness et al. (2016) vary by a few tenths of a M_⊙ at most.

We conclude that the inverse correlation between metallicity and multiplicity fraction seen in APOGEE must be real. Our result is in qualitative agreement with the works of Gao et al. (2014, 2017) that used low-resolution spectra of MS stars, and the findings of Yuan et al. (2015) from photometric studies, but it is in apparent contradiction with the results of Hettinger et al. (2015) from nearby F-type dwarfs. Moe & Di Stefano (2013) determined that the binary fraction at very short periods (log P < 1.3) for 10 M_⊙ primaries does not vary by more than ∼30% across the [Fe/H] range between −0.7 and 0.1 by studying early-type eclipsing binaries in the Milky Way and the Magellanic Clouds. Older studies based on small samples of spectroscopic orbits by Latham et al. (2002) and Carney et al. (2005) also reported no significant correlation between metallicity and multiplicity. However, caution must be used again when comparing such different surveys. Because they are not looking at the same stars, and they are not sensitive to the same periods, it is hard to establish to what degree these observational studies might be mutually inconsistent. To date, ours is the only study based on high-resolution spectra that span all evolutionary phases between the MS and core He-burning for Sun-like stars in a large fraction of the Milky Way disk, with a sample size orders of magnitude larger than any previous attempt. The recent re-analysis of the MS spectra in APOGEE by El-Badry et al. (2018b) has provided an independent confirmation of the inverse correlation between metallicity and multiplicity that we report here.

A detailed interpretation of the metallicity effect that we have discovered is outside the scope of the present work. Naively, one might say that our results support the conclusions of Machida (2008) and Machida et al. (2009) regarding the impact of metallicity on star formation, and that studies where metallicity has little or no effect on stellar multiplicity at birth can be ruled out. (e.g., Bate 2014). These multiplicity statistics set "at birth" will be subject to dynamical evolution through orbital capture and disruption, which can be quite different in different environments such as the thin disk, the thick disk, and the halo. However, numerical simulations imply that the hard binaries we study here (log P ≤ 3.3) should be relatively unaffected by this (Kroupa & Burkert 2001; Goodwin & Kroupa 2005; Kroupa & Petr-Gotzens 2011). More complex processes involving stellar mergers can also affect the relationship between multiplicity, metallicity, and stellar age (e.g., Jofré et al. 2016; Izzard et al. 2018).

Regardless of its origin, a metallicity-dependent multiplicity fraction has profound implications for the rates of interacting binaries like SNe Ia and neutron star mergers, which are observed in host galaxies with a wide range of redshifts and therefore metallicities (Zahid et al. 2013). Binary population synthesis calculations for these objects often assume that the multiplicity statistics in the MS (i.e., the initial conditions) are not metallicity-dependent (e.g., Claeys et al. 2014; de Mink & Belczynski 2015), an assumption that should be revised in light of our results. The dynamical stability of circumbinary planets will also be affected by a higher fraction of short-period binaries at low metallicities (Jaime et al. 2014), which could be related to the well-established paucity of planets around low-metallicity stars (Johnson et al. 2010).

4.2.1. Future Prospects