BINARY POPULATIONS IN MILKY WAY SATELLITE GALAXIES: CONSTRAINTS FROM MULTI-EPOCH DATA IN THE CARINA, FORNAX, SCULPTOR, AND SEXTANS DWARF SPHEROIDAL GALAXIES

Recent progress in hydrodynamical simulations has allowed, for the first time, numerical simulations of star formation that include all the relevant physics down to solar system scales (Bate 2009; Offner et al. 2009). Star formation is now understood to occur due to the gravitational collapse and fragmentation of a turbulent molecular cloud, generally in the presence of radiative feedback (Bate 2012) and magnetic fields (Price & Bate 2010). Simulations of star formation that incorporate these physical effects now offer detailed predictions about the statistical properties of stellar systems that must be tested against observational data to arrive at a complete theory of star formation. A major component of such predictions involve the statistical properties of binary star systems. While approximately half of solar-type field stars in the Milky Way are known to have binary companions, the fraction of young pre-main sequence stars in binary or higher-order star systems is much higher and possibly greater than 90% (Leinert et al. 1993; Kohler & Leinert 1998; Patience et al. 2002). This binary fraction is later reduced as many wide binaries are subsequently disrupted by cluster dynamics or orbital decay. Thus, binary systems are likely the dominant mode of star formation, and the ability to successfully predict their statistical properties will be an essential component of a complete theory of star formation.

At present, a complete census of the statistical properties of binary star systems exists only in the solar neighborhood, at distances out to ≈30 pc from the Sun (Raghavan et al. 2010; Duquennoy & Mayor 1991). While a number of visual binary studies have been made in nearby open clusters (Brandner & Koehler 1998; Patience et al. 2002), these studies are limited to wide binaries with separations ≳10 AU. Moreover, such wide binaries are likely to have been "dynamically processed" by stellar encounters within the cluster, and thus may not represent the primordial binary populations present shortly after star formation has occurred. Likewise, spectroscopic and photometric studies have been made in globular clusters (Yan & Cohen 1996; Milone et al. 2012; Sollima et al. 2007, 2012), but in many of these systems even close binaries have been dynamically processed in the dense core of the cluster over cosmological timescales. To test theories of star formation, it is important to study binary populations in diffuse clusters or field populations over a range of separations for which the binary systems can be considered primordial, and compare their statistical properties to that of simulations.

Dwarf spheroidal (dSph) galaxies are excellent objects for this purpose. Because they are more diffuse than large globular clusters, stellar encounters are sufficiently rare in dSphs that binary stars with separations less than 10 AU should be largely unaltered since their formation. Dwarf spheroidals represent a variety of star forming environments, ages, and metallicities, against which theories of star formation can be tested. Furthermore, as simulations advance and our understanding of star formation is refined, dSphs may ultimately become laboratories for a new brand of "stellar archaeology": using the statistical properties of their binary star populations to discern the initial conditions under which the galaxy originally formed. The extent to which binary properties are sensitive to initial conditions during star formation is still an open question. While differences between binary populations of different clusters have been observed (King et al. 2012), these differences may be largely due to subsequent dynamical processing with the cluster (Marks & Kroupa 2012; Kroupa 1995). If binary star properties do indeed have a significant dependence on initial conditions (e.g., cloud temperature, density, magnetic field, metallicity), they may eventually become a tool for understanding star formation histories.

Another reason for constraining binary properties in dwarf galaxies is to understand the effect of binary orbital motion on the stellar velocity distribution, and hence the estimated mass distribution, of the galaxies (Olszewski et al. 1996; Minor et al. 2010). The spatial mass distribution of dark matter within a galaxy, including the existence and size of a central core, is determined by properties of the dark matter particle itself. However, because there exists a degeneracy between mass and velocity anisotropy (Wolf et al. 2010), it is not possible to infer the dark matter density profile in a straightforward way. One way to overcome this degeneracy is to make use of higher moments in the velocity distribution to distinguish between radial and tangential orbits, and hence infer the velocity anisotropy (Łokas et al. 2005; Łokas 2009; Richardson & Fairbairn 2013). However because binary motion can also contribute significantly to higher moments in the velocity distribution, it is important to correct for binaries when determining the velocity anisotropy. Therefore, quantifying the effect of binary motion on a galaxy's velocity distribution can lead to a better understanding of the distribution of dark matter in dwarf galaxies.

Finding detailed constraints on the binary populations in these galaxies is a difficult task. In the most sensitive spectroscopic surveys of dwarf spheroidals to date, typical line-of-sight velocity measurement errors (≈1–3 km s⁻¹ for red giant stars) are large enough to render binary motion in systems with periods longer than 10 yr unobservable. Since binaries with such long periods cannot be directly observed in these galaxies, any statistical information about long-period binaries can only be inferred by extrapolation. Furthermore, among the red giant population, information about the shortest-period binaries is lost because close binary systems are destroyed when the red giant star overflows its Roche lobe (Paczyński 1971). Another complication is that any component of a star's velocity due to binary motion must be separated from the measurement error, and therefore deriving accurate and reliable measurement errors is essential before obtaining binary constraints.

In this paper, we show that all of these difficulties can be surmounted by a likelihood approach, and demonstrate the method using multi-epoch data from the Magellan/MMFS sample of Walker et al. (2009a) in the Carina, Fornax, Sculptor, and Sextans dwarf spheroidal galaxies. In Section 2 we derive a multi-epoch likelihood for binary and non-binary stars, followed by a discussion of binary model uncertainties in Section 3. In Section 4, we show how this likelihood can incorporate an error model to derive accurate measurement errors in the Magellan/MMFS sample, which we apply to the data in Section 5. To ensure the error model has properly reproduced the measurement errors, we then test the model on simulated data in Section 6. Next, the binary fraction in each galaxy is inferred in Section 7, under the assumption of a Milky Way-like period distribution, and our main results are plotted in Figure 5, wherein probability distributions in binary fraction are plotted for each galaxy. The inferred binary fraction of the combined sample of all four galaxies is plotted in Figure 6, and the best-fit binary fractions in each galaxy are listed in Table 2. In Section 8 we find more general constraints on the binary fraction, mean period, and spread of periods in each galaxy. In Section 9 we apply our methodology to simulated data to investigate the prospects for finding more precise binary constraints in a larger multi-epoch dataset. Finally, we discuss implications of our results in Section 10 and conclude in Section 11.

In order to constrain binary properties from radial velocity measurements, we must find the likelihood of our binary model parameters for a given set of measured velocities. If the fraction of stars in binary systems B is taken as a model parameter, then the likelihood will take the form $\mathcal {L} = \mathcal {L}_{nb} + B\mathcal {L}_b$, where $\mathcal {L}_{nb}$ and $\mathcal {L}_b$ represent the likelihood for non-binary and binary stars, respectively. In addition, we will choose a set of parameters $\mathcal {P}$ that specify the distributions of binary properties. These binary properties may include the orbital periods, mass ratios, and orbital eccentricities. The assumed distributions in these properties and our chosen model parameters will be described in Section 3.

In principle, even in a sample of stars with only one velocity measurement each, binary properties could be inferred from the line-of-sight velocity distribution of the sample, where the presence of binary stars is evident from a non-Gaussian tail representing short-period binaries with high orbital velocities. However, inferring binary properties from the velocity distribution is fraught with difficulties. Because the classical dwarf spheroidal galaxies have velocity dispersions greater than 7 km s⁻¹, only stars with radial velocities that differ from the galaxy's systemic velocity by at least ≳14 km s⁻¹ will be evident as binary stars, and only a small fraction of binaries have velocity amplitudes this large. For this reason, extremely large samples would be required to have any hope of constraining binary properties from single-epoch velocity measurements alone. Contamination from Milky Way stars must also be accounted for, as such stars may be mistaken for binaries. Finally, non-Gaussianities may be present in the galaxy's velocity distribution due to velocity anisotropy and must also be distinguished from binary motion.

Because this work is primarily concerned with finding binary constraints, we can avoid the difficulties inherent in using the galaxy's velocity distribution by focusing on the stars with repeat measurements. Suppose a star has n radial velocity measurements {v_i} = {v₁, ⋅⋅⋅, v_n} and measurement errors {σ_i} taken at the corresponding times, {t_i}. For stars with at least two measurements (n ⩾ 2), we shall derive a likelihood in velocity differences Δv₁, Δv₂, and so on, where we define Δv_i = v_{i + 1} − v₁. For the sake of readability, we will henceforth suppress the brackets denoting sets of measurements (e.g., P({v_i}) → P(v_i)).

Our desired likelihood is derived as follows. We first write down the probability, for single stars, of obtaining a set of radial velocity measurements {v_i} assuming that it has intrinsic velocity v_true. Assuming the measurement errors are Gaussian, this likelihood can be expressed as follows:

$$\begin{eqnarray} P(v_i|v_{{\rm true}},\sigma _i) & = & \prod _{i=1}^n \frac{e^{-(v_i-v_{{\rm true}})^2/2\sigma _i^2}}{\sqrt{2\pi \sigma _i^2}} \nonumber \\ & = & \mathcal {N}(v_i,\sigma _i) \frac{e^{-(v_{{\rm true}} - \langle v \rangle)^2/2\sigma _m^2}}{\sqrt{2\pi \sigma _m^2}}. \end{eqnarray} \tag{ 1 }$$

Here, 〈v〉 and σ_m are the variance-weighted average velocity and measurement error,

$$\begin{equation} \langle v \rangle = \sigma _m^2 \sum _{i=1}^n \frac{v_i}{\sigma _i^2}, \end{equation} \tag{ 2 }$$

$$\begin{equation} \sigma _m^2 = \left(\sum _{i=1}^{n}\frac{1}{\sigma _i^2}\right)^{-1}, \end{equation} \tag{ 3 }$$

and a tedious algebraic calculation shows that the factor $\mathcal {N}(v_i,\sigma _i)$ is given by

$$\begin{eqnarray} \mathcal {N}(v_i,\sigma _i) &=& \frac{\sqrt{2\pi \sigma _m^2}}{\prod _{i=1}^n \sqrt{2\pi \sigma _i^2}} \nonumber \\ && \times \exp \left\lbrace \frac{-1}{4}\sum _{i,j=1}^n\frac{(v_i-v_j)^2}{\sigma _i^2 + \sigma _j^2 + \sigma _i^2 \sigma _j^2\left(\sum _{k \ne i,j}\frac{1}{\sigma _k^2}\right) } \right\rbrace. \nonumber\\&& \end{eqnarray} \tag{ 4 }$$

The final term in the denominator of the exponent is explicitly zero when n = 2. Note that this factor depends only on velocity differences v_i − v_j; it has no dependence on the star's intrinsic velocity v_true.

To derive a likelihood in velocity differences Δv, we change variables in Equation (1) from (v₁,v₂,v₃,...) to (〈v〉,Δv₁,Δv₂,...) and integrate over 〈v〉. Then the Gaussian factor in Equation (1) integrates to 1, and we are left with

$$\begin{eqnarray} \mathcal {L}_{nb}(\Delta v_i|\sigma _i) &&= \frac{\sqrt{2\pi \sigma _m^2}}{\prod _{i=1}^n \sqrt{2\pi \sigma _i^2}} \nonumber \\ &&\times \exp \left\lbrace \frac{-1}{4}\sum _{i,j=1}^n\frac{(\Delta v_{i-1}-\Delta v_{j-1})^2}{\sigma _i^2 + \sigma _j^2 + \sigma _i^2 \sigma _j^2\left(\sum _{k\ne i,j}\frac{1}{\sigma _k^2}\right) } \right\rbrace. \nonumber\\&& \end{eqnarray} \tag{ 5 }$$

This is our desired likelihood for non-binary stars; it is essentially the same factor $\mathcal {N}(v_i,\sigma _i)$ above, but expressed in the new variables Δv_i. As an aside, we note that it is possible to derive a likelihood in all the velocities v_i using this formalism; this is the approach taken in Martinez et al. (2011). Since we are only interested in binary constraints rather than the galaxy's intrinsic velocity distribution, only a likelihood in velocity differences is considered in the present work.

The corresponding likelihood for binary stars cannot be found analytically, so it must be calculated by running a Monte Carlo simulation for a large number of simulated stars and binning in the velocity differences Δv_i (the details of the Monte Carlo simulation are discussed in Minor et al. 2010). The likelihood depends not only on the velocities and measurement errors, but also the time intervals between epochs and the absolute magnitude M_V of the star, from which we determine the mass and radius of the primary star via a stellar population synthesis model. The likelihood will also depend on the set of parameters $\mathcal {P}$ used to characterize the binary population which we discuss further in Section 3. In terms of all the relevant model parameters, the likelihood for each star is therefore

$$\begin{eqnarray} &&\mathcal {L}(\Delta v_i|\sigma _i,t_i,M;B,\mathcal {P}) = (1-B)\mathcal {L}_{nb}(\Delta v_i|\sigma _i) \nonumber\\ &&+\,m B\mathcal {L}_b(\Delta v_i|\sigma _i,t_i,M;\mathcal {P}). \end{eqnarray} \tag{ 6 }$$

We can put this in a somewhat simpler form by noting that the likelihood for single stars does not depend on any model parameters, and thus can be divided out of Equation (6). We then have

$$\begin{equation} \mathcal {L}(\Delta v_i|\sigma _i,t_i,M;B,\mathcal {P}) \propto (1-B) + B J(\mathcal {P}), \end{equation} \tag{ 7 }$$

where

$$\begin{equation} J(\mathcal {P}) = \frac{\mathcal {L}_b(\Delta v_i|\sigma _i,t_i,M;\mathcal {P})}{\mathcal {L}_{nb}(\Delta v_i|\sigma _i)}. \end{equation} \tag{ 8 }$$

For non-binary stars or binaries with small velocity amplitudes compared to the measurement errors, the J-factor in Equation (7) will be of order 1 or less, since the binary likelihood will typically be slightly smaller than the non-binary likelihood. However, if a binary star exhibits large velocity changes compared to the measurement errors, the J-factor will be larger than 1, possibly by many orders of magnitude, since the non-binary likelihood $\mathcal {L}_{nb}$ in the denominator will be very small.

For each star, the factor $J(\mathcal {P})$ can be calculated over a grid of values in the binary parameters $\mathcal {P}$, after which the J-factors can then be found by interpolation for any values of the binary parameters $\mathcal {P}$. After calculating the J-factors, the binary fraction B and parameters $\mathcal {P}$ can then be inferred either by a maximum-likelihood or Bayesian analysis.

Here we confront the issue of which binary population model parameters to constrain and how to deal with uncertainties in these parameters. In addition to the binary fraction, binary star populations are described by distributions in three properties: the orbital period P, mass ratio q, and eccentricity e. (We have assumed a uniform distribution of orbital inclinations, as is observed in Milky Way binaries; see Minor et al. 2010.) Without a large set of measurements, orbital eccentricities are difficult to constrain because binaries with high eccentricity spend only a small fraction of their orbital period near periastron where their observed velocities are high. For this reason, we will assume the distribution of eccentricities to have the form given in Minor et al. (2010), with parameters fixed to those observed in solar neighborhood binaries.

Similarly, several measurements are typically required to constrain the mass ratio of a binary separately from its orbital period. On the other hand, theoretical considerations and open cluster surveys suggest that while the period distributions of binary populations may vary significantly depending on initial conditions (Fisher 2004; Brandner & Koehler 1998), the distribution of mass ratios may be approximately universal in form. This is known to be the case for binaries with long periods (P > 1000 days), where the distribution of mass ratios traces the Salpeter initial mass function to good approximation for q ≳ 0.5 (assuming the primary star masses fall within a restricted range, as it does for stellar populations older than a few Gyr; see Duquennoy & Mayor 1991). The mass ratio distribution observed in shorter-period binaries is closer to being uniform (Goldberg et al. 2003; Mazeh et al. 1992), although it is uncertain whether this distribution is universal in form. In light of the above considerations, we assume the distribution of mass ratios to take a fixed form similar to that observed in the solar neighborhood, described in Minor et al. (2010), with a Gaussian form for long-period binaries and a flat distribution for short-period binaries.

An important caveat is that if the mass distribution differs from our assumed form, this can introduce bias in the results. For example, if the actual mass distribution for short-period binaries rises steeply at high mass ratios, instead of being flat, this would result in greater binary velocity variations. Therefore our analysis would infer a somewhat higher binary fraction than is actually present under the assumption of a flat mass ratio distribution. In principle, we could also vary the parameters describing the mass ratio distribution to account for this, but this is computationally intensive; the binary likelihoods will have to be computed over a grid of binary parameters, including the parameters used to describe the mass ratio distribution. In addition, without a large number of epochs to separate each star's binary period from its velocity amplitude, these parameters will be highly degenerate with other binary parameters (specifically, the binary fraction and period distribution parameters) and thus we are unlikely to infer anything about the behavior of the mass-ratio distribution. This degeneracy is similar to that of the period distribution parameters and binary fraction, which will be discussed in detail in Section 9.

By way of analogy to Milky Way field binaries, we assume the period distribution of dwarf spheroidal galaxies to have a log-normal form (Duquennoy & Mayor 1991; Fischer & Marcy 1992; Mayor et al. 1992; Raghavan et al. 2010). Initially we shall fix the mean period μ_log P and dispersion of periods σ_log P to the values observed in solar neighborhood binaries (where μ_log P = 2.24, σ_log P = 2.3), while in later sections (8–9) these parameters will be allowed to vary. In general, we therefore have three binary parameters that will be constrained: the binary fraction B, the mean log-period μ_log P, and log-spread of periods σ_log P. It should be borne in mind that for the dwarf galaxy samples considered in this paper, binary motion in stars with periods ≳10 yr will be unobservable due to the measurement error (the correspondence between observable binary periods and measurement error will be explored in detail in Section 9.3). Thus we can only directly probe the period distribution for periods shorter than ∼10 yr, so the shape of the period distribution at longer periods must be inferred by extrapolation. This extrapolation will hold only to the extent that the period distribution takes a symmetric log-normal form, which we take as a working hypothesis in this paper.

Since any velocity variations beyond what is indicated by the measurement errors may be interpreted as binary motion, to constrain binary properties it is essential to estimate the measurement errors as well as possible. In the data from Walker et al. (2009a), measurement errors are determined from repeat velocity measurements by applying an error model, which will be reviewed and extended in this section. It is important to note, however, that if measurement errors are estimated by means other than velocity variability (e.g., from signal-to-noise (S/N) properties), an error model such as the one considered here need not be applied before obtaining binary constraints.

In the error model from Walker et al. (2009a), Gaussian measurement errors σ are estimated using the model

$$\begin{equation} \sigma = \sqrt{\sigma _0^2 + \sigma _{{\rm CCF}}^2}, \end{equation} \tag{ 9 }$$

$$\begin{equation} \sigma _{{\rm CCF}} = \frac{\alpha }{(1+R)^x} \end{equation} \tag{ 10 }$$

where σ₀ is a baseline error, σ_CCF is the spectrum cross-correlation error, and R_i is called the Tonry–Davis R-value for the measurement (Tonry & Davis 1979), defined as the height of the maximum peak in the spectrum cross-correlation function divided by the average peak height. This error model is an extension of the error model from Tonry & Davis (1979), which takes x = 1. The MIKE spectrograph has two channels, "red" and "blue," which operate over different wavelength ranges, so we must have two sets of error model parameters representing each channel. Our error model then has four parameters given by the set $\mathcal {E} = \lbrace \alpha _{{\rm red}}$, α_blue, x_red and x_blue}. The baseline error also depends on the channel color, and was determined to be σ_{0, red} = 0.6 km s⁻¹ and σ_{0, blue} = 0.26 km s⁻¹.

The implementation of the error model described above in Walker et al. (2007) has a few disadvantages. First, their method of constraining the error model parameters in makes certain assumptions about the intrinsic velocity of each star before deriving measurement errors. Since the star's intrinsic velocity cannot be well determined from repeat measurements without knowing the measurement errors in the first place, this may lead to bias in the error model.

In contrast, since our multi-epoch likelihood is expressed solely in terms of velocity differences, we can naturally incorporate the error model with no assumptions about intrinsic velocity. If the star's measurements all have time intervals of less than a week between them, binary motion can be neglected and we can incorporate the error model into our likelihood by expressing the measurement error in terms of the error model parameters. Combining Equations (9) and (10), we can write the error in the ith measurement of the star as

$$\begin{equation} \sigma _i = \sqrt{\sigma _{0,c_i}^2 + \frac{\alpha _{c_i}^2}{(1+R_i)^{2x_{c_i}}}}, \end{equation} \tag{ 11 }$$

where c_i is the channel color of the ith measurement (either "red" or "blue"). This expression is then inserted into the likelihood from Equation (5), which we now write as

$$\begin{eqnarray} \mathcal {L}(\Delta v_i|R_i&,&c_i; \mathcal {E}) = \frac{\sqrt{2\pi \sigma _m^2}}{\prod _{i=1}^n \sqrt{2\pi \sigma _i^2}} \nonumber \\ &&\times \exp \left\lbrace \frac{-1}{4}\sum _{i,j=1}^n\frac{(\Delta v_{i-1}-\Delta v_{j-1})^2}{\sigma _i^2 + \sigma _j^2 + \sigma _i^2 \sigma _j^2\left(\sum _{k\ne i,j}\frac{1}{\sigma _k^2}\right) } \right\rbrace.\nonumber\\&& \end{eqnarray} \tag{ 12 }$$

More importantly, the error model in Walker et al. (2007) does not take binary motion into account, which may bias the errors if the time between measurements is longer than a few weeks. We shall address this by incorporating binarity into the likelihood of each star if at least one of the time intervals between measurements is long enough for binary velocity variation to be observed. To determine this minimum time interval, we note that for intervals smaller than ∼10 days, the distribution of observed velocity changes is well-approximated by a Gaussian; thus we choose 10 days as our minimum threshold for observed binary motion. Since we are now using the binary likelihood, we must constrain the binary fraction B and the error model parameters simultaneously.

At this stage, the reader may be wondering why our binary model is being applied twice: first to find the measurement errors, then afterward to constrain the binary fraction B and binary model parameters $\mathcal {P}$. Why not constrain all the error model and binary parameters simultaneously? The answer is that this would be extremely computationally expensive, because the binary likelihoods must be recalculated every time the error model parameters are varied. Fortunately, if we are (for the moment) only interested in deriving the best possible measurement errors, it is sufficient to assume a Milky Way-like period distribution and constrain only the binary fraction B in order to distinguish binary velocity variation from measurement error. Furthermore, as will become evident shortly, the binary likelihoods become computationally expensive to calculate while varying the error parameters if a star contains more than one "long" time interval (≳10 days) between measurements. Thus, any extra measurements occurring at later time intervals longer than 10 days will be discarded for the purpose of determining the error model parameters, although they will be used later to help determine the best possible binary constraints.

In summary, while it may appear cumbersome to apply the binary model twice, the benefit of finding more detailed binary constraints after determining the measurement errors is that ultimately the binary parameters $\mathcal {P}$ will be constrained in addition to the binary fraction B, and all the available measurements will be used to make this determination.

To incorporate binary motion into our likelihood, we need to be able to quickly calculate the binary likelihood for a given set of assumed measurement errors, which we continually update as the error model parameters are changed. To accomplish this, we first run a Monte Carlo simulation to generate the binary likelihood in log |Δv| with zero measurement error, which we denote by L_b(log |Δv|; Δt, M_V) where M_V is the absolute magnitude of the star and Δt is the time interval between measurements. The likelihood is generated in log |Δv|, rather than Δv, because the binary likelihood becomes singular as Δv  →  0. The binary likelihood is generated over a table of absolute magnitude values M_V and time intervals Δt. In practice the long time intervals are typically close to a multiple of one year, so for our purposes it is sufficient to pick time intervals of 1 yr, 2 yr, and so on. For each star with a long time interval between measurements, we interpolate in magnitude and choose the time interval (in multiples of one year) closest to the actual interval between measurements. It should be emphasized that in later sections, we will use the exact time intervals between measurements when deriving detailed binary constraints; however, for the purpose of deriving measurement errors, quantizing the time intervals in this way is less computationally intensive, and does not impact the binary modeling or the measurement error results significantly.

Now suppose that the star has n velocity measurements, and the long time interval is between measurement k and k + 1. To find the binary likelihood in the presence of measurement error, we must find the convolution of the binary likelihood and the likelihood from Equation (12) (representing measurement error) in the velocity difference Δv_k. By switching to an integral over $\log |\Delta v_k^{\prime }|$, this can be shown to equal

$$\begin{eqnarray} \mathcal {L}_b(\Delta v_i|R_i&,&c_i,\Delta t_i,M;\mathcal {E}) \nonumber \\ = && \frac{1}{2}\int _{-\infty }^{\infty } L_b(\log |\Delta v_k^{\prime }|;\Delta t_k,M) \nonumber \\ && \times \left \lbrace \mathcal {L}(\cdots,\Delta v_k - |\Delta v_k^{\prime }|,\cdots |R_i,c_i;\mathcal {E})\right. \nonumber \\ &&\left. +\, \mathcal {L}(\cdots,\Delta v_k + |\Delta v_k^{\prime }|,\cdots |R_i,c_i;\mathcal {E})\right \rbrace d(\log |\Delta v_k^{\prime }|) \nonumber \\ && \end{eqnarray} \tag{ 13 }$$

where $\mathcal {L}(\Delta v_1,\cdots,\Delta v_{n-1}|R_i,c_i;\mathcal {E})$ is the likelihood in Equation (12). This is our binary likelihood expressed in terms of the error model parameters, which must be calculated by performing the integral every time the error model parameters are varied. By now it should be evident why only one long time interval (>10 days) between measurements is allowed: multiple long time intervals would require multiple integrations to account for binary variability over each interval, and this would be prohibitively expensive to compute while varying the error parameters.

Finally, for each star with a time interval longer than a week, we express the likelihood in terms of the binary fraction B as

$$\begin{eqnarray} \mathcal {L}(\Delta v_i|R_i&,&c_i,\Delta t_i,M;B,\mathcal {E}) = (1-B)\mathcal {L}_{nb}(\Delta v_i|R_i,c_i;\mathcal {E}) \nonumber\\ &&+\, B\mathcal {L}_b(\Delta v_i|R_i,c_i,\Delta t_i,M;\mathcal {E}). \end{eqnarray} \tag{ 14 }$$

We shall therefore constrain the four error model parameters $\mathcal {E}$ and binary fraction B simultaneously. We emphasize again that once the measurement errors have been determined by this procedure, we shall go back and find more precise binary constraints by including all the measurements. The maximum likelihood error model parameters obtained by this procedure will be used to determine the measurement errors for each star before proceeding with the full binary analysis.

We now apply the error model described in the previous section to derive measurement errors in the Carina, Fornax, Sculptor, and Sextans dSph galaxies. Before plunging into the analysis, we note that there are potential sources of non-Gaussian velocity variability besides binary orbital motion. The most common of these are due to false peak selection in the spectrum cross-correlation function, which occurs at low Tonry–Davis R-values and low S/Ns. In Walker et al. (2009a), measurements with R-values less than 4 are discarded for this reason. We take the additional step of discarding measurements with low S/N, which can also exhibit non-Gaussian error—this is apparent from the fact that the distribution of velocity changes Δv is noticeably non-Gaussian with the inclusion of low S/N measurements, even over timescales too small (∼ a few days) for binary velocity changes to be observed. On the other hand, the derived measurement errors are not significantly affected for S/N thresholds greater than ∼1, except for cuts well above 2 when the sample size is reduced significantly. For this reason, in the following analyses we shall discard measurements with S/N less than 1.2. While this S/N threshold might appear somewhat low, it should be borne in mind that the original R-value cut in Walker et al. (2009a) already removed many low S/N measurements with poor cross-correlation fits, since these are generally correlated with low R-values. However, there are a large number of measurements with S/N between 1 and 2 that have high R-values (up to ∼30), and these are very unlikely to be spurious measurements due to noise. Thus, we find that a S/N cut at 1.2 is sufficient to eliminate most non-Gaussian error without sacrificing an unnecessarily large number of genuine measurements. We do find, however, that measurements using the red channel with R-values less than 7 correlate strongly with low S/N (≲1), due to the red channel's lower spectral resolution compared to the blue channel. Therefore, as a further precaution these measurements in the red channel are also discarded to avoid non-Gaussian errors.

When generating the binary likelihoods, the star's magnitude is used to infer its mass and stellar radius under the assumption that it lies on the red giant branch (see Minor et al. 2010 for details); thus it is essential to exclude horizontal branch stars from each sample. In the Carina and Sculptor samples we discard stars with magnitudes m_V > 20.3 to ensure that horizontal branch stars are not included. The Fornax and Sculptor samples do not extend to sufficiently faint magnitudes to include horizontal branch stars, so no magnitude cuts are necessary in those samples.

As outlined in Section 4, our error model presently allows only one "long" time interval (>10 days) between measurements—any further measurements occurring at long time intervals are discarded, although they will be included later when deriving detailed binary constraints in Sections 7 and 8. After making the aforementioned S/N, R-value, velocity, and magnitude cuts in the Fornax sample, there remain 458 stars with repeat measurements, the largest of any galaxy in the survey. Of these, however, only 201 stars have time intervals of 1 yr or longer between any pair of measurements. Many of the latter subset have measurements at three epochs, with time intervals of approximately a week and one or two years between measurements. Since none of the stars have multiple time intervals greater than 10 days, no further measurements were discarded for the purpose of determining measurement errors in Fornax.

After applying our error model to the data, the results for the Fornax dwarf galaxy are plotted in Figure 1(a). For all repeat measurements in the sample, we plot the cumulative distribution of Δv/σ_2e, where Δv = |v₂ − v₁| is the difference between two successive line-of-sight velocity measurements and σ_2e is the measurement error in Δv, given by $\sigma _{2e}^2 = \sigma _1^2 + \sigma _2^2$ where σ₁, σ₂ are the individual errors on each velocity measurement. If the velocity changes are only due to Gaussian measurement errors, in the limit of a large sample this distribution should approach the complimentary error function ${\rm erfc}[\Delta v/\sqrt{2}\sigma _{2e}]$ (solid dark line).

Zoom In Zoom Out Reset image size

From the figure it is evident that Gaussian measurement errors are well-reproduced in Fornax for Δv ≲ 1.2σ_2e, comprising ∼80% of the stars in the sample, while binary orbital motion associated with large velocity changes becomes manifest above this threshold. The maximum likelihood binary fraction in Fornax is B  ≈  0.5, while the maximum likelihood error parameters are $\mathcal {E} = \lbrace \alpha _{{\rm red}}=36.9$, α_blue = 55.2, x_red = 1.06 and x_blue = 1.37}. To test this, we simulated a dataset with the same derived measurement errors, epochs, and magnitudes as in the real data, and assigned new velocity values for each measurement assuming a star population with binary fraction B = 0.5. The distribution for a typical realization is plotted in Figure 1(a) (dashed line), and is quite consistent with the data given the scatter in samples of this size. By comparison, we plot the distribution using the published measurement errors from Walker et al. (2009a), from which it is evident that the measurement errors are underestimated even for small velocity variations. This discrepancy may be not only due to binaries, but also the presence of low S/N measurements biasing the error model. We thus conclude that the published measurement errors for the Fornax dSph are significantly underestimated; in fact the median measurement error in the sample using the published errors is 1.1 km s⁻¹, while the median measurement error from our model is 1.7 km s⁻¹, differing by a factor of ≈55%. The median measurement errors in each galaxy are listed in Table 1.

In the case of the Carina and Sextans samples, contamination by foreground Milky Way stars is a major concern because of the relatively low surface brightness of these galaxies. In addition, many stars selected for repeat measurements have velocities more than 3σ away from the galaxy's systemic velocity, where σ is the galaxy's velocity dispersion. Since binary velocity variations greater than 50 km s⁻¹ are extremely unlikely, noting that each galaxy has a velocity dispersion of order ∼10 km s⁻¹, we therefore discard stars with velocities >60 km s⁻¹ away from the galaxy's systemic velocity as probable nonmember stars. (It should be noted that a small subset of these stars do indeed show velocity variations beyond the measurement error, but with an amplitude significantly smaller than 50 km s⁻¹. Given that nearly all have high metallicities, such stars are most likely Milky Way binaries.) In the Carina sample, there are 257 stars with repeat measurements; after making the rough velocity cut described above, together with a magnitude cut at m_V = 20.3 to get rid of horizontal branch stars, only 107 stars remain in the sample. Similarly, the Sextans sample contains 203 stars with repeat measurements, among which only 134 remain after making velocity and magnitude cuts. By contrast, the Sculptor sample contains 198 stars with repeat measurements, of which 190 remain after velocity and magnitude cuts, indicating relatively little contamination by nonmember stars.

Applying our error model to the remaining sample in the Carina dSph, we obtain an extraordinarily low best-fit binary fraction, B  =  0.05, the value of which will be refined in Section 7. Using our derived measurement errors, the resulting distributions in Δv/σ_2e are plotted in Figure 1(b). Because of the relatively small size of the sample, the distribution can differ significantly from the expected large-sample limit, but the absence of binary velocity variation is immediately apparent in the figure (we will return to this in detail in Section 7). To test our model, we again simulate a dataset with the same derived measurement errors, epochs, and magnitudes as the real data, and assign new velocity values for each measurement assuming our best-fit binary fraction for Carina. Two random realizations are plotted; although the distributions can vary considerably because of the small sample size, we see that the actual distribution using our measurement errors is indeed consistent with the model.

Finally, we apply our error model to the Sculptor and Sextans dSph samples. After discarding measurements in the red channel with low R-values, we find that the remaining Sextans sample contains only 16 stars with repeat measurements in the red channel, considerably less than the other galaxies in the sample. As this is an insufficiently small number of repeat measurements to derive reliable measurement errors for the red channel, we discard the red channel measurements in Sextans altogether for the remainder of this work. After making the final cuts, we plot distributions in Δv/σ_2e for Sculptor and Sextans in Figures 1(c) and (d) respectively. For each galaxy, we plot distributions using our derived measurement errors, along with distributions using the published measurement errors. From the figure it is evident that for both galaxies, the published data underestimates the measurement error for large velocity differences, while in Sextans the error is overestimated for small velocity changes. By contrast, the Gaussian error is well-accounted for in our derived measurement errors, with binary velocity variation showing up for velocity changes greater than ≈1.4 km s⁻¹.

While the principle goal of our analysis in this section has been to derive measurement errors, a rough estimate of the binary fraction in each galaxy has been derived simultaneously. However, the binary constraints obtained in this section are only approximately valid, for several reasons: first, to save time during the likelihood calculation, for each star with a long time interval (>10 days) between measurements, the long time interval was rounded to a multiple of 1 yr—thus the binary likelihoods were generated for multiples of 1 yr without the need for additional interpolation. For some stars the rounding error is as long as three months, so the binary likelihood is only approximately accurate in those cases. More importantly, as already discussed, certain repeat measurements were discarded if the star already contained multi-epoch measurements with one long time interval. Some of those discarded measurements may indicate binary velocity variation and thus are important for finding binary constraints. Finally, by inferring the binary fraction, we have assumed thus far that the distribution of orbital periods in each galaxy is identical to that of Milky Way field binaries, which may not be the case. Ultimately, it would be preferable to constrain not only the binary fraction, but also the period distribution parameters in each galaxy. The full binary analysis will be performed in Sections 7 and 8.

To verify that we have reproduced the Gaussian measurement errors well, we shall apply our error model to a series of simulated datasets with a distribution of Gaussian measurement errors similar to that obtained in the Magellan/MMFS survey of Walker et al. (2009a). This is accomplished by drawing a random Tonry–Davis R-value for each measurement from a Gaussian distribution whose mean and spread are similar to those observed in the real data, while the spectrograph channel used to obtain each measurement is chosen to be either red or blue with equal probability. The "true" error model parameters are chosen to be identical to those obtained from the Fornax sample (which have been obtained in Section 5). First we consider samples with measurements at two epochs, separated a year apart, with sample sizes of 200 and 500 stars. In each case, we simulate a sample with a binary fraction of zero and apply our error model without correcting for binaries. We also simulate a sample with a binary fraction B = 0.5 and apply our error model with the binary correction. In all simulations, an absolute magnitude limit M_V = 1 is assumed.

The results are displayed in Figure 2. As in the previous section, we plot the cumulative distribution of Δv/σ_2e, where Δv = |v₂ − v₁| is the difference between two successive line-of-sight velocity measurements and $\sigma _{2e} = \sigma _1^2 + \sigma _2^2$ where σ₁, σ₂ are the individual errors on each velocity measurement. In Figures 2(a) and (b) we consider two typical random realizations of a 200-star sample, where each star has measurements at two epochs separated by a year, and plot distributions using the "true" measurement errors for each simulation, as well as the derived measurement errors by applying our error model to the simulation. In each realization, the distribution using the derived errors follows the distribution using the "true" errors for both the B = 0 and B = 0.5 populations; the percent difference between the derived and true errors is smaller than 10% for at least ≈80% of the stars in the sample, and smaller than 30% for ≈95% of the stars in the sample. The Gaussian measurement errors are thus reasonably well-reproduced by the error model for most stars in the sample. For velocity changes Δv/σ_2e < 1, the distributions follow the large-sample limit, indicating the distribution of velocities below this threshold is Gaussian and binary behavior is not distinguishable from measurement error for such small velocity changes. As can be seen from Figures 2(c) and (d), the measurement errors are very well-reproduced in a 500-star sample regardless of the realization—in both cases shown, the percent difference between the derived and true errors is smaller than 10% for at least 97% of the stars in the sample. Since the distributions are cumulative, if a few measurement errors are reproduced incorrectly, the distribution using the derived errors can be offset from the true distribution; note however that in Figure 2(c) the distribution for B = 0.5 parallels the true distribution, indicating most of the measurement errors are indeed accurate.

Zoom In Zoom Out Reset image size

As a final sanity check, we verify that our error model works for a larger number of epochs. In Figure 3 we consider a sample of 200 stars with measurements at four epochs, where the first time interval between measurements is equal to 1 yr, and the remaining time intervals is equal to 1 day. Since time intervals shorter than a week are too short for binary orbital motion to be evident, we expect the additional epochs to aid in constraining the measurement errors, although without helping to constrain the binary fraction directly. We find that applying our error model to the four-epoch sample reproduces the measurement errors very well, although the sample with B = 0.5 is not as well reproduced as in the 500-star sample because the binary fraction is not as well constrained.

Zoom In Zoom Out Reset image size

7.1. Sample Selection Criteria; Contamination by Milky Way Stars

We now adopt a Bayesian approach to infer properties of the binary populations in each galaxy. For simplicity, in this section we will assume the distribution of binary periods in each galaxy is identical to that of Milky Way field binaries (with a mean period of ∼180 yr and log-spread of periods equal to 2.3; see Duquennoy & Mayor 1991), an assumption that will be relaxed in Section 8.

For reasons discussed in Section 5, we make a magnitude cut to exclude horizontal branch stars in the Carina and Sextans galaxies. Because these galaxies have relatively small member samples, however, it would be preferable to make a more specific color–magnitude cut to exclude the horizontal branch without sacrificing too many red giant branch stars. To accomplish this, we cut out stars with m_V > 20.4 for both galaxies, but make an additional color–magnitude cut to exclude the bright end of the horizontal branch. For Carina, stars with m_V > 20.35 and color V − I < 0.85 are probable red clump horizontal branch stars and are therefore excluded; likewise, in Sextans, a similar cut is made for m_V > 20.3 and V − I < 0.8. (For color–magnitude diagrams of the red giant branch/horizontal branch region in each galaxy, we refer the reader to Walker et al. 2007). We also discard measurements with low S/N to avoid non-Gaussian errors that may bias the inferred binary fraction. As we will show specifically for Carina in Section 7.3, our results are essentially unchanged for S/N thresholds between 1.2 and 2, except that the constraints become weaker for higher thresholds as a greater fraction of stars are being removed from the sample. Therefore, as in Section 5, we will discard measurements with S/N smaller than 1.2 in the calculations that follow. To eliminate obvious nonmember stars, we make a rough velocity cut as described in Section 5 for each galaxy. Finally, measurements that occur within 10 days of each other are averaged together (using Equations (2) and (3)), since binary velocity variation cannot be observed on shorter timescales given the measurement errors. Only stars with repeat measurements at two or more epochs (after this averaging) are included in the final sample. Color–magnitude diagrams for our final sample in each galaxy are shown in Figure 4.

Zoom In Zoom Out Reset image size

Another potential bias in the inferred binary fraction is contamination by foreground Milky Way binary stars. The majority of foreground stars are K-dwarfs in the Milky Way disk, many of which can be expected to have binary companions. We can expect larger velocity variations in these binaries compared to binaries in the background dSph galaxy, for two reasons: first, since the observed nonmember star lies on the main sequence, the secondary star may not be significantly dimmer than the primary star. Thus the measured radial velocity may be that of the less massive secondary star, for which larger velocity amplitudes are expected. Second, since the observed member stars lie on the red giant branch, very close binaries with high velocity amplitudes will have been destroyed by Roche-lobe overflow (Paczyński 1971), whereas this is not necessarily the case for the foreground Milky Way stars since the observed star lies on the main sequence. Thus if there is significant contamination by Milky Way binaries, we can expect that it will most likely bias the inferred binary fraction to a high value.

We can cut down on contamination significantly by performing a rough velocity cut as outlined in Section 5. In addition, since the metallicity distribution is well-determined in each galaxy in terms of the Mg i triplet (λ ∼ 5170 Å) absorption line strength, we can use metallicity information to help distinguish between member and non-member stars. Following Walker et al. (2009b), we model the distribution of Mg-triplet pseudo-equivalent widths (or magnesium strengths) for both the member and nonmember stars as a Gaussian, so that the member star population has mean magnesium strength $\bar{w}$ and dispersion in magnesium strength σ_w; likewise, the nonmember stars have similar parameters $\bar{w}_{{\rm MW}}$, σ_{w, MW}. The values of each of these parameters characterizing the metallicity distribution in each sample are determined in Walker et al. (2009b). For each individual star we calculate the average magnesium strength w by performing a weighted average over all measurements w_i, and likewise we calculate the average measurement error in magnesium strength, e_w. Assuming Gaussian measurement errors in w, we then define the likelihood in magnesium strength for member stars as

$$\begin{equation} \mathcal {L}_w(w_i|e_{w,i}) = \frac{1}{\sqrt{2\pi \left(e_w^2 + \sigma _w^2\right)}}\exp \left[\frac{-(w-\bar{w})^2}{2\left(e_w^2 + \sigma _w^2\right)}\right]. \end{equation} \tag{ 15 }$$

The corresponding likelihood for Milky Way stars is identical, except with the Milky Way magnesium strength distribution parameters $\bar{w}_{{\rm MW}}$, σ_{w, MW}. For the member stars, we use the same velocity likelihood as in Equation (6), while for the nonmember stars we use the nonbinary velocity likelihood in Equation (5). Assuming a fraction F of stars in the sample are members, our likelihood can be written

$$\begin{eqnarray} &&\mathcal {L}(\Delta v_i, w_i|e_{w,i},\sigma _i,t_i,M; B, \mathcal {P})\nonumber\\ & = & (1-F)\mathcal {L}_{{\rm MW}}(\Delta v_i|\sigma _i)\mathcal {L}_{w,{\rm MW}}(w_i|e_{w,i}) \nonumber\\ & & +\, F\mathcal {L}(\Delta v_i|\sigma _i,t_i,M;B,\mathcal {P})\mathcal {L}_w(w_i|e_{w,i}). \end{eqnarray} \tag{ 16 }$$

Since we are not modeling the binary population of the Milky Way stars, there is still the possibility that a binary nonmember star may be counted as a member if it exhibits a large velocity change. However, the stellar metallicities in the Sculptor, Sextans, and Carina galaxies are significantly lower than that of the nonmember stars, so the metallicity likelihoods will significantly reduce the impact of these stars on the inferred binary fraction. The Fornax dSph has a metallicity distribution similar to that of the nonmember stars, so the method outlined above cannot resolve nonmember contamination in Fornax; fortunately however, because of its high surface brightness, the fraction of nonmember stars in the Fornax sample is expected to be very low.

7.2. Binary Fraction Constraints

Adopting the likelihood in Equation (16), we use a nested sampling routine (Skilling 2004; Feroz et al. 2009) to obtain marginal posterior probability distributions in the binary fraction B and member fraction F in each galaxy, assuming a flat prior in each parameter. With the exception of Sextans, we find the inferred member fraction in each galaxy is greater than 0.9 to within 99% confidence limits, indicating our rough velocity cut (discussed in Section 5) successfully removed most of the Milky Way stars. In Sextans, however, the member fraction is constrained to be $0.75_{-0.05}^{+0.04}$, indicating contamination may still be an issue, although metallicity information can still distinguish between member and nonmember stars in most cases due to Sextans' low mean magnesium line strength ($\bar{w}\, {=}\, 0.36 \ {\rm \mathring{\rm{A}}}$, compared to $\bar{w}_{{\rm MW}} = 0.84 \ {\rm \mathring{\rm{A}}}$ for the foreground Milky Way stars). Since the member fraction in each galaxy is well-constrained by the metallicities in any case, the result is not strongly dependent on our prior in the member fraction F.

In Figure 5, we plot posterior probability distributions in the binary fraction B after marginalizing over the member fraction F. The most probable inferred binary fractions in each galaxy are listed in Table 2, where the errors are given by the 15.87% and 84.13% percentiles of the probability distributions (corresponding to 1σ error bars in a Gaussian distribution).

Zoom In Zoom Out Reset image size

Table 2. Best-fit Binary Fractions, Assuming Milky Way Period Distribution

Galaxy	Binary Fraction
Carina	$0.14^{+0.28}_{-0.05}$
Fornax	$0.44^{+0.26}_{-0.12}$
Sculptor	$0.59^{+0.24}_{-0.16}$
Sextans	$0.69^{+0.19}_{-0.23}$
Combined	$0.46^{+0.13}_{-0.09}$

Download table as:  ASCIITypeset image

To evaluate whether the inferred binary fractions in Figure 5 are consistent with that of Milky Way field binaries, we note that binaries with G-dwarf primaries in the solar neighborhood are found to have a binary fraction of ≈0.5 (Duquennoy & Mayor 1991); however, binary fraction is known to be a function of mass, with smaller-mass primaries having smaller binary fractions (Fischer & Marcy 1992; Raghavan et al. 2010). In the solar neighborhood, Raghavan et al. (2010) finds that primary masses in the approximate range of 0.75–1 M_☉ correspond to an average binary fraction of ∼0.4, while 1–1.4 M_☉ primaries correspond to an average binary fraction of ∼0.5. The galaxies in this study have multiple star populations with widely varying ages, and thus a range of primary masses on the red giant branch. While photometry shows the Sculptor and Sextans dSphs have predominantly very old star populations with ages ≳10 Gyr (de Boer et al. 2012a; Lee et al. 2009), the Fornax dSph shows a wide spread of stellar ages from 2 Gyr up to 13 Gyr (de Boer et al. 2012b), while Carina's dominant stellar population has ages of 4–7 Gyr (Hurley-Keller et al. 1998). On the red giant branch, this corresponds to stellar masses ranging from 0.8 M_☉ up to approximately 1.2 M_☉. We therefore assume an inferred binary fraction in the range 0.4–0.6 to be roughly consistent with that of Milky Way field binaries.

By this standard, Fornax, Sculptor, and Sextans all have binary fractions consistent with Milky Way field binaries to within 68% confidence limits. Carina, by contrast, has a most probable inferred binary fraction of 0.14, and its inferred binary fraction is less than 0.29 to within 68% confidence limits, and less than 0.57 to within 95% confidence limits. It must be borne in mind that here we have assumed Carina's period distribution to be similar to that of the Milky Way. More generally, the results imply that Carina is largely devoid of binaries with periods (≲10 yr) that would produce observable velocity changes over the relevant timescale of a few years. As will be made explicit in Section 8, this could also be explained if the mean period of Carina's binary population is longer than that of the Milky Way, or has a smaller spread of periods, rather than having a very small binary fraction. Whatever the reason, Carina's lack of short-period binaries appears to be inconsistent with having a binary population similar to that of Milky Way field binaries. We will discuss the possible implications of this result for star formation in Section 10.

Continuing under the assumption of Milky Way-like period distributions in each galaxy, we can ask, what is the binary fraction of the combined sample of all four galaxies? There are 621 stars in the combined sample, affording much better constraints on binary fraction compared to the sample in each galaxy separately. The inferred binary fraction of the combined sample is plotted in Figure 6 (solid curve), and its best-fit value is given in the last entry of Table 2. The most probable inferred binary fraction is 0.46, and the 68% confidence interval lies in the range 0.37–0.59, which is consistent with the Milky Way field binary fraction ≈0.5 for solar-type stars. Finally, since Carina has an anomalously low apparent binary fraction, we also plot the combined sample of Fornax, Sculptor, and Sextans without including the stars from Carina (dotted curve). The most inferred probable binary fraction without Carina is 0.55, with the 68% confidence interval lying in the range 0.44–0.70, somewhat high but still consistent with Milky Way field binaries.

Zoom In Zoom Out Reset image size

7.3. Low Binary Fraction in Carina dSph?

The apparent shortage of short-period binaries in the Carina dSph cannot be easily explained by sources of bias in the inferred binary fraction. Systematic errors such as contamination by foreground Milky Way stars, non-Gaussian errors due to low S/N or other factors, or the presence of variable stars all tend to increase the non-Gaussian velocity variability of the stars, and thus would inflate the inferred binary fraction, rather than reduce it. Indeed, such biases might be a factor in the somewhat high inferred binary fractions of Sextans and Sculptor. In particular, contamination by Milky Way binaries could be inflating the binary fraction of Sextans, which has a lower member fraction (≈0.75) in our selected sample compared to the other galaxies.

We can, however, investigate whether the Carina results are sensitive to the assumed S/N cutoff. In Figure 7 we plot posteriors in the binary fraction of Carina for S/N cutoffs ranging from 1 to 1.6. We find that increasing the cutoff from 1.2 does not change the inferred binary fraction significantly, although the constraint is weakened slightly because the sample size is reduced for higher cutoffs. If the S/N cutoff is reduced to 1, however, the most probable inferred binary fraction jumps from 0.14 to 0.24. This is due to the inclusion of a single star, labeled Car-1543, with a velocity change Δv = 15.7 km s⁻¹ and Δv/σ_2e = 3.6. While it may appear surprising that a single star in a sample of ∼100 stars could affect the best-fit binary fraction so dramatically, one must keep in mind that only a small fraction of the stars exhibit velocity changes beyond the measurement error. This accounts for the Poisson-like probability distribution in the binary fraction evident in Figure 7. If Car-1543 is indeed a binary, then the inferred binary fraction is less than 0.43 at 68% CL, and less than 0.74 at 95% CL; this increases somewhat the probability that Carina's binary population is consistent with Milky Way field binaries, although the best-fit binary fraction is still low. However, the radial velocity measurements for this star, 206.6 km s⁻¹ and 222.3 km s⁻¹ (compared to Carina's systemic radial velocity of 223 km s⁻¹), each have relatively low S/N, equal to 1.09 and 1.19 respectively, and thus may be suspect.

Zoom In Zoom Out Reset image size

A further complication arises in that there are seven stars in the Carina sample at magnitudes m_V ≳ 20.3 (placing them on or near the horizontal branch) that show significant velocity variations, greater than 15 km s⁻¹ and as high as 40 km s⁻¹. Only one of these stars lies within the color–magnitude cut mentioned above, and this star is a clear non-member in view of its low velocity compared to Carina's systemic velocity. However, although none of these stars fall within our sample, such large velocity variations at faint magnitudes call for an explanation. None of the stars have colors that place them in the RR Lyrae instability strip, so they are unlikely to be variable stars. Instead, four of them lie in the red clump horizontal branch region of the color–magnitude diagram, two are clear Milky Way stars due to their low velocities, and only one member star lies on the red giant branch. However, three stars have S/N less than 1.1, and all seven stars have measured velocities with low Tonry–Davis R-values, between 4 and 5 (measurements with R-values less than 4 were discarded from the published sample to avoid non-Gaussian error due to false-peak selection; see Walker et al. 2009a). Since the velocity variations are quite high compared to what might be expected from binaries, we find it likely that these stars show velocity variations because of non-Gaussian error related to their relatively faint magnitudes. While the other galaxies in the sample contain relatively few stars showing variability at these magnitudes, we note that the Carina sample contains 60 multi-epoch stars with m_V > 20.3, while the other three galaxies combined only have 19 multi-epoch stars in this magnitude range. Thus it is perhaps not surprising that Carina has a disproportionately high number of stars with large velocity variations at the faint magnitude end of the sample. We consider it unlikely that the velocity variations of stars in Carina beyond m_V > 20.3 are primarily due to binary motion.

In the previous section we derived binary fraction constraints by assuming the period distribution of each galaxy is identical to that of Milky Way field binaries, which takes the form of a log-normal distribution with mean μ_log P = 2.24 and dispersion σ_log P = 2.3. Here we shall likewise assume the period distribution to have a log-normal form, but will allow μ_log P and σ_log P to vary. Thus the likelihood in each star is given by Equation (6), where our set of binary parameters $\mathcal {P} = \lbrace \mu _{\log P},\sigma _{\log P}\rbrace$. As in Section 7.2, we also use magnesium strength information to help determine membership and constrain the member fraction F. Thus our four parameters being varied are the member fraction F and the set of binary parameters {B, μ_log P, σ_log P}.

As in the previous section, we choose a flat prior in the binary fraction B over the interval from 0 to 1. Choosing a sensible prior in the period distribution parameters μ_log P and σ_log P, however, requires more careful examination. Binary populations in nearby open clusters have been shown to display period distributions that are more peaked compared to that of field binaries (Brandner & Koehler 1998; Scally et al. 1999). This suggests that the observed wider period distribution of Milky Way field binaries may represent a superposition of narrower distributions resulting from a variety of initial star-forming conditions. However, the period distributions observed in open clusters still cover multiple decades of period, and that of classical dSphs might be expected to be wider still. To be conservative, we shall consider a uniform prior in the mean log-period μ_log P over the range [0, 4], where the lower limit in this range corresponds to a mean period of one year; as we shall see shortly, a binary population with such a short mean period would produce velocity variations in excess of what is observed in each galaxy, unless the binary fraction is quite small. We will also choose a conservative prior over the period distribution width σ_log P over the range [1, 3]. The lower limit σ_log P = 1.0 is narrower than any open cluster yet observed, and the field binary populations in larger dwarf galaxies might be expected to be broader than that of clusters; on the other hand, given the relatively small size of dwarf galaxies, it would seem unlikely for the period distribution to be considerably broader than that of the Milky Way. In practice, the range of the priors will only be relevant when marginalizing over the allowed range of σ_log P to obtain probability distributions in B and μ_log P, as a more conservative prior would produce weaker constraints in these parameters.

When attempting to constrain the period distribution of a binary population, it becomes desirable to include as many repeat measurements as possible for each star. Unfortunately however, calculating the binary likelihoods becomes computationally expensive for more than two epochs, because the binary tail in the likelihood becomes spread out over the velocity space (Δv₁, Δv₂,...) and thus a large number of Monte Carlo points are required to generate a smooth likelihood, especially for stars with large velocity changes in which case the likelihood is quite small. By running our code in parallel over a large number of processors, we can generate likelihoods for stars with measurements over as many as four epochs. Only the Sextans sample contains stars with measurements at five epochs or more, albeit a relatively small number (seven stars, assuming the epochs to be separated by more than two days). For these stars, we discard any epochs occurring after the fourth measurement.

For each galaxy, we calculate the likelihoods for each star in a table of values over the period distribution parameters and subsequently interpolate to find the binary likelihood for any combination of binary parameters μ_log P, and σ_log P. After generating the likelihoods for each star, we obtain joint posterior probability distributions in μ_log P versus B, and μ_log P versus σ_log P by marginalizing over the other parameters. Contour plots of the resulting posteriors for each galaxy are shown in Figure 8. For simplicity, each distribution is normalized so the highest peak value is equal to 1. For each galaxy, we find that the resulting posteriors are degenerate in the three binary parameters—for example, a population with a low binary fraction and short mean period, and a population with a high binary fraction and long mean period, are equally allowed by the data. These degeneracies result from the fact that a large fraction of short-period binaries can be obtained by having a large binary fraction, a short mean period (small μ_log P), or a wide period distribution (large σ_log P). The degeneracy between the binary fraction and the period distribution parameters has been investigated in detail in Minor et al. (2010). In the next section we shall discuss what would be required to break this degeneracy and obtain strong independent constraints on the binary fraction and period distribution parameters.

Zoom In Zoom Out Reset image size

By way of comparison, Milky Way field binaries in the solar neighborhood have binary parameters B = 0.5, μ_log P = 2.24, σ_log P = 2.3. This set of parameters lies squarely within the allowed parameter region for the Fornax, Sculptor, and Sextans dSphs, indicating that a Milky Way-like binary population is compatible with the data in those galaxies. On the other hand, the Milky Way parameter set is an outlier for Carina, falling on the edge of its allowed parameter space. Compared to Milky Way field binaries, Carina's binary population likely has either a smaller binary fraction (B < 0.5), a longer mean period (μ_log P > 2.3), or a combination of these.

For each galaxy, we see that parameter regions where the mean period is short and the binary fraction is sufficiently high are ruled out, because such binary populations would produce greater velocity variations than are observed in the data. On the other hand, regions where the mean period is long and binary fraction is low are also ruled out, since they would produce fewer velocity variations than are observed. Once again the Carina dSph stands out, in that it allows small binary fractions even if the mean period is relatively long (up to 10,000 yr). This reflects the lack of velocity variations observed in Carina. However, in the region where the mean period is longer than that of the Milky Way (μ_log P > 2.24), the binary fraction is essentially unconstrained, for a simple reason: even if the binary fraction is high, binary motion can be "hidden" if the width of the period distribution is small, restricting all binaries to long periods; conversely, if the period distribution is broadened, the binary fraction can be reduced so that very few short-period binaries would be observed. Thus there is a degeneracy between binary fraction and the width of the period distribution that prevents us from constraining the binary parameters if the mean period is assumed to be similar to the Milky Way or longer. We can say, however, that if the mean period is of order 30 yr (μ_log P ≈ 1.5 or shorter, the binary fraction must be smaller than 0.4 regardless of the width of the period distribution, and is likely smaller than 0.2. We emphasize again that these results are somewhat prior-dependent, in that a more conservative prior on σ_log P would widen the allowed parameter region somewhat.

9.1. Degeneracy between Binary Fraction and Period Distribution Parameters

Since the likelihood method demonstrated in previous sections can constrain the parameters characterizing a binary population, it is natural to ask: how large a sample, and how many epochs, would be required to obtain strong and independent constraints on the binary fraction and period distribution parameters? To give some idea of this, we first simulate a binary population with a binary fraction B = 0.5, and period distribution parameters identical to that of Milky Way field binaries (μ_log P = 2.24, σ_log P = 2.3). We assume a measurement error of 1.5 km s⁻¹, which is attainable from multi-object spectrographs, especially if measurements are averaged together over short timescales. After using a Monte Carlo simulation to generate mock stellar line-of-sight velocities in a large sample of 1500 stars, we then obtain posterior probability distributions using the same procedure as in Section 8 with identical priors. Posteriors in B versus μ_log P and μ_log P versus σ_log P are plotted in Figure 9.

Zoom In Zoom Out Reset image size

Figure 9(a) shows that while the allowed parameter space in B and μ_log P has narrowed considerably compared to the data (Figure 8), there is a clear degeneracy between these parameters: a high binary fraction and long mean period can produce the same observed binary variation, to within the measurement errors, as a low binary fraction and short mean period. A similar degeneracy exists between the binary fraction and spread of periods σ_log P. The reason for this degeneracy, which is explored in detail in Minor et al. (2010), is that only binaries with periods in the short-period tail of the period distribution have observable variations beyond the measurement error. Thus the shape of the period distribution is not well constrained, as is evident in Figure 9(b), and adjusting the mean period or spread of periods can be largely compensated for by changing the binary fraction. We find that this degeneracy is still present even if measurements are taken at three or four epochs, indicating that the measurement error is simply too large to constrain the shape of the period distribution independently of binary fraction. We note, however, that if the mean period of a binary population is shorter than that of Milky Way binaries, then the period distribution can be constrained more easily since velocity variation is observable for binaries with periods covering a larger portion of the period distribution.

Apart from the aforementioned degeneracy in binary fraction and mean period, the most probable point in the parameter space of Figure 9(a) is at a low binary fraction and short mean period compared to that of the input population. Rather than being largely statistical in nature, this is actually a consequence of the prior chosen. We chose a flat prior in the width of the period distribution σ_log P extending to 3; this allows for periods as short as a few days. Binaries with such short periods would exhibit large velocity variations, and hence for large period distribution widths, the binary fraction would have to be quite small to be in accord with the (simulated) data. Thus, by marginalizing over our prior in σ_log P, we are integrating over regions of parameter space where the binary fraction must be small, which biases the binary fraction to a low value. The lesson is that one must be careful not to choose an overly conservative prior: by choosing a prior that allows for too large a period distribution width, one may be inadvertently biasing the results.

9.2. Breaking the Degeneracy between Binary Fraction and Period Distribution Parameters in a Simulated Dataset

Before we tackle the question of what would be required to constrain the binary fraction and period distribution parameters independently in dSphs, it would be reassuring to verify that this is even possible in principle using our likelihood approach. To show this, we simulate a binary population with a mean period of 10 yr (μ_log P = 1), a spread of periods σ_log P = 1.5, and a binary fraction B = 0.5. We choose a measurement error of 0.5 km s⁻¹ and generate simulated samples of 500, 1000, 2000, and 5000 stars with this binary population. Since binary periods are better constrained with a large number of epochs, we assume measurements taken at four epochs (the largest number allowable at present; computations become prohibitively expensive beyond four epochs), each spaced a year apart.

From the simulated line-of-sight velocity data, we obtain posterior probability distributions via a nested sampling routine in each parameter B, μ_log P, σ_log P, as shown in Figure 10. Even for a sample of 500 stars, the degeneracy between binary fraction B and mean log-period μ_log P is being broken, as evidenced by the peaks appearing in the posteriors near their correct values. The width of the period distribution σ_log P is more difficult to constrain, and only becomes well-constrained beyond 1000–2000 stars. The degeneracy between all three parameters is clearly broken beyond 2000 stars. Note that once the degeneracy is broken, the error in the maximum-probability values in each parameter (i.e., the difference between the peak value and the "true" value) is considerably smaller than the width of the posterior distributions, which suggests that the Bayesian error may be overly conservative in this case.

Zoom In Zoom Out Reset image size

How important is the number of epochs for constraining the period distribution parameters? In Figure 11 we plot posterior distributions in B for two-, three-, and four-epoch samples of 2000 and 5000 stars. If measurements are taken at only two epochs, the degeneracy is not broken even for a 5000-star sample. Likewise, if measurements are taken at three epochs, the degeneracy is not broken for a 2000-star sample, although it is broken for a sample as large as 5000 stars, as evidenced by the peak appearing in the posterior. However, the degeneracy is much easier to break in a four-epoch sample; the binary fraction is better constrained in a four-epoch sample of 2000 stars compared to a three-epoch sample of 5000 stars. We thus conclude that if the measurement error is small enough to allow the period distribution to be constrained, three- or four-epoch samples of at least several hundred (and possibly thousands) of stars are required to obtain strong and independent constraints on the binary parameters—in other words, to break the degeneracy between the binary fraction and period distribution parameters.

Zoom In Zoom Out Reset image size

9.3. Required Measurement Error for Obtaining Independent Constraints on Binary Parameters

For a binary population with a given mean period, what measurement error would be required to constrain the period distribution independently of the binary fraction? In principle, even if only binaries whose periods lie on the short-period tail of the period distribution have observable velocity changes, the shape of the period distribution could be inferred if a sufficient number of binaries are observed whose periods are well-known. This only holds to the extent that the assumption of a log-normal period distribution is accurate, so that the shape of the period distribution can be extrapolated to longer periods. Unfortunately however, with only two to four epochs, binary periods of individual stars cannot be determined with great accuracy. Furthermore, the likelihood approach is not computationally feasible at present with more than four epochs, because of the enormous number of Monte Carlo points required to generate smooth likelihoods. It is thus reasonable to expect that to constrain the period distribution well using the likelihood approach, observing only binaries on the short-period tail of the period distribution is not sufficient. We make the assumption that binary periods nearly as long as the mean period (≈180 yr for Milky Way field binaries) should be observable, i.e., should have observable velocity variations, in order to constrain the period distribution parameters well. We can get an approximate idea of the measurement error required for this by considering typical velocity amplitudes associated with binaries with a given orbital period. In Minor et al. (2010), a formula is derived for the velocity amplitude of a binary as a function of period, which is averaged over the possible orientations, mass ratios, and eccentricities of binary systems (assuming Milky Way-like probability density functions in these parameters):

$$\begin{equation} v_{{\rm max}} \approx (5.7 \ {\rm km\ s^{-1}}) \left(\frac{M/M_\odot }{P/{\rm year}}\right)^{\frac{1}{3}} \end{equation} \tag{ 17 }$$

where M is the mass of the primary star. If we assume M ≈ 0.8 M_☉ which is typical of old stellar populations in dwarf spheroidals, we find that binaries with periods longer than a few decades will yield velocity amplitudes less than 2 km s⁻¹, and thus would be difficult to distinguish from a measurement error of 1–2 km s⁻¹. Binaries with a period of 180 yr would have velocity amplitudes of ≈0.9 km s⁻¹ on average. However, there is a further difficulty: if measurements are taken over the course of a few years, binaries of such long periods would traverse only a small fraction of their orbits and thus would exhibit velocity changes much smaller than their overall velocity amplitude. As a rough approximation, if we assume a circular orbit, it can be shown that the rms change in line-of-sight velocity after a time Δt is given by the simple equation

$$\begin{equation} \Delta v_{{\rm rms}} = v_{{\rm max}} \left|\sin \frac{2\pi \Delta t}{P}\right|, \end{equation} \tag{ 18 }$$

where we have averaged over all possible orbital angles. In an approximate sense, velocity changes will be observable if they are larger than the measurement error between two pairs of velocity measurements. If the measurement error σ_m is roughly the same for either measurement, then the error in the velocity difference Δv is $\sigma _{2e} = \sqrt{2}\sigma _m$. Using this together with Equations (17) and (18), we find the measurement error must be smaller than

$$\begin{equation} \sigma _m \approx 4.0 \ {\rm km}\ {\rm s}^{-1}\left(\frac{M}{P}\right)^{\frac{1}{3}}\left|\sin \frac{2\pi \Delta t}{P}\right|. \end{equation} \tag{ 19 }$$

This is the maximum (approximate) measurement error with which one can observe binary motion for a binary with period P over a timescale Δt. In Figure 12 we plot the required measurement error as a function of log of period over timescales of 1 yr and 5 yr. Note that the required measurement error is very small near periods for which the timescale Δt is an integer multiple of one-half the period, since the velocity change is typically small for these timescales. However if the period is long enough, the required measurement error decreases monotonically. For binaries with periods equal to the Milky Way's mean period (log P = 2.24), however, the maximum measurement error for observing binary motion over 5 yr is ∼0.1 km s⁻¹, and much smaller over a timescale of 1 yr. We can therefore expect that if a binary population has a mean period similar to that of Milky Way field binaries, it would be extremely difficult to constrain its period distribution parameters independently of binary fraction if the measurement error is significantly larger than 0.1 km s⁻¹, unless the measurements are taken over timescales considerably longer than 5 yr. The prospects are improved if the mean period is shorter; for example if the mean period is 10 yr, Figure 12 suggests that the period distribution can be constrained independently of binary fraction if the measurement error is less than roughly 1 km s⁻¹.

Zoom In Zoom Out Reset image size

It should be emphasized that Equation (19) gives only an approximate idea of the measurement error required to break the degeneracy between binary fraction and the period distribution parameters; it says nothing about the number of stars and epochs that would be required to break this degeneracy, even if the measurement error is small enough. As we explored in Section 9.2 for the case of a relatively short mean period, at least three or four epochs and several hundred up to several thousands of stars would be required, but this requirement would likely vary depending on the actual period distribution and measurement errors in question.

Finally, we note that the task of constraining the binary period distribution would be considerably easier if the binary fraction were known a priori. This is possible in globular clusters, where the binary fraction can be constrained by photometry (Sollima et al. 2007, 2012). Because the stars in the cluster have nearly identical ages and metallicities, the main sequence is well-defined and binaries will lie off the main-sequence for single stars due to the combined light of the primary and secondary star. After the binary fraction is constrained from photometry, this constraint could be included as a prior in a likelihood analysis following a radial velocity variability survey to constrain the period distribution in the cluster. Unfortunately in dwarf spheroidal galaxies, constraining the binary fraction through photometry is a much more difficult task due to the complicated distribution of stellar ages and metallicities in each galaxy.

We demonstrated in Section 7.2 that the Carina dSph sample is nearly devoid of short-period binaries, with a binary fraction less than 0.5 to within 90% confidence limits if a Milky Way-like period distribution is assumed. While the reason for Carina's apparent shortage of close binaries is unclear, several factors suggest themselves. Photometric studies show Carina to have a bursty star formation history (Smecker-Hane et al. 1996; Hurley-Keller et al. 1998), with two predominant star populations: approximately 25% of its stars are older than 10 Gyr, while the remaining majority are intermediate age stars with ages in the range of 4–7 Gyr. The question is therefore, what initial conditions during its relatively short period of star formation might have prevented close binaries from forming? Potential factors are the initial temperature of the star-forming clouds prior to collapse (Sterzik et al. 2003), as well as the presence of magnetic fields and strong radiative feedback, which can all play a role in inhibiting binary star formation (Price & Bate 2010). While Carina's low metallicity compared to the Milky Way might conceivably affect radiative feedback and hence suppress binary star formation, we note that the Sculptor and Sextans dSphs also have relatively low mean metallicities, and their best-fit binary fractions are both higher than 0.5. Thus, its low metallicity alone is unlikely to explain the apparent lack of binaries in the Carina dSph.

It is worth noting that a recent spectroscopic study of metal abundances in Carina (Lemasle et al. 2012) showed little evidence of iron enrichment from type Ia supernovae (SNe) during its intermediate period of star formation, even though a few Gyr should have been a long enough timescale for SNe Ia to contribute iron to the interstellar medium while star formation was still taking place. While there is still some uncertainty in the timescale for star formation of Carina's intermediate age population, this lack of evidence for SNe Ia would be expected if Carina is indeed deficient in short-period binaries, as the present study suggests.

As we cautioned in Section 7.3, however, the possibility that Carina's apparently low fraction of binaries may be a statistical coincidence increases somewhat if the star Car-1543 is a binary, although the best-fit binary fraction is still low. However, this star's observed velocity variation may be explained by low S/N measurements (S/N ∼ 1.1). Furthermore, stars near the horizontal branch region also exhibit velocity variations, although most or all of these are probably explained by poor spectrum cross-correlation fits related to their faint magnitudes.

We have developed a general methodology for constraining properties of binary star populations from multi-epoch line-of-sight velocity measurements in dwarf galaxies. This method has been applied to data from the Magellan/MMFS sample of Walker et al. (2009a) in the Carina, Fornax, Sculptor, and Sextans dwarf spheroidal galaxies to find constraints in their binary fraction and period distribution parameters. To obtain the best possible binary constraints, we have also re-derived the measurement errors in the sample by extending the error model of Walker et al. (2007) to account for binary orbital motion. The best-fit binary fractions in each galaxy, assuming a Milky Way-like period distribution, are listed in Table 2, with probability distributions in binary fraction plotted in Figure 5; more general probability distributions in the period distribution parameters and binary fraction are plotted in Figure 8. We conclude with the following points:

• 1.  
  If a Milky Way-like period distribution is assumed in each galaxy, the Fornax, Sculptor, and Sextans dSph galaxies have binary fractions that are roughly consistent with that of Milky Way field binaries, whereas the Carina dSph is apparently deficient in binaries compared to the Milky Way field. Carina's inferred binary fraction is less than 0.5 at the 90% confidence level, with a best fit value of $0.14^{+0.28}_{-0.05}$; thus a Milky Way-like binary fraction of ≈0.5 in Carina is statistically unlikely. Relaxing the assumption of a Milky Way-like period distribution, the lack of observed binary velocity variation in Carina could be explained by either a small binary fraction, or a long mean period, compared to the other galaxies in the sample (see Figure 8).
• 2.  
  With the exception of Carina, the published measurement errors in the Magellan/MMFS sample of Walker et al. (2009a) underestimate the Gaussian measurement error somewhat, as can be seen in Table 1. This is most striking in the case of the Fornax dSph, where the median measurement error is larger than the published value by a factor of ≈ 55%. While this difference of ≈0.6 km s⁻¹ in the median error is small compared to Fornax's velocity dispersion of ≈12 km s⁻¹, the extra measurement error should be taken into consideration when applying mass models to Fornax. This may be particularly of concern if higher moments in the velocity distribution are used, since stars lying on the tail of Fornax's velocity distribution may have unaccounted-for measurement error if the published errors are used.
• 3.  
  While multi-epoch surveys can in principle produce independent constraints on a galaxy's binary fraction and the period distribution parameters, these parameters are unfortunately degenerate given the present measurement errors. For example, a binary population with either a low binary fraction and short mean period, or a high binary fraction and long mean period, can both produce the same velocity variations observed in the data equally well. If a galaxy's mean binary period is roughly comparable to that of the Milky Way, strong independent constraints on the binary parameters can only be obtained with a measurement error of order ∼0.1 km s⁻¹ or smaller, obtainable only by a high-resolution spectrograph. In the near future, the best possible binary constraints may be obtained in globular clusters where the binary fraction can be estimated independently through photometry, followed by a multi-epoch radial velocity survey to constrain the binary period distribution.
• 4.  
  Although the degeneracy between the binary parameters cannot be broken with measurement errors presently obtainable by multi-object spectrographs, meaningful comparisons between different binary populations can still be made. This was demonstrated in Figure 8: although independent constraints on the binary fraction and mean period cannot be obtained, it is evident that the Carina dSph occupies a different region of parameter space compared to the other galaxies in the sample. Larger multi-epoch surveys in dwarf galaxies will narrow the parameter space further to allow even more robust comparisons between galaxies. By comparing the allowed parameter space for different galaxies, we can ultimately address the question of whether field binary populations are universal in nature, or vary depending on the initial conditions present during the galaxy's epoch of star formation. Large multi-epoch surveys, in concert with statistical methods like that demonstrated here, will provide a powerful test of star formation theories in the future.

The author would like to thank Matt Walker for generously providing the error model data without which this work would not have been possible. I would also like to thank Erik Tollerud for providing useful comments on the manuscript, and Manoj Kaplinghat, James Bullock, Erik Tollerud, and Greg Martinez for providing valuable feedback and support throughout the course of this project.

This research was supported in part by a grant of computer time from the City University of New York High Performance Computing Center under NSF Grants CNS-0855217, CNS-0958379 and ACI-1126113. This work was supported in part by NSF grant AST-1153335.