Building and Calibrating the Binary Star Population Using Kepler Data

Galaxia takes a variety of parameters that fine tune the generated population, as shown in Table 1. Of all these parameters, we will only focus on the magnitude limits and sky position as these are unique to our process. We generate stars down to a magnitude of r_SDSS = 25. While this is roughly 5 orders of magnitude fainter than the magnitude limit of Kepler, it is required in order to generate the low-mass, and subsequently faint, stellar components that will become secondaries. We create circular regions that are centered on, and circumscribe, each of the 22 Kepler modules. These synthetic stellar populations are then trimmed using the Python package K2fov ⁴ (Mullally et al. 2016) to identify which objects actually fall on the active silicon of the detector. The 22 stellar populations, now trimmed, are combined and constitute our population of single stars (see Figure 1).

Zoom In Zoom Out Reset image size

The first step in synthesizing the binaries is to select a primary star for each system. Primaries are selected from the single star population in accordance with their respective multiplicity fraction. The multiplicity fraction is defined as the ratio:

$$\begin{eqnarray}&&{f}_{m}\equiv \displaystyle \frac{B+T+\ldots }{S+B+T+\ldots },\end{eqnarray} \tag{ 1 }$$

where S, B, T, ... are the numbers of single, binary, triple, and higher-order multiple systems, respectively. For a star of mass M in the stellar population, we say it has a probability of f_m to be drawn as a primary star. This relationship depends on mass and we adopt the following analytical relationship from Arenou (2011):

$$\begin{eqnarray}&&{f}_{m}(M)={c}_{1}\tanh ({c}_{2}M+{c}_{3}),\end{eqnarray} \tag{ 2 }$$

where c₁, c₂, and c₃ are free parameters and M is the stellar mass. We modify the original coefficients given by Arenou (2011) to better fit data reported by Duchêne & Kraus (2013) and Raghavan et al. (2010) (see the top panel of Figure 2). The coefficients given by Arenou (2011) are c₁ = 0.8388, c₂ = 0.688, and c₃ = 0.079, while our modified values are c₁ = 1, c₂ = 0.31 ± 0.05, and c₃ = 0.18 ± 0.04. We have fixed c₁ to unity as only 6₋₃ ⁺⁶% ${6}_{-3}^{+6}{\rm{ \% }}$ of O-type stars are singles (Moe & Di Stefano 2017).

Zoom In Zoom Out Reset image size

Due to the inherent simplicity of the analytical model assumed here, our derived relationship potentially overestimates the multiplicity rates above 3 M_⊙. We expect this to be a marginal issue for the Kepler field as approximately only 400 single stars, out of 28.7 million, have masses above 3 M_⊙, In addition, even if all of these systems, single stars and binaries with primaries above 3 M_⊙, made it into the target list, this would only make up approximately 0.2% of the observed systems.

For higher-mass populations, a better analytical model will need to be considered. After the primary stars have been drawn, we turn our attention to the orbital parameters.

Binary star orbits are characterized by the mass ratio (q = M₂/M₁), orbital period (log P with P in days), orbital eccentricity (e), semimajor axis (a), inclination (i), and the argument of periastron (ω). It is necessary to place strict limits on mass ratio and period to ensure that only viable systems are generated. Mass ratio, along with the mass of the primary, determines the mass of the secondary. The minimum mass for a star to sustain nuclear fusion is taken to be 0.07 M_⊙, chosen to agree with Galaxia. This sets a limit on the mass ratio, as a function of primary mass, given by:

$$\begin{eqnarray}&&{q}_{\min }({M}_{1})=\displaystyle \frac{0.07\,{M}_{\odot }}{{M}_{1}}.\end{eqnarray} \tag{ 3 }$$

The relationship between mass ratio and primary mass is illustrated in Figure 3. The truncation of the distribution is clearly visible and is more pronounced for lower-mass stars. After the mass ratio has been drawn, we can draw secondary stars.

Zoom In Zoom Out Reset image size

In addition to the mass requirement set by mass ratio, we also require the stars to be coeval. To determine the most appropriate secondary star from the synthetic star population to serve as a secondary we use the following merit function for each binary

$$\begin{eqnarray}&&S(M,\mathrm{log}\,{A}_{})=\sqrt{{\left(\displaystyle \frac{M-{{qM}}_{1}}{{M}_{\max }-{M}_{\min }}\right)}^{2}+{\left(\displaystyle \frac{\mathrm{log}{A}_{}-\mathrm{log}{A}_{1}}{\mathrm{log}{A}_{\max }-\mathrm{log}{A}_{\min }}\right)}^{2}},\,\end{eqnarray} \tag{ 4 }$$

where S is our secondary star merit function and M_min, M_max, $\mathrm{log}\,{A}_{\min }$, and $\mathrm{log}\,{A}_{\max }$ are the minimum and maximum masses (in solar units) and ages (in gigayears) of the synthetic star sample, respectively. We divide by the difference of the extrema to place the parameters on a normalized scale. A secondary is selected if the value of S does not exceed 0.01 and if M < M₁. We include the condition M < M₁ to ensure that our primaries are the more massive members of the system. The value of 0.01 corresponds to a combined normalized parameter deviation of 1% from the desired mass of qM₁ and $\mathrm{log}\,{A}_{1}$. The minimum and maximum values of mass and age used are ${M}_{\min }=0.07\,{M}_{\odot }$, ${M}_{\max }\,=7.49\,{M}_{\odot }$, $\mathrm{log}\,{A}_{\min }=4.94$, and $\mathrm{log}\,{A}_{\max }=10.1$. Therefore the selected secondary is guaranteed to have $\left|{M}_{2}-{{qM}}_{1}\right|\leqslant 0.07\,{M}_{\odot }$ and $\left|\mathrm{log}\,{A}_{2}-\mathrm{log}\,{A}_{1}\right|\leqslant 0.05$.

The drawn secondary adopts the nonstellar properties of the primary such as location in the sky, distance, and interstellar extinction. This process does not preserve stellar mass content as secondaries are drawn with replacement. The resulting synthesized populated of single stars and binaries is inflated relative to the model stellar population produced with Galaxia.

Once secondaries have been drawn for each system, we can draw period $(\mathrm{log}P)$ and eccentricity (e). To model the $\mathrm{log}P$ distribution, we construct a binned model, with n bins of equal width and an additional overflow bin. The n bins span the observable eclipsing binary $\mathrm{log}P$ range of −0.64 (5.5 hr) to 2.86 (730 days). Of course binaries exist outside of this range and, to account for this, we have included an additional bin which ranges up to 8.0 (≈270 kyr). While the n bins that fall within the observable range are fit explicitly, the overflow bin is fit by requiring the integral of the model density to be unity. Hence, the overflow bin scales the rest of the model and is an estimate of how many binaries are above the maximum observable period.

Before we can draw from the period model, we must first determine the minimum allowed period, ${\rm{l}}{\rm{o}}{\rm{g}}{P}_{min}$. To ensure that we only generate detached binaries, we filter the separation of the system at periastron with:

$$\begin{eqnarray}&&a(1-e)\gt s\left({R}_{1}+{R}_{2}\right),\end{eqnarray} \tag{ 5 }$$

where a is the semimajor axis, s is the separation factor, R₁ is the radius of the primary, and R₂ is the radius of the secondary. The separation factor in Equation (5) must be set to a value above 1 to ensure that adequate separation between the two stars is maintained. A star that is critically rotating, such that it cannot spin faster without ejecting mass, will have its equatorial radius equal to $\tfrac{3}{2}$ of its polar radius (see Prša 2018 for a complete derivation). We have chosen to set $s=\tfrac{3}{2}$ to correspond to the situation in which we have both stars spinning at their critical breakup rotation. The minimum separation of the system:

$$\begin{eqnarray*}&&{a}_{\min }=\displaystyle \frac{3({R}_{1}+{R}_{2})}{2(1-e)}\end{eqnarray*}$$

is used, along with Kepler's third law, to obtain:

$$\begin{eqnarray*}&&{P}_{\min }=3\pi {\left(\displaystyle \frac{3}{2G({M}_{1}+{M}_{2})}\right)}^{\tfrac{1}{2}}{\left(\displaystyle \frac{{R}_{1}+{R}_{2}}{1-e}\right)}^{\tfrac{3}{2}}.\end{eqnarray*}$$

The period is drawn from the distribution with an additional requirement that ${\rm{l}}{\rm{o}}{\rm{g}}P\geqslant {\rm{l}}{\rm{o}}{\rm{g}}{P}_{min}$.

Eccentricity is uniformly drawn between 0 and ${e}_{\max }$ and occurs only after period has been drawn. We use two criteria to determine ${e}_{\max }$: the first given by Equation (5) and the second from Equation (3) of Moe & Di Stefano (2017). The criterion from Moe & Di Stefano (2017),

$$\begin{eqnarray*}\begin{array}{r}\begin{array}{c}{{e}_{max}}_{,{\rm{M}}{\rm{D}}{\rm{S}}}=\left\{\begin{array}{cc}0 & {\rm{f}}{\rm{o}}{\rm{r}}\,P\leqslant 2\,{\rm{d}}{\rm{a}}{\rm{y}}{\rm{s}}\\ 1-{\left({\textstyle {\textstyle \tfrac{P}{2{\rm{d}}{\rm{a}}{\rm{y}}{\rm{s}}}}}\right)}^{-2/3} & {\rm{f}}{\rm{o}}{\rm{r}}\,P\gt 2\,{\rm{d}}{\rm{a}}{\rm{y}}{\rm{s}},\end{array}\right.\end{array}\end{array}\end{eqnarray*}$$

assumes everything below 2 days is circularized and guarantees that the binary components do not fill their Roche lobes by more than 70%. Combining this with Equation (5), we have:

$$\begin{eqnarray}&&{e}_{\max }=\min \left(1-\displaystyle \frac{1.5({R}_{1}+{R}_{2})}{a},{{e}_{\max }}_{,\mathrm{MDS}}\right)\end{eqnarray} \tag{ 6 }$$

where we use whichever criterion provides the more conservative value for ${e}_{\max }$.

The final properties to be drawn are inclination and argument of periastron. Inclination is sampled uniformly in terms of $\cos i$ for i between 0° and 180°. The argument of periastron, ω, is drawn uniformly between 0° and 360°.

Once we have a sample of synthesized binary systems, we need to determine which of those will eclipse. We compute the projected separation of the system and use

$$\begin{eqnarray*}&&{R}_{1}+{R}_{2}\gt r(\nu )\cos i\end{eqnarray*}$$

as the criterion for an eclipse where r(ν) is the instantaneous separation of the stellar centers, as a function of true anomaly, given by

$$\begin{eqnarray*}&&r(\nu )=\displaystyle \frac{a(1-{e}^{2})}{1+e\cos (\nu )}.\end{eqnarray*}$$

For circular orbits, there are two critical points, ν_crit, that correspond to the minimum projected separation of the system. These points occur at ${\nu }_{\mathrm{crit}}=\tfrac{\pi }{2}-\omega$ and ${\nu }_{\mathrm{crit}}=\tfrac{3\pi }{2}-\omega$. For eccentric orbits the point of projected closest approach cannot be solved analytically. The criterion we use to determine if a system is a synthetic eclipsing binary is

$$\begin{eqnarray}&&\begin{array}{l}{R}_{1}+{R}_{2}\gt \displaystyle \frac{a(1-{e}^{2})}{1+e\cos ({\nu }_{\mathrm{crit}})}\\ \times \,\sqrt{{\cos }^{2}({\nu }_{\mathrm{crit}}+\omega )+{\sin }^{2}({\nu }_{\mathrm{crit}}+\omega ){\cos }^{2}i},\end{array}\end{eqnarray} \tag{ 7 }$$

where the right side of Equation (7) is the projected separation of the system in terms of the orbital parameters. We evaluate Equation (7) at both values of ν_crit to determine if an eclipsing event occurs.

Currently, we do not differentiate between binaries which have only a single eclipsing component versus both components. Our eclipse criterion assumes spherical stars, which is appropriate given that all systems have been generated to ensure that they are detached. In addition, we do not consider the depth, nor the profile, of the eclipsing signal which would require a Kepler light-curve noise model and a limb-darkening model, respectively. The systems that satisfy the geometry constraint, Equation (7), are considered to be eclipsing.

Due to telemetry restrictions, fewer than 200,000 targets had data collected. In addition, the process of target selection was not random, but was specifically chosen to optimize the detection of Earthlike planets about Sunlike stars. This biased the target list toward FGK-type main-sequence stars. Magnitudes and colors obtained from surveys such as the Sloan Digital Sky Survey (SDSS) and the Two Micron All Sky Survey (2MASS) were used to perform the selection. To simulate the complex target selection process we draw systems from our synthetic sample (single and binary systems) that are as similar as possible, in terms of magnitude and color, to the Kepler Stellar Properties Catalog (KSPC; Mathur et al. 2017).

The KSPC contains, in addition to the broadband visible Kepler magnitude, K_p , values for the 2MASS infrared J, H, and K_s bands. Galaxia provides absolute magnitudes for the SDSS g- and r- bands and the 2MASS J, H, and K_s infrared bands along with the extinction coefficient E(B − V) for each star. We compute the apparent magnitudes for each of these bands, while also determining the combined magnitude for each binary treating all binaries as unresolved by Kepler. To determine the Kepler magnitudes we use the SDSS g- and r-bands to determine K_p via the relationships provided by Brown et al. (2011):

$$\begin{eqnarray}{K}_{p}=\left\{\begin{array}{ll}0.2g+0.8r & \mathrm{for}\ (g-r)\leqslant 0.8\\ 0.1g+0.9r & \mathrm{for}\ (g-r)\gt 0.8.\end{array}\right.\end{eqnarray} \tag{ 8 }$$

We determine the synthetic target list using the following merit function:

$$\begin{eqnarray}&&{S}_{{ts}}=\sqrt{{\left({K}_{p}^{* }-{K}_{p}\right)}^{2}+{\left({J}^{* }-{K}_{s}^{* }-J+K\right)}^{2}+{\left({H}^{* }-H\right)}^{2}},\end{eqnarray} \tag{ 9 }$$

where we denote the parameters from the synthetic population with an asterisk while the parameters without an asterisk are the parameters from the KSPC. Systems are drawn from the synthetic population, without replacement, for each target in the KSPC.

The J − K_s versus H distribution shows the presence of a large population of giants. In the left panel of Figure 4, the strip of systems with J–K_s ≈ 0.8 are red giants. The right middle panel of Figure 4 shows that the catalog contains a small clump of faint red giants, J–K_s > 0.8 and H > 13, that is not reproduced by our simulated selection process.

Zoom In Zoom Out Reset image size

Using the conversions provided by Caldwell et al. (1993) and the approximation by Ballesteros (2012), we find that a star with J–K_s ≈ 1.5 would have an effective temperature of 3200 K. The relationships by Caldwell et al. (1993) do not hold for values of J–K_s > 1.5 but it is safe to say that these would correspond to nonstellar sources. The relationships used to derive the Kepler magnitude have systematic errors as high as 0.6 mag toward fainter Kepler magnitudes for cooler stars and could account for the clustering features seen in the difference between our synthetic target list and the observed target list. Looking at the two main selection breaks that occur at 14th and 16th mag in Figure 5, we can see that our process underselects brighter targets and overselects fainter systems. For a given H magnitude the corresponding K_p will be overestimated for cooler systems, which explains the discrepancies in Figure 5.

Zoom In Zoom Out Reset image size

The instrument also suffered from gaps in the data. Some of these gaps were scheduled, such as rolling of the spacecraft to reorient the solar panels toward the Sun and to download data, while others, like the failure of two CCDs, were not. An empirical two-eclipse detection efficiency model is provided by Kirk et al. (2016), see Figure 6.

Zoom In Zoom Out Reset image size

Their empirical detection efficiency model was obtained by scanning all of the Kepler light curves over a range of observable test periods. We modify their empirical relationship by first considering the scenario of a periodic duty cycle.

The probability of observing exactly k eclipses out of N eclipsing events is given by the binomial distribution

$$\begin{eqnarray*}&&{p}_{k}=\left(\displaystyle \genfrac{}{}{0em}{}{N}{k}\right){f}_{{\rm{dc}}}^{k}{\left(1-{f}_{{\rm{dc}}}\right)}^{N-k}\end{eqnarray*}$$

where f_dc is the duty cycle. The number of eclipses that can occur for a given period is

$$\begin{eqnarray*}&&N=1+\lfloor (1460\,\mathrm{days})/P\rfloor \end{eqnarray*}$$

where 1460 days is the length of the original Kepler mission. The KEBC only includes systems that have had three eclipses observed. The probability of at least three observed eclipses is the same as

$$\begin{eqnarray*}&&{p}_{k\geqslant 3}=1-{p}_{0}-{p}_{1}-{p}_{2}\end{eqnarray*}$$

where p₀, p₁, and p₂ are the probabilities of observing no eclipses, only one eclipse, and exactly two eclipses, respectively. Our modified three-eclipse detection efficiency function is, using f_dc = 0.92,

$$\begin{eqnarray*}&&{p}_{k\geqslant 3}={f}_{{\rm{emp}}}-0.5N(N-1){\left(0.92\right)}^{2}{\left(0.08\right)}^{(N-2)},\end{eqnarray*}$$

where we have replaced the uniform probability of observing at least two eclipses with the empirical relationship from Kirk et al. (2016), f_emp.

The process of fitting the $\mathrm{log}P$ distribution begins by vastly overestimating the total number of systems. Figure 8 shows the SEBC after each run.

Zoom In Zoom Out Reset image size

Each bin is updated after every iteration and after roughly 15 runs the total number of synthetically observed systems matches the total number from the KEBC. We run 250 iterations to ensure that the model has burned in about the solution. To get our resulting model, we perform an additional 50 iterations from which we compute the mean bin values and associated uncertainties.

The resulting binary parameter distributions are shown in Figure 9. The mass distribution shows a decline at low masses. This drop off is a result of using a magnitude-limited sample. The mass distribution of the target-selected sample peaks around a solar mass. This is in line with the target selection process and shows that our process of simulated target selection captures the distribution of targets appropriately. The selected eclipsing population is dominated by binaries with M₁ > 0.5 M_⊙. This is due to a combination of effects. The first is the declining fractional eclipse probability for smaller, and typically less massive stars, which can be seen if one compares the shape of the eclipsing binaries to the binaries. Eclipse probability is proportional to the sum of fractional radii so smaller stars will present eclipses less frequently than larger stars. When combined with the target selection bias, eclipsing binaries below 0.5 M_⊙ and 0.4 R_⊙ are effectively suppressed.

Zoom In Zoom Out Reset image size

The actual target list was not free of giants, but included roughly 5000 stars with $\mathrm{log}g\lt 3.5$. The synthetic target-selected radius distribution is in agreement, having less than 2000 binaries with R₁ > 10 R_⊙. There appears to be a discrepancy between the target-selected eclipsing binary sample and the observed eclipsing binaries which widens for larger stars. This follows directly from our minimum period constraint in Equation (5): systems with larger stars are more likely to have values of $\mathrm{log}P\gt 2.86$ and are unable to have three eclipses observed.

The $\mathrm{log}P$ distribution appears to be relatively flat above 0.5 (roughly 3.2 days). We estimate that the period range covered by Kepler contains 12.0(2)% and 83.2(2)% of all binaries and eclipsing binaries in the field, respectively. The fraction of eclipsing binaries to binaries declines at longer periods because of the larger separation of the system which reduces eclipse efficiency. The discrepancy between the selected and observed eclipsing binary samples at longer $\mathrm{log}P$ is due to the lower detection efficiency. The obtained mass ratio distribution, M₂/M₁, reflects the fixed target mass ratio distribution q.

Figure 10 shows the input distributions that were held fixed throughout the synthesis. The eccentricity distributions of the simulated underlying binary population and target-selected binary sample show a slight abundance in the first bin due to circularization. These distributions show that eclipse detection is boosted for circularized systems and high-eccentricity systems. Circularized systems tend to have shorter periods than eccentric systems and thus have a larger eclipse detection efficiency. High-eccentricity systems will have a boost in detection efficiency because of the closer approach of the systems near periastron. The discrepancy in the eccentricity distributions for selected eclipsing binaries and the observed eclipsing binaries is related to the constraint from Equation (5). Higher eccentricities require longer periods making these systems more likely to have values of $\mathrm{log}P$ outside of the observable window. The decline of the number of systems for smaller mass ratios is due to the minimum mass ratio constraint provided in Equation (3). In addition, smaller and larger mass ratios are overrepresented in the target list as compared to more moderate mass ratios. This is due to the shape of the mass ratio distribution with respect to primary mass, presented in Figure 3. Smaller mass ratios are going to be present only in systems with more massive primaries making these more likely to make it into the target list. Higher mass ratio systems will tend to be brighter and will be slightly favored for a given primary mass.

Zoom In Zoom Out Reset image size

The i and ω parameters pertain to the viewing angle and show no unexpected features for the binary population. The i distribution shows that eclipses only occur between 45° and 135° which is precisely what one would expect for spherical stars. The ω distribution shows that eclipses are more likely when periastron is aligned with the plane of the sky, which occurs at 90° and 270° as opposed to when periastron is perpendicular to the plane of the sky, which occurs at 0° and 180°.

In order to validate our resulting $\mathrm{log}P$ distribution, we compare it to the volume-limited sample of Raghavan et al. (2010), as reported in Moe & Di Stefano (2017). While the Raghavan et al. (2010) sample is composed of solar-like stars found within 25 pc of the Sun, it still provides a useful point of comparison for our results. Figure 11 shows that the $\mathrm{log}P$ distribution of our model is consistent with Raghavan et al. within uncertainty.

Zoom In Zoom Out Reset image size

A power law of the form f(P) = kP^α was fitted to the derived distribution over the range of 2–730 days, with k = 0.044 ± 0.003 and α = 0.00 ± 0.02. This range was chosen to avoid circularization effects (see Equation (6)) at shorter periods. The resulting fit shows a flat curve consistent with Raghavan et al. (2010) within 1σ uncertainty.

We investigated how our choice of input distributions, namely mass ratio and eccentricity, affected $\mathrm{log}P$. Figure 12 shows three separate input mass ratio distributions and the resulting $\mathrm{log}P$ distribution for each. We conclude from the difference plot that the resulting $\mathrm{log}P$ distribution is only marginally coupled to the mass ratio distribution (up to 2σ deviations). The overflow bin percentage also shows the relatively low coupling between mass ratio with less than a percent variation across all three distributions. Figure 13 shows the three separate input eccentricity distributions and the resulting $\mathrm{log}P$ distribution for each.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Steeper eccentricity distributions result in an attenuation of the $\mathrm{log}P$ distribution in the observable period range by shuffling the systems to even longer periods. The difference plot shows that these eccentricity distributions do result in distinct period distributions. Figure 14 shows the simulated underlying binary distribution of e as a function of $\mathrm{log}P$.

Zoom In Zoom Out Reset image size

These distributions have identical boundaries because of the restriction on ${e}_{\max }$. The non-flat distributions contain very fewer low eccentricity systems above $\mathrm{log}P\gt 0.5$. Ongoing work (A. Prša et al. 2021, in preparation) will provide accurate eccentricities for the majority of the KEBC.