COMPANIONS TO APOGEE STARS. I. A MILKY WAY-SPANNING CATALOG OF STELLAR AND SUBSTELLAR COMPANION CANDIDATES AND THEIR DIVERSE HOSTS

To fully understand the types of companions to which we are sensitive, we need a clear understanding of dependencies of the RV precision on stellar parameters. Therefore, we created an empirical model of the RV precision based on the primary derived stellar parameters (${T}_{{\rm{eff}}}$, $\mathrm{log}g$, [Fe/H]) and the S/N for each visit of the star:

$$\begin{eqnarray}\mathrm{log}{\sigma }_{v} & = & 1.56+(4.87\times {10}^{-5}){T}_{{\rm{eff}}}+0.135\mathrm{log}g\\ & & -0.518[\mathrm{Fe}/{\rm{H}}]-(5.55\times {10}^{-3}){\rm{S}}/{\rm{N}},\end{eqnarray} \tag{ 1 }$$

where S/N is the signal-to-noise ratio of the visit spectrum from which the RV measurement was derived, and ${\sigma }_{v}$ is the RV measurement error in m s⁻¹. This model was determined by fitting a linear function of each parameter of interest using all APOGEE stars with at least 8 visits, excluding stars used as telluric standards and stars that have unreliable stellar parameters. The left panel of Figure 2 displays two of the stronger effects on RV error: [Fe/H] and S/N per visit. The effects of $\mathrm{log}g$ and ${T}_{{\rm{eff}}}$ are illustrated in the right panel. These effects are closely related to the strength and number of absorption lines in the spectra. For a typical solar metallicity ([Fe/H] = 0) giant (${T}_{{\rm{eff}}}=4000\;{\rm{K}}$, $\mathrm{log}g=3$) and typical solar metallicity dwarf (${T}_{{\rm{eff}}}=5000\;{\rm{K}}$, $\mathrm{log}g=4.5$) stars with ${\rm{S}}/{\rm{N}}=10$, we derive a typical RV precision of ∼130 m s⁻¹ and ∼230 m s⁻¹, respectively per visit. These are the random RV uncertainties reported by the APOGEE pipeline, and are likely to be underestimates of the true uncertainty (see Appendix B.1).

Zoom In Zoom Out Reset image size

RV measurements from observations with ${\rm{S}}/{\rm{N}}\lt 5$, as well visits that produced failure conditions in the RV pipeline, were not included in the final RV curves submitted to the orbit fitter. This reduced the number of stars for which Keplerian orbits could be attempted from 14,840 to 9454 stars.

Likely RV variable stars were selected using the following statistic:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{RV}}=\mathrm{stddev}\left(\displaystyle \frac{{\boldsymbol{v}}-\tilde{v}}{{{\boldsymbol{\sigma }}}_{v}}\right)\geqslant 2.5,\end{eqnarray} \tag{ 2 }$$

where ${\boldsymbol{v}}$ and ${{\boldsymbol{\sigma }}}_{v}$ are the RV measurements and their uncertainties, and $\tilde{v}$ is the median RV measurement for the star. The criterion was motivated by the false positive analysis presented in Appendix A.1.2. There are also several additional pieces of information that we used to pre-reject stars that would have resulted in poor or erroneous Keplerian orbit fits. Therefore we also removed stars with the following criteria:

• 1.  
  The system's primary must be characterized with reliable stellar parameters (${T}_{{\rm{eff}}},\mathrm{log}g$, [Fe/H]), so the ASPCAP STAR_BAD flag must not be set for the star. Derivations of the RVs and the physical parameters of the system both rely on reasonable estimates of the stellar parameters of the host star.
• 2.  
  The star cannot have been used as a telluric standard. These stars are selected for APOGEE observation for their nearly featureless spectra, so it is likely that RVs derived for these stars are unreliable and would lead to false positive signals.
• 3.  
  The combined spectrum from which the stellar parameters and RVs were derived cannot be contaminated with spurious signals due to poor combination of the visit spectra, so the SUSPECT_RV_COMBINATION flag must not be set for the star. This criterion also catches the double-lined spectroscopic binaries (SB2s) that would have resulted in poor stellar parameters, RVs, or orbital parameters from our current pipelines.

This preselection reduced the number of stars for which Keplerian orbit fits were attempted from 9454 to 907. This is not to say the stars excluded do not have any sort of RV variation, but the false positive interpretation cannot be ruled out for these stars, so we elected not to include them.

Before any further stellar properties are estimated, we divide the stars in this sample into 5 classes defined by the following crteria:

• 1.  
  Pre-main Sequence (PMS): Stars flagged in APOGEE as young stellar cluster members (IC 348 and Orion).
• 2.  
  Red Clump (RC): Stars in the APOGEE RC Catalog (Bovy et al. 2014).
• 3.  
  Red Giant (RG): Stars not selected as RC or PMS stars with

  $$\begin{eqnarray*}{T}_{{\rm{eff}}} & \lt & 5500\;{\rm{K}},\\ \mathrm{log}g & \lt & 3.7+0.1[\mathrm{Fe}/{\rm{H}}].\end{eqnarray*}$$

  The second relation was derived by mapping the $\mathrm{log}g$ of the base of the giant branch as a function of [Fe/H] from Dartmouth isochrones (Chaboyer et al. 2008) for typical ages expected of APOGEE giants.
• 4.  
  Subgiant (SG): Stars not selected as RC or PMS stars with

  $$\begin{eqnarray*}{T}_{{\rm{eff}}} & \gt & 4800\;{\rm{K}},\\ \mathrm{log}g & \geqslant & 3.7+0.1[\mathrm{Fe}/{\rm{H}}],\\ \mathrm{log}g & \leqslant & 4-(7\times {10}^{-5})({T}_{{\rm{eff}}}-8000\;{\rm{K}}).\end{eqnarray*}$$

  The second relation only applies for ${T}_{{\rm{eff}}}\lt 5500\;{\rm{K}}$. The third relation was determined by the $\mathrm{log}g$ at the highest ${T}_{{\rm{eff}}}$ of Dartmouth isochrones at a variety of ages and [Fe/H], roughly mapping the main-sequence turnoff (MSTO), and fitting a liner function to these points.
• 5.  
  Dwarf (MS): Any star that does not fit into any of the above categories are classified as MS stars.

These classifications are saved for the catalog, and illustrated in Figure 4.

Zoom In Zoom Out Reset image size

In addition to stellar parameters, we need an estimate of the stars' bolometric magnitudes to compare to the bolometric luminosities we calculate and use in the following derivations of the masses and radii of the primary stars. We adopt the extinction coefficient, A_K from the APOGEE targeting data (Zasowski et al. 2013). If the APOGEE targeting A_K is not populated or is less than zero, then we adopt the WISE all-sky K-band extinction. In the rare case ($\lt 1\%$ of stars run through the apOrbit pipeline) that neither quantity is available, we assume A_K = 0, and flag the star. The extinction-corrected K_s magnitude is then ${K}_{0}={K}_{s}-{A}_{K}$. We derived the bolometric correction to the 2MASS K_s band from Dartmouth isochrones:

$$\begin{eqnarray}{\mathrm{BC}}_{K} & = & (2.7+0.15[\mathrm{Fe}/{\rm{H}}])\\ & & -(25+0.5[\mathrm{Fe}/{\rm{H}}]){X}^{2-0.1[\mathrm{Fe}/{\rm{H}}]}{e}^{-X}\end{eqnarray} \tag{ 4 }$$

for PMS, dwarf and SG stars, where $X=\mathrm{log}{T}_{{\rm{eff}}}-3.5$, and

$$\begin{eqnarray}&&{\mathrm{BC}}_{K}=(6.8-0.2[\mathrm{Fe}/{\rm{H}}])(3.96-\mathrm{log}{T}_{{\rm{eff}}})\end{eqnarray} \tag{ 5 }$$

for RG and RC stars. This correction yields the bolometric magnitude of the star: ${m}_{\mathrm{bol}}={K}_{0}+{\mathrm{BC}}_{K}$.

For stars selected as dwarf and subgiant stars, we adopted the Torres et al. (2010) relations to estimate the mass and radius of the primary star:

$$\begin{eqnarray}\mathrm{log}{M}_{\star } & = & {a}_{1}+{a}_{2}X+{a}_{3}{X}^{2}+{a}_{4}{X}^{3}\\ & & +{a}_{5}{(\mathrm{log}g)}^{2}+{a}_{6}{(\mathrm{log}g)}^{3}+{a}_{7}[\mathrm{Fe}/{\rm{H}}],\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}\mathrm{log}{R}_{\star } & = & {b}_{1}+{b}_{2}X+{b}_{3}{X}^{2}+{b}_{4}{X}^{3}\\ & & +{b}_{5}{(\mathrm{log}g)}^{2}+{b}_{6}{(\mathrm{log}g)}^{3}+{b}_{7}[\mathrm{Fe}/{\rm{H}}],\end{eqnarray} \tag{ 7 }$$

where $X=\mathrm{log}{T}_{{\rm{eff}}}-4.1$ and the coefficients, a_i and b_i are given in Table 4 of Torres et al. (2010). This empirical relationship has a scatter of $6.4\%$ in mass and $3.2\%$ in radius, so for dwarfs and subgiants, we adopt ${\sigma }_{{M}_{\star }}=0.064{M}_{\star }$ as the uncertainty in the mass, and ${\sigma }_{{R}_{\star }}=0.032{R}_{\star }$. This information allows one to estimate the luminosity, L_⋆, as well as the distance, d, to these stars:

$$\begin{eqnarray}&&{L}_{\star }=4\pi {R}_{\star }^{2}{\sigma }_{\mathrm{SB}}{T}_{{\rm{eff}}}^{4}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{M}_{\mathrm{bol}}=4.77-2.5\mathrm{log}\left(\displaystyle \frac{{L}_{\star }}{{L}_{\odot }}\right)\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&d={10}^{1+0.2({m}_{\mathrm{bol}}-{M}_{\mathrm{bol}})},\end{eqnarray} \tag{ 10 }$$

where M_bol is the star's absolute bolometric magnitude. Uncertainty for these parameters are also derived through normal propagation of uncertainties, which yields a $13.5\%$ typical distance uncertainty for dwarfs and subgiants. A total of 340 of the 907 stars for which fitting was attempted used this prescription.

Unfortunately, the Torres et al. (2010) relations are not applicable to giant and PMS stars. For example, using the Torres et al. (2010) relations to derive the mass of Arcturus (${T}_{{\rm{eff}}}=4286$ K, $\mathrm{log}g=1.66$, [Fe/H] = −0.52) yields a mass of $3.5\;{M}_{\odot }$ compared to the accepted mass of $1.08\;{M}_{\odot }$ (Ramírez & Allende Prieto 2011). Therefore, we must resort to alternate methods for estimating the mass of the primary.

Efforts are currently underway to compile all published (or soon-to-be published) distance measurements to APOGEE stars. For stars selected as RG and RC stars, we employ a preliminary version of this distance catalog as the basis for our mass derivation. The most accurate distances for APOGEE stars are those derived from asteroseismic parameters from the APOGEE-Kepler catalog (APOKASC; Pinsonneault et al. 2014). These distances were given first priority because they only have $\sim 2\%$ random errors (Rodrigues et al. 2014). Unfortunately, no stars in this sample matched APOKASC stars with distance measurements, but we include it in the pipeline in hopes that future versions of the APOKASC catalog will overlap with future versions of this catalog. Our second choice, if the star is a RC star, is to use distances derived from the APOGEE RC catalog. These distances are cited to have 5%–10% random errors, and 71 stars of the 907 run through the apOrbit pipeline are RC stars. If the star has neither of the above distances available, we adopt the spectrophotometic distance estimates derived by Santiago et al. (2015), Hayden et al. (2015), or Schultheis et al. (2014), based on which estimate has the lowest error. These distances generally have <15%–20% uncertainties, and for most of the RG stars run through the apOrbit pipeline (489 stars), we adopt these distances. The six PMS stars in this sample are located in the young cluster IC 348 ($d=316\pm 22$ pc; Herbig 1998), so we adopt the distance to this cluster as the approximate distance to these stars. From the adopted distance, d, we estimate the luminosity of the star, and thus its radius and mass:

$$\begin{eqnarray}&&{M}_{\mathrm{bol}}={m}_{\mathrm{bol}}-5\mathrm{log}(d)+5\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{L}_{\star }={10}^{-0.4({M}_{\mathrm{bol}}-4.77)}{L}_{\odot }\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{R}_{\star }=\sqrt{\displaystyle \frac{{L}_{\star }}{4\pi {\sigma }_{\mathrm{SB}}{T}_{{\rm{eff}}}^{4}}}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{M}_{\star }=\displaystyle \frac{{10}^{\mathrm{log}g}{R}_{\star }^{2}}{G}.\end{eqnarray} \tag{ 14 }$$

Following typical propagation of uncertainties, these techniques produce a mass uncertainty floor of 26% due to the uncertainty in $\mathrm{log}g$. The median of mass uncertainties for these techniques is around 28%.

If a giant star has no distance measurement available, we adopt a characteristic mass from a TRILEGAL (Girardi et al. 2005) simulation using parameters typical of APOGEE giants. The median mass for all stars in this simulation with $\mathrm{log}g\lt 3.8$ and 3500 K $\lt \quad {T}_{{\rm{eff}}}\quad \lt$ 5000 K in the direction of Galactic Coordinates $({\ell },b)=(0,40)$ is ${M}_{\star }=1.6\pm 0.6\;{M}_{\odot }$ (∼40% mass uncertainty), which we adopt as the typical mass for all giant stars without a distance measurement. From this we derive ${R}_{\star }={({{GM}}_{\star }/{10}^{\mathrm{log}g})}^{1/2}$ and d, as for the dwarfs, both with typical estimated uncertainties of 25%. Fortunately, we only need to adopt this type of mass estimate for one star run through the apOrbit pipeline.

We employ the Fast ${\chi }^{2}$ Period Search ($F{\chi }^{2}$) algorithm (Palmer 2009) to search for periodic signals. This algorithm chooses the period based on the largest reduction in ${\chi }^{2}$ between a sinusoidal fit employing the first n_h harmonics of a fundamental period, p_i, compared to a global n_d-degree polynomal fit. The $F{\chi }^{2}$ algorithm uses harmonics of the fundamental period in its fits, which produces improved performance with non-circular orbits compared to the traditional Lomb–Scargle algorithm (Scargle 1982). Another advantage of the $F{\chi }^{2}$ algorithm is a built-in avoidance of periodic signals introduced by the cadence of the data, i.e., inputting data taken every n days will not return a n-day period as the best fit.

For our purposes, we employ three harmonics (n_H = 3), execute a search in four (logarithmic) period bins (0.3–3 days, 3–30 days, 30–300 days, and 300–3000 days), and oversample ten times the default frequency sampling such that the frequency step is ${\rm{\Delta }}f=1/(10{n}_{h}{\rm{\Delta }}T)$, where ${\rm{\Delta }}T$ is the longest temporal baseline of the observations. The search is executed once with a constant (n_d = 0) fit and once with a linear fit (n_d = 1). The periods in each bin, p_j that produce the greatest reduction in ${\chi }^{2}$, ${\rm{\Delta }}{\chi }_{\mathrm{max}}^{2}$, are then assessed for their significance using the following criterion:

$$\begin{eqnarray}&&{P}_{n-2{n}_{h}}({\rm{\Delta }}{\chi }_{\mathrm{max}}^{2})\geqslant 0.997,\end{eqnarray} \tag{ 15 }$$

where ${P}_{n-2{n}_{h}}({\rm{\Delta }}{\chi }_{\mathrm{max}}^{2})$ is the probability for a ${\chi }^{2}$ distribution with $n-2{n}_{h}$ degrees of freedom, and n is the number of RV epochs. The above limit is the equivalent of a $3\sigma$ detection. Periods that are not deemed significant by this metric are not used for full Keplerian orbit fitting. The significant periods (p_j) and their harmonics (1/3, 1/2, 2, and 3 times each value of p_j) are then each used for Keplerian orbit fitting.

Once the best periods are identified, Keplerian models with those periods are fit to the RV measurements using the MPFIT algorithm (Markwardt 2009). MPFIT is a Levenberg–Marquardt nonlinear least squares fitter implemented in IDL. This code is wrapped in an IDL code MP_RVFIT used in the MARVELS survey (De Lee et al. 2013). MP_RVFIT takes the input period and searches parameter space of the other Keplerian orbital parameters ($K,e,{\rm{\Omega }},{T}_{p}$, and global velocity trends) and returns the Keplerian model that satisfies the period with the lowest ${\chi }^{2}$.

Having a precise period is extremely important for acquiring an accurate Keplerian model, and simply submitting the periods from the period-finding algorithm to MP_RVFIT often leads to unsatisfactory results. Here we describe the bisector method implemented to achieve the best possible period. We initially submit the periods described above to MP_RVFIT, and keep the three periods (${p}_{k,0}$) that produce the best fits based on the modified reduced-chi-squared goodness of fit statistic, ${\chi }_{\mathrm{mod}}^{2}$, described in Section 3.3.4. For each of these periods we implement a bisector method to narrow in on the exact period. For each p_k we run MP_RVFIT with three periods: ${p}_{k,0}$ and ${p}_{k,0}\pm {\rm{\Delta }}{p}_{0}$, where ${\rm{\Delta }}{p}_{0}=0.5{p}_{k,0}$. We then compare the ${\chi }_{\mathrm{mod}}^{2}$ for the best fits for the three periods, and update p_k and ${\rm{\Delta }}p$ accordingly:

$$\begin{eqnarray}&&\mathrm{If}\ {\chi }_{{p}_{k,i}}^{2}\leqslant {\chi }_{{p}_{k,i}\pm {\rm{\Delta }}{p}_{i}}^{2}:\\ &&\qquad {p}_{k,i+1}={p}_{k,i},{\rm{\Delta }}{p}_{i+1}={\rm{\Delta }}{p}_{i}/2,\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\mathrm{If}\ {\chi }_{{p}_{k,i}\pm {\rm{\Delta }}{p}_{i}}^{2}\lt {\chi }_{{p}_{k,i}}^{2}:\\ &&\qquad {p}_{k,i+1}={p}_{k,i}\pm {\rm{\Delta }}{p}_{i},{\rm{\Delta }}{p}_{i+1}={\rm{\Delta }}{p}_{i}.\end{eqnarray} \tag{ 17 }$$

For the ${\chi }_{{p}_{i}-{\rm{\Delta }}{p}_{i}}^{2}\lt {\chi }_{{p}_{i}}^{2}$ case, if ${p}_{i}-2{\rm{\Delta }}{p}_{i}\lt 0.1$, then we use ${\rm{\Delta }}{p}_{i}={\rm{\Delta }}{p}_{i}/2$ for the next update. This iteration is performed until the change in ${\chi }_{\mathrm{mod}}^{2}$ is less than 0.01 or ${n}_{\mathrm{iter}}=50$ iterations are reached. The distribution of the required number of iterations for systems in the final sample had a median of 15 with few systems above 25. Therefore, the choice to terminate systems on their 50th iteration is more than justified as these systems are unlikely to converge in a timely manner. These systems are also not included in the final catalog (see Section 4.2). The final values of ${p}_{k,{n}_{\mathrm{iter}}}$ are then submitted to MP_RVFIT one final time, and the results saved for the catalog. The data saved are described in Section 4. For a few example Keplerian orbit models see Figure 5.

Zoom In Zoom Out Reset image size

Directly from the orbital parameters, we can calculate the projected semimajor axis of the primary star:

$$\begin{eqnarray}&&{a}_{\star }\mathrm{sin}i=\displaystyle \frac{\mathrm{KP}}{2\pi }\sqrt{1-{e}^{2}}.\end{eqnarray} \tag{ 18 }$$

From this measurement we can define the mass function of the system:

$$\begin{eqnarray}&&f(m,{M}_{\star })=4{\pi }^{2}\displaystyle \frac{{({a}_{\star }\mathrm{sin}i)}^{3}}{{\mathrm{GP}}^{2}}=\displaystyle \frac{{(m\mathrm{sin}i)}^{3}}{{({M}_{\star }+m)}^{2}}.\end{eqnarray} \tag{ 19 }$$

This quantity is saved in the catalog, but we also attempt to estimate the secondary mass directly:

$$\begin{eqnarray}&&m\mathrm{sin}i={[f(m,{M}_{\star }){M}_{\star }^{2}{(1+(m/{M}_{\star }))}^{2}]}^{1/3}.\end{eqnarray} \tag{ 20 }$$

The general case of this equation cannot be solved analytically, but often when dealing with planetary companions, we can make the assumption that $m\ll {M}_{\star }$, and thus can make the approximation $m\mathrm{sin}i\approx {(f(m,{M}_{\star }){M}_{\star }^{2})}^{1/3}$. For companions with $m\mathrm{sin}i\lt 0.1{M}_{\star }$, this approximation is accurate to within 10%, but this sample contains higher-mass companions for which we want reasonable mass estimates. In these cases we solve the above equation iteratively, initially assuming m = 0, returning the above estimate, and iterating until $m\mathrm{sin}i$ changes by $\lt {10}^{-4}\;{M}_{\odot }$. Since we are interested in estimating the minimum mass of the companion, we solve for the $\mathrm{sin}i=1$ case, and thus use $m\approx m\mathrm{sin}i$ after the first iteration. This iterative method for determining m was tested for a variety of mass ratios and a variety of starting points for m (not just m = 0). From these tests, we have found this method to be quite robust.

Finally, from the estimate of $m\mathrm{sin}i$, we provide an estimate of the semimajor axis, a, of the secondary:

$$\begin{eqnarray}&&a={a}_{\star }\mathrm{sin}i\displaystyle \frac{{M}_{\star }}{m\mathrm{sin}i}.\end{eqnarray} \tag{ 21 }$$

Finally we compile the three best models from the run with no global linear fit and the three best models from the linear fit run, and compare them to select the best overall fit. Ideally the phase and velocity coverage of the model are uniformly sampled by the data, and we aimed to preferably select models that are as close to this ideal as possible. A useful way to quantify the phase coverage of the data is the uniformity index (Madore & Freedman 2005):

$$\begin{eqnarray}&&{U}_{N}=\displaystyle \frac{N}{N-1}\left[1-\displaystyle \sum _{i=1}^{N}{({\phi }_{i+1}-{\phi }_{i})}^{2},\right]\end{eqnarray} \tag{ 22 }$$

where the values ${\phi }_{i}$ are the sorted phases associated with the corresponding Modified Julian Date (MJD) of the measurement i, and ${\phi }_{N+1}={\phi }_{1}+1$. This statistic is normalized such that $0\leqslant {U}_{N}\leqslant 1$, where U_N = 1 would indicate a curve evenly sampled in phase space. Using a similar derivation, we also define an analogous "velocity" uniformity index with the same properties as U_N:

$$\begin{eqnarray}&&{V}_{N}=\displaystyle \frac{N}{N-1}\left[1-\displaystyle \sum _{i=1}^{N}{({\nu }_{i+1}-{\nu }_{i})}^{2}\right].\end{eqnarray} \tag{ 23 }$$

We define a "velocity phase," ${\nu }_{i}=({v}_{i}-{v}_{\mathrm{min}})/({v}_{\mathrm{max}}-{v}_{\mathrm{min}})$, to have the same properties as ${\phi }_{i}$ above, where the values $\nu =0$ and $\nu =1$ indicate the minimum and maximum velocities of the model, v_min and v_max. The values of v_i are the RV measurements, sorted by their value, with the adopted global velocity trend subtracted. For models that do not apply a global linear trend, the trend subtracted is the average of the raw velocities: ${v}_{i}={v}_{\mathrm{raw},i}-{\bar{v}}_{\mathrm{raw}}$. Measured velocities below the minimum or above the maximum are assigned $\nu =0$ and $\nu =1$, respectively. The purpose of this metric is to prevent the pipeline from selecting an extremely eccentric orbit when the data do not support such a model. Values of U_N and V_N are given in the example RV curves of Figure 5.

Combining the above statistic with the traditional reduced ${\chi }^{2}$ goodness-of-fit statistic (${\chi }_{\mathrm{red}}^{2}$), we define the modified ${\chi }^{2}$ statistic,

$$\begin{eqnarray}&&{\chi }_{\mathrm{mod}}^{2}=\displaystyle \frac{{\chi }_{\mathrm{red}}^{2}}{\sqrt{{U}_{N}{V}_{N}}},\end{eqnarray} \tag{ 24 }$$

by which the models are ranked. In the case that U_N = 0 or V_N = 0, ${\chi }_{\mathrm{mod}}^{2}$ would be recorded as a floating-point infinity and automatically be ranked below all other fits. However, there are some conditions where the fit is unacceptable, but still may be selected as the best fit using the above metric. Therefore, we defined criteria that split the fits into "good" and "marginal" fits. Any of the following criteria would warrant a "marginal" classification:

• 1.  
  Periods within 5% of 3, 2, 1, 1/2, or 1/3 day,
• 2.  
  Periods, P, longer than twice the baseline, $2{\rm{\Delta }}T$,
• 3.  
  Extremely eccentric solutions ($e\gt 0.934$)²⁶ ,
• 4.  
  Orbital solutions that send the companion into the host star: $a(1-e)\lt {R}_{\star }$,
• 5.  
  Poor phase and velocity coverage (${U}_{N}{V}_{N}\lt 0.5$).

The good and marginal fits are ranked by ${\chi }_{\mathrm{mod}}^{2}$ separately, and the best fit is the good fit with the lowest ${\chi }_{\mathrm{mod}}^{2}$. If all of the fits were deemed marginal, then the best fit is the marginal fit with the lowest ${\chi }_{\mathrm{mod}}^{2}.$ For more details on the verification and performance of the apOrbit pipeline, see Appendix A.

The top panel of Figure 6 indicates that this sample reproduces the BD desert, but only for orbits with $a\quad \lt$ 0.1–0.2 AU ($P\quad \lt$ 10–30 days), which is significantly less than the 3 AU extent of the desert as stated in Grether & Lineweaver (2006). However, their sample mostly considered solar-like dwarf hosts, while this sample contains stars with a variety of spectral types, as well as many evolved stars. From the top panel of Figure 7, it appears that the relative number of BD companions decreases as host mass increases for MS hosts. Likely M dwarfs (MS with ${M}_{\star }\lt 0.6\;{M}_{\odot }$) have roughly equal numbers of BD and stellar-mass companions, while K dwarfs (MS with $0.6\lt {M}_{\star }/{M}_{\odot }\lt 0.85$) have roughly half the number of BD candidate companions as stellar-mass candidate companions. The G dwarfs (MS with $0.85\lt {M}_{\star }/{M}_{\odot }\lt 1.1$) show a similar relative number of BD companions compared to stellar-mass companions, but they are less uniformly distributed throughout the BD mass regime compared to the lower mass BD candidate hosts, suggesting a higher probability that many of these BD candidates are scattered into the BD mass regime by inclination effects. These results leads one to believe the interpretation of Duchêne & Kraus (2013) that the BD desert is simply a special case for solar-mass stars of a more general lack of extreme mass ratio ($q\lesssim 0.1$) systems. For example, if, in general, systems with $q\lt 0.08$ are rare (i.e., a BD companion around a 1 M_⊙ companion), then a relatively high-mass BD companion ($m\gt 0.04\;{M}_{\odot }$) orbiting a $0.5\;{M}_{\odot }$ star should be a more common occurrence.

Zoom In Zoom Out Reset image size

Out of the 112 BD companion candidates in this sample, 71 orbit evolved stars. All but two of the giant (RC and RG) hosts have masses $\gt 0.8\;{M}_{\odot }$ and only one of the SG hosts has a mass $\lt 1\;{M}_{\odot }$. Considering that stars like the Sun lose up to a third of their mass on the red giant branch (RGB), it is a reasonable assumption that a vast majority of the evolved stars in this sample descended from main-sequence F (or earlier) dwarfs. As can be seen from the bottom panel of Figure 7, the evolved stars have roughly half the number of BD candidate companions as stellar-mass candidate companions, and the BD-mass candidates are distributed throughout the BD-mass regime, similar to the K dwarf distribution. If the evolved stars are indeed evolved F dwarfs, and we follow the progression from above, one would expect these stars to have a smaller relative number of BD companions compared to even the G dwarfs. However, it has been previously suggested that the BD desert observed for Solar-like stars may cease to exist for F dwarf stars (Guillot et al. 2014). Their proposed explanation of this effect is that G dwarfs are more efficient at tidal dissipation. In general, compared to Jupiter-mass planets, more massive small separation companions undergo stronger tidal interaction with their host star through angular momentum exchange. Stellar-mass companions, however, have sufficient orbital angular momentum to remain in a stable orbit, which explains the demise of small separation BD-mass but not stellar-mass companions. However, F (and earlier) dwarfs are known to remain rapid rotators (${v}_{\mathrm{rot}}\;\sim$ 20–100 km s⁻¹) throughout their main-sequence lifetimes due to their smaller outer convective zones leading to weaker magnetic breaking. This means F dwarfs are also less efficient at extracting angular momentum from an orbiting companion. Therefore, rapid rotators such as F dwarfs inhibit tidal dissipation, which explains this "F dwarf oasis" for BD companions. The dynamical model presented in Figure 4 of Guillot et al. (2014) shows that a companion in the BD-mass regime on an initial 3 day orbit around a $1\;{M}_{\odot }$ star will survive for $\lt 40\%$ of the star's main sequence lifetime ($\lesssim 4$ Gyr), while the same companion around a host star with ${M}_{\star }\gt 1.2\;{M}_{\odot }$ will survive for at least the entirety of the host star's main sequence lifetime (∼6.5 Gyr for a $1.2\;{M}_{\odot }$ star). The presence of a large number BD companions orbiting the evolved stars in this sample strongly supports this "F dwarf oasis" hypothesis.

However, the tidal effects explanation would only strongly affect the closest-in companions. Since the rotation period of a G dwarf is ${P}_{\star }=30$ days (compared to a few days for an F dwarf), tidal dissipation could only explain BD companions with orbital periods less than 30 days, and the majority of the BD candidate companions in this sample have periods significantly greater than that. Therefore, tidal dissipation can only explain the BDs (or lack thereof) with orbits within 0.2 AU. Curiously, this sample reproduces the BD desert out to approximately 0.2 AU, suggesting this mechanism may indeed play a role in shaping the BD desert. Another possible explanation for the presence of BD candidate companions is Roche lobe overflow of the star as it evolves off the main sequence onto an orbiting planetary-mass candidate, allowing it to grow to BD mass as the star evolves up the RGB. Eggleton (1983) gives the following approximation for the Roche lobe of a primary donor star with mass M₁ orbited by a companion with M₂:

$$\begin{eqnarray}&&\displaystyle \frac{{r}_{1}}{a}=\displaystyle \frac{0.49{q}^{-2/3}}{0.6{q}^{-2/3}+\mathrm{ln}(1+{q}^{-1/3})},\end{eqnarray} \tag{ 31 }$$

where, $q={M}_{2}/{M}_{1}$, a is the separation of the two bodies, and r₁ is the Roche lobe radius of the potential donor. In the bottom panel of Figure 6, we mark the Roche lobe as a function of $\langle q\rangle$. As an interesting note, it appears there are seven stars in this sample that are currently at or near Roche lobe overflow, three of which are currently of planetary mass, and three of which are BD mass. These systems will all be the subject of further scrutiny. In general, for a 1–10 M_Jup planet to cause a $\sim 1\;{M}_{\odot }$ primary to overflow its Roche lobe, the radius of the primary would have to exceed ∼70%–80% of the separation between the two bodies. This would not be an unreasonable expectation for a companion orbiting within 1 AU, as solar-mass stars can achieve radii approaching 1 AU at the tip of the RGB. Therefore, this mechanism may be a way to explain the relatively large number of BD companion candidates orbiting the evolved stars in this sample. Overall, this catalog's large number of systems with short-period BD companion candidates challenges the notion of the BD desert as we know it, and certainly warrants further investigation.

In Figure 6, we also see the distribution of orbital eccentricities. As expected, the smallest-separation ($a\lt 0.1$ AU) stellar-mass companions all have circular orbits. The circularization cutoff period increases with the age of the system with 5–10 Gyr systems having cutoff periods of 12–20 days (Mathieu et al. 2004). All of the binary companions in this catalog with $a\lt 0.1\;{\rm{AU}}$ have $P\lt 20$ days. Therefore, the distribution of eccentricities for the binary systems in this sample, with circular orbits at small separation, and eccentric orbits a large separations is not unexpected. The closest ($a\lt 0.01$ AU) planetary-mass companions appear to have also circularized, as expected, but a surprising result is the relatively large fraction of eccentric orbits for relatively close-in planetary-mass candidate companions. For the RG and RC hosts, one interpretation of these eccentricities is ongoing tidally induced migration (see Section 5.2.1 for further discussion of this). However, the majority of the small-separation planetary and BD candidate companions orbit dwarf and SG stars. For these systems, their higher eccentricities may be further evidence for the mechanism suggested by Tsang et al. (2014) whereby stellar illumination heating a gap cleared by a forming planet may excite the eccentricity of the planet in the gap.

Of this catalog's candidate companion systems, there are 50 systems with a mass ratio $\langle q\rangle =\langle m\rangle /{M}_{\star }\geqslant 0.5$ (see bottom panel of Figure 6). One would expect that these systems would manifest themselves as SB2s, but these systems show no strong indication of such behavior in their APOGEE spectra. Of these, 24 are RG stars, which would explain their lack of SB2 behavior, as their companion is likely still on the main sequence, and thus the flux ratio would be too large. However, this still leaves 26 MS and SG hosts, of which one explanation is that they host massive compact objects, such as stellar remnants. These seven companions have $0.3\;{M}_{\odot }\lt \langle m\rangle \lt 1.2\;{M}_{\odot }$, which would indicate these systems might host white dwarf companions, eight of which may be low-mass ($\langle m\rangle \lt 0.45\;{M}_{\odot }$) He-core white dwarfs (Liebert et al. 2005). Furthermore, two of the systems with a RG host have $\langle q\rangle \gt 1$, indicating the companion has already completed its evolution, and the recovered $\langle m\rangle$ of the companions (2.8 and $1.6\;{M}_{\odot }$) indicates they may be neutron stars.

Tidal dissipation is thought to play an important role in the destruction of planetary systems as a star evolves off the main sequence and expands (Penev et al. 2012). This sample contains 225 $a\lt 3\;{\rm{AU}}$ candidate companions to evolved stars, indicating either many initial small separation companions survive engulfment or farther-orbiting planets undergo increasing tidal migration as its host ascends the giant branch, bringing the companion closer to its host star. The nine candidate planetary-mass ($\langle m\rangle \lt 0.013\;{M}_{\odot }$) companions orbiting giant stars in this sample would be a 20% increase in the number of currently known giant stars hosting a planet (∼50 according to the tabulation by Jones et al. 2014a). As Jones et al. (2014a) mentions, there is a small separation cut-off for RG hosts. The current record-holder for smallest separation of an RV-detected planet RG host is HIP 67851b with $a=0.539\;{\rm{AU}}$ (Jones et al. 2014b). The shortest period planet orbiting a giant star, Kepler 91b, is on a 6 day orbit (Lillo-Box et al. 2014). Most of the candidate planets orbiting giants lie between these two systems, with a few candidates closer than Kepler 91b.

Of the evolved stars in this sample, 23 are verified RC stars (Bovy et al. 2014). RC stars are metal-rich stars which have passed through the tip of the RGB and have contracted due to the ignition of core helium burning. It is expected that stars like the Sun may reach radii up to 1 AU when they reach the tip of the RGB. Therefore, the presence of companion candidates orbiting RC stars at $a\lt 1\;\;\mathrm{AU}$ in this catalog (see Figure 9) is a surprising discovery. To investigate this further we compared the RC stars to RGs in this sample, but we consider only RGs with [Fe/H] > −0.42 ([Fe/H] of the most metal-poor RC in this sample) and $2.4\leqslant \mathrm{log}g\leqslant 3.3$ (the $\mathrm{log}g$ range of the RC stars) to eliminate possible effects from RV sensitivity issues. This also allows us to compare stars approximately half way up the giant branch to stars that have already passed through the tip of the RGB, and have achieved their largest extent. These 92 RG stars have 62 stellar-mass, 25 BD-mass, and 5 planet-mass companion candidates compared to 18, 5 and 0 for the 23 RC giants.

Zoom In Zoom Out Reset image size

A cursory look at these numbers (and Figure 9) shows a lack of smaller companions for RC stars, as well as a companion candidates found at smaller separations for the 92 RG stars when compared to RC stars (0.07 AU versus 0.2 AU at the low-mass end). It is also interesting to note that no companion candidates have circular orbits among RC hosts (smallest e = 0.284). This all points to the role of the tidal migration and destruction of companions, particularly planetary-mass companions. However, any tidally induced migration of companions will be much weaker than when the star was in the RGB phase. Therefore any companions with $a\lt 1\;{\rm{AU}}$ around an RC star likely would have to have survived inside the star's envelope during its RGB phase. Most of the RC hosts with $a\lt 1\;{\rm{AU}}$ candidate companions are likely post-common envelope systems, and thus may have experienced drag-induced migration to bring them to their current orbit. These systems certainly warrant deeper investigation.

According to the compilation of exoplanets.org (Han et al. 2014), of the confirmed planet-hosting stars with metallicity measurements, only 15 have [Fe/H] < −0.5. This sample has 41 stars with [Fe/H] < −0.5, and of these, two host candidate planetary-mass companions and 14 host candidate BD companions. The most metal-poor stars in this sample approach [Fe/H] = −2. While there are no candidate planets among the most metal-poor ([Fe/H] < −1) hosts (5 stars), there are two companions in the BD mass regime. The smaller fraction of the lowest-mass companions detected among the most metal-poor stars in this sample is not surprising as the RV uncertainties are higher for metal-poor stars, as described in Equation (1). Also, it is not too surprising to find metal-poor stars hosting binary companions, given the Carney et al. (2003) result. However, finding a population of metal-poor stars potentially hosting BD companions is surprising in the context of the core accretion model of companion formation, and may suggest an alternate formation mechanism for these companions.

We simulated 9000 planetary systems with random characteristics. The masses of the primary stars were drawn from the distribution of estimated masses for the actual candidate substellar hosts in the APOGEE data. The companion masses and periods were drawn from the distributions specified by Tabachnik & Tremaine (2002), using a mass range of 1–100 M_Jup and periods from 0.1 to 2000 days. Eccentricities were drawn from a uniform distribution with a maximum of e = 0.934, which corresponds to the eccentricity of HD 80606b, the largest eccentricity in the exoplanets.org. database. Companions with $P\lt 5$ days were assumed to have circular orbits. The radii of the APOGEE candidate host stars were also estimated, and planets with orbital separations less than 5 ${R}_{\star }$ were considered unphysical because the tidal decay of planetary orbits becomes relevant at such small separations. The longitude of periastron and the orbital phase of periastron passage relative to a reference date were drawn from uniform distributions.

With the orbital characteristics of the simulated companions defined, we simply used the helio_rv code in the IDL astronomy users library²⁷ to calculate the measured heliocentric RV for each system on a set of observation dates. The observation dates for each system were designed to mimic the way the survey proceeded. The observations for each star were spread randomly over a 3.2 year time period assuming the telescope was on-sky for 15 days followed by 14 days off sky since APOGEE observed primarily during bright time. The simulated measured RVs consisted of the actual motion of the star at the time of observation plus two sources of noise, drawn from Gaussian distributions. The first is simply measurement noise, which nominally has ${\sigma }_{v}=100\;{\rm{m}}\;{{\rm{s}}}^{-1}$ but is increased to ${\sigma }_{v}=130\;{\rm{m}}$ s⁻¹ for 20% of the visits to simulate poor observing conditions. The second noise source is intrinsic stellar atmospheric RV jitter, with an amplitude drawn from the distribution in Frink et al. (2001).

A second data set of 9000 simulated system were generated with much of the same parameters as the first, except that it had mass ratios approaching one, all orbital parameters were drawn from a uniform distribution, and it was much sparser in the lower-mass companion regime. We combined these two data sets to obtain complete coverage of parameter space. From the combined data set, we generated RV curves with 9, 12, 16, and 24 visits selected from the 24-visit parent sample with the 100 m s⁻¹ RV uncertainty level, as well as a set where the base uncertainty level is inflated to 1 km s⁻¹ to emulate RV measurements from the metal-poor host stars in this sample.

The full test suite of simulated systems was run through the apOrbit pipeline each time an update to the fitting algorithms was implemented. Many of these updates were inspired by the simulated systems for which the pipeline failed to reproduce the correct orbit in the previous run. In addition, many of the criteria used to select the RV variable and gold candidate sample, described in Section 4, were inspired by these results. Notable failures in previous runs that led to new selection criteria for candidate companions included:

• 1.  
  Long-Period Systems: the longest-period simulated systems demonstrated the largest scatter in their results. This inspired the use of the phase uniformity index (see Section 3.3.4), as well as a procedure to reject any solutions for which the period was longer than twice the baseline (Section 4.2).
• 2.  
  Highly Eccentric Systems: the code had the most difficulty reproducing the orbital parameters of systems with high eccentricity ($e\gt 0.9$). However, these systems are extremely rare, so this is not a major issue. Nevertheless, this result still led to the decision to reject all orbital solutions with $e\gt 0.934$ (Section 3.3.4), which is the planetary system with the largest known eccentricity anyway. Even with this cut, however, the more eccentric the system, the more trouble the code had in recovering the correct orbital parameters. In particular, systems with low numbers of visits had the most issues. This result inspired the use of the velocity uniformity index in the fitting code (Section 3.3.4), and it led to the decision to implement more stringent significance cuts for eccentric systems (Section 4.1).
• 3.  
  One-day Aliased Systems: early tests of the code on simulated systems revealed a tendency for solutions to cluster around integer fractions of one day, despite the initial period selection avoidance of such periods. This inspired the decision to reject any periods within 5% of 1/3, 1/2,1, 2, or 3 days (Section 3.3.4).

These simulations were also used to understand how RV noise from the star or measurement error can potentially lead to a false positive candidate companion. To accomplish this we ran the simulated systems through the apOrbit pipeline following the procedures laid out in Section 3.3 using just the RV signals from the star's atmospheric jitter and random measurement errors. We ran the simulations using four different numbers of visits (9, 12, 16, and 24) and two uncertainty levels (0.1 and 1 km s⁻¹), and selected candidates using the criteria described in Section 4. The results of the false positive tests are summarized in Figure 11. For systems with 24 visits, out of 18,000 simulated systems only two systems at the 1 km s⁻¹ uncertainty level registered as false positives. The raw number of false positives increased dramatically from 24 to 16 visits, and the higher uncertainties lead to a higher rate of false positives. Most false positives clustered around the sensitivity limit corresponding to the RV measurements from which they were derived (see Equation (32) below), and are generally assigned eccentric orbits ($e\gt 0.5$). This information, along with visual inspection of these fits, led to a pre-cut based on the velocity variations of the star which was based on the value of the statistic described in Section 3.1.2 for these systems.

Zoom In Zoom Out Reset image size

In the results presented here, we ran the 18,000 simulated systems through the apOrbit pipeline as described in Section 3.3.2 for the four visit levels (${n}_{\mathrm{RV}}=9,12,16,24$) and two uncertainty levels (${\sigma }_{v}=0.1,1$ km s⁻¹), and selected companion candidates in the same manner as the gold sample, as described in Sections 3.1.2, 3.3.4, and 4. The systems correctly recovered (P and K recovered within 10%, and e recovered within 0.1) are shown in Figure 12. Systems with large K but small errors can lead to larger values of ${\chi }^{2}$ for a fit that still produces the correct orbital parameters. Unfortunately, removing the ${\chi }_{\mathrm{mod}}^{2}$ constraint allows many systems with incorrect solutions to pass, so we err on the side of caution and keep it in place. Due to limits of the period search and the other constraints on the period, a, and ${\chi }_{\mathrm{mod}}^{2}$ of the fit described in Section 4.2, we expect to be able to recover orbits for companions having $0.01\;\mathrm{AU}\lesssim a\lesssim 3\;\mathrm{AU}$, depending on the baseline, number of visits, and the RV uncertainty level (see Figure 12). By fitting a trendline to the simulated systems with $2.8\leqslant K/{\tilde{\sigma }}_{v}\leqslant 3$ for various values of ${\sigma }_{v}$, we also find that the lower limit on detectability, which we will refer to as the sensitivity function (SF), of $m\mathrm{sin}i$ can be written as:

$$\begin{eqnarray}&&\mathrm{log}(m\mathrm{sin}i)=0.48\mathrm{log}(a)-C,\end{eqnarray} \tag{ 32 }$$

where the constant offset, C, depends on the sensitivity level (which we interpret as the median RV uncertainty, ${\tilde{\sigma }}_{v}$ ):

$$\begin{eqnarray}&&C/\mathrm{log}({M}_{\odot })=2.0-0.3{\mathrm{log}}_{2}({\tilde{\sigma }}_{v}/100\;{\rm{m}}\;{{\rm{s}}}^{-1}).\end{eqnarray} \tag{ 33 }$$

Zoom In Zoom Out Reset image size

We then compared each observed/recovered orbital parameter, X_o, with the actual parameters for the system, X, to determine how accurately the parameters are recovered as a function of parameter space. These results are summarized in Figure 13. For a large portion of the parameter space, the selected candidates reproduce the correct orbital parameters quite well. The systems that give the most trouble appear to be the low-mass companions, companions at large separations, and companions with highly eccentric orbits. Unsurprisingly, parameter recovery is overall worse for stars with fewer visits and higher RV uncertainties, but the drop in performance was not as dramatic between the 24 visit and 16 visit simulations as it was between the 16 visit to 9 visit simulations. However, from these results we can still conclude that in almost all regimes, the recovered orbital parameters are at least characteristic of the true values for the system.

Zoom In Zoom Out Reset image size

Finally, we construct the recovery rates across the parameter space covered by this catalog. These are summarized in Figure 14. Unsurprisingly, recovery rate drops as n_RV decreases, and higher RV uncertainties lead to lower recovery rates in general.

Zoom In Zoom Out Reset image size