Search
Close
Search
Advanced Search
Search Menu
AI Discovery Assistant

Abstract

We have developed a novel direct algorithm to derive the mass-ratio distribution (MRD) of short-period binaries from an observed sample of single-lined spectroscopic binaries (SB1s). The algorithm considers a class of parametrized MRDs and finds the set of parameters that best fits the observed sample. The algorithm consists of four parts. First, we define a new observable, the modified mass function, that can be calculated for each binary in the sample. We show that the distribution of the modified mass function follows the shape of the underlying MRD, making it more advantageous than the previously used mass function, reduced mass function or reduced mass function logarithm. Second, we derive the likelihood of the sample of modified mass functions to be observed given an assumed MRD. A Markov chain Monte Carlo search enables the algorithm to find the parameters that best fit the observations. Third, we propose that the unknown MRD is expressed by a linear combination of a basis of functions that spans the possible MRDs, and we suggest two such bases. Fourth, we show how to account for the undetected systems that have a radial-velocity amplitude below a certain threshold. Without the correction, this observational bias suppresses the derived MRD for low mass ratios. Numerous simulations show that the algorithm works well with either of the two suggested bases. The four parts of the algorithm are independent, but the combination of these four parts makes the algorithm highly effective at deriving the MRD of the binary population.

methods: data analysis, methods: statistical, binaries: spectroscopic, stars: statistics

1 INTRODUCTION

The study of the mass-ratio distribution (MRD) of binaries, short-period binaries in particular, has a long history. This is because the MRD plays a key role in various aspects of the theory of binary formation and evolution. First, it provides one of the very few ways to challenge the theories of binary formation using observations (e.g. Bate & Bonnell 1997; Satsuka et al. 2017). Secondly, the primordial MRD is a one of the few major inputs for population syntheses of binaries (e.g. Toonen, Nelemans & Portegies Zwart 2012; Yungelson & Kuranov 2017), which try to predict, for example, the rate of supernova explosion and black hole mergers. Thirdly, the binary fraction and the MRD have been shown to play important roles in star cluster evolution (e.g. Hut et al. 1992; Benacquista & Downing 2013). Finally, an understanding of the low end of the MRD is crucial for the determination of the borders of the brown dwarf desert (e.g. Mazeh et al. 2003; Grether & Lineweaver 2006) that separates exoplanets from low-mass stellar secondaries.

It is therefore not surprising that quite a few studies have tried to derive the MRD of binaries, using short-period spectroscopic binaries (SBs) in particular. Reviews of the early studies can be found in Trimble (1990) and Mazeh & Goldberg (1992). Because some of these studies (e.g. Lucy & Ricco 1979; Duquennoy & Mayor 1991; Tokovinin 2000; Goldberg, Mazeh & Latham 2003; Fisher, Schröder & Smith 2005; Raghavan et al. 2010; Boffin 2010, 2015; Curé et al. 2015) yielded conflicting results, the shape of the MRD of short-period binaries continues to be an open question.

Derivations of the MRD should be based on a complete sample of binaries discovered by a systematic radial-velocity (RV) search. Generally, MRDs can depend on the mass of the primary star (Kouwenhoven et al. 2009), so analysed samples should be restricted to some narrow range of spectral types. In an ideal world, with spectra of unlimited resolution and signal-to-noise ratio, each SB would be a double-lined spectroscopic binary (SB2), with mass ratio derived directly from the ratio of the RV amplitudes of the two components. In reality, the derivation of the MRD of short-period binaries is based on samples for which most of the binaries observed are single-lined spectroscopic binaries (SB1s), where only the RVs of the primary star can be measured. Even after obtaining observations of a large sample of SB1s, the derivation of the MRD is hampered by the fact that for each of the SB1s, the orbital solution cannot yield the mass ratio itself but merely the mass function, which is a combination of two unknowns: the mass ratio and the plane-of-motion inclination angle. Therefore, a statistical approach must be applied to the observational results, assuming random distribution of the orbital inclination of the sample as a whole.

Two main approaches have been used to disentangle the MRD from the orbital inclination (e.g. Heacox 1995). In the inverse approach, we consider the sample of derived mass functions and work back, usually iteratively (see the classical work of Lucy 1974), to the underlying MRD of the binary population (e.g. Lucy & Ricco 1979; Mazeh & Goldberg 1992; Boffin 2010; Curé et al. 2015).

In the direct approach, however, we assume a certain MRD of the binary population, calculate the resulting expected distribution of some observable O, and compare it with the set of {Oi} obtained from the sample of SB1s, where the ith binary is represented by Oi. Then, we find the MRD that best fits the observed set {Oi} (e.g. Hogeveen 1992; Tokovinin 1992; Ducati, Penteado & Turcati 2011). The comparison between the expected distribution and the observed sample is usually done by comparing histograms, not a very powerful approach, which does not allow statistical derivation of the best parameter(s) of the MRD and its (their) confidence intervals.

In previous studies, the observable used was the reduced mass function y, obtained by dividing the mass function by the mass of the primary star; however, see Lucy & Ricco (1979) and Boffin (2010, 2015), who promoted the use of a logarithm of the observable instead. However, we will show that the expected distributions of both y and log y for very different MRDs are similar, making the derivation of the true MRD difficult.

Here we present a novel algorithm to solve the problem, which consists of four parts. First, we introduce a new observable S, which we call the ‘modified mass function’, that is derived for each SB1. We then show that the shape of the distribution of the obtained S for a binary sample is similar to that of the underlying MRD of the population. Even more importantly, different MRDs result in different S-distributions.

Secondly, we suggest a comparison of the expected distribution of S with the observed sample by deriving the likelihood of the observed set of {Si}, given the assumed MRD (see also Tokovinin 1992). We can then find the parameters of the MRD that maximize the likelihood of observing the sample, using a Markov chain Monte Carlo (MCMC) approach, for example.

Thirdly, we suggest that the unknown MRD is expressed by a linear combination of a basis of functions, with some unknown coefficients. We then search for the coefficients that maximize the likelihood of the observed sample.

Fourthly, following Mazeh & Goldberg (1992) we show how to account for the undetected binaries that have an RV amplitude below a certain threshold. Without the correction, this observational bias suppresses the derived MRD for low mass ratios.

In Section 2, we introduce the modified mass function. In Section 3, we detail the search for the best set of parameters using the MCMC process, and we suggest two sets of functions. In Section 4, we detail two simulated examples, which demonstrate that the algorithm works well. In Section 5, we present our correction function. In Section 6, we briefly summarize this work and lay out possible further refinements of the algorithm. In the following series of papers, we will apply the algorithm to various SB1 samples.

2 THE MODIFIED MASS FUNCTION

The mass ratio of a binary system is defined as q ≡ m2/m1, where m1 and m2 are the stellar masses of the primary and secondary, respectively. For an SB1, only the spectrum of m1 is seen in the spectrum of the system, and therefore only the primary RVs are obtained. When enough measurements are secured, an orbital solution is derived, yielding the orbital period P, eccentricity e and primary RV semi-amplitude K1. The binary mass function is defined as
(1)
where i is the orbital inclination. An estimate of m1 is often available from the binary spectra and can be factored out, leaving a reduced mass function y, with only two unknown parameters, q and i:
(2)
Henceforth, we assume that the primary is also the more massive star in the binary system, namely 0 < q ≤ 1. Under this assumption, the reduced mass function becomes bounded as well, 0 < y ≤ 0.25. The relation between q, y and i is plotted in Fig. 1. We choose to work in the (1 − cos i, q) plane, as the distribution of 1 − cos i is uniform for random orientation of the orbits. The plot shows the possible values of q for a binary with y = 10−2. The grey area in the plot shows all the possible cases with y ≤ 10−2.
The (1 − cos i, q) parameter space, with an SB1 with y = 10−2. The grey area shows the corresponding S value. The red dot marks the minimum value of q (i.e. $\mathcal {Q}_{y}$), which is 0.25.
Figure 1.

The (1 − cos i, q) parameter space, with an SB1 with y = 10−2. The grey area shows the corresponding S value. The red dot marks the minimum value of q (i.e. |$\mathcal {Q}_{y}$|⁠), which is 0.25.

Open in new tabDownload slide
Notably, each value of y is associated with a minimal possible q value, which can be determined by setting the inclination angle i in equation (2) to be 90°. This q minimum, |$\mathcal {Q}_{y}$|⁠, is the only real root of the polynomial Py(q),
(3)
for which an explicit expression was given by Heacox (1995) (see also equation A3). The value of 0.25 for y = 10−2 is plotted as a red dot in Fig. 1.

Previous techniques used the observables y or log y as tools for deriving the MRD. For the direct method, it is necessary to obtain the probability density function (PDF) of y, fy or flog y, given the PDF of the MRD, fq. To obtain fy(y; fq), we have to calculate the probability of obtaining a value between y and y + dy over the parameter space of Fig. 1, given fq(q). This is done in Appendix A.

We seek a transformation |$S=\mathbb {S}(y)$| such that significant functional properties of the MRD, fq, will be qualitatively demonstrated by its resulting PDF, fS. This results in three requirements. First, a uniform fq should yield a uniform fS. Secondly, the transformation |$\mathbb {S}$| is required to be strictly increasing and continuous. Finally, for fS to be comparable with fq, the range of |$\mathbb {S}$| is required to be the [0, 1] interval.

These requirements are uniquely met by the cumulative distribution function (CDF) of y, assuming a uniform distribution of q. For example, |$\mathbb {S}(y=10^{-2})$| is simply the grey area of Fig. 1 for y = 10−2. The area can be written as the integral
(4)
where the integrand is the height above the curve of Fig. 1. The relation between S, hereafter called the modified mass function, and y is demonstrated in Fig. 2.
Modified mass function S as a function of the reduced mass function y.
Figure 2.

Modified mass function S as a function of the reduced mass function y.

Open in new tabDownload slide

The modified mass function, by its definition, resembles a copula (e.g. Nelsen 2013). While copulas are widely used in many fields, especially in quantitative finance, its astrophysical applications are rare (e.g. Scherrer et al. 2010). A detailed derivation of |$\mathbb {S}$| appears in Appendix B.

2.1 Distribution of the modified mass function

To derive the PDF of S, we note that S is defined as a function of y, and therefore
(5)
where y(S) is the inverse of |$\mathbb {S}$|⁠. Because the modified mass function S is, by definition, the CDF of y for a uniform MRD, the last factor in equation (5) is
(6)
Therefore, we obtain
(7)
An explicit expression for fy is developed in Appendix A, and shown in equation (A7). Inserting the two expressions (i.e. the PDF of y for the actual MRD and for flat distribution), we finally obtain
(8)
where
(9)
(see Appendix A). Equation (8) is effectively a weighted average of fq for a given S value, taken over the allowed q range, |$[\mathcal {Q}_{y(S)},1]$|⁠, and weighted by the assumed distribution of isotropic inclination angles.

Fig. 3 shows the derived fy and fS for three different simple fq functions, demonstrating how, unlike fy or flog y, fS captures the shape of the underlying MRD, fq.

Three distributions of q and their corresponding y, log y and S distributions. The top panel shows the uniform (black line), linearly increasing (brown dashed line) and decreasing (green dotted line) distributions of q. The three lower panels show the corresponding y, log y and S distributions, with the same line types and colours.
Figure 3.

Three distributions of q and their corresponding y, log y and S distributions. The top panel shows the uniform (black line), linearly increasing (brown dashed line) and decreasing (green dotted line) distributions of q. The three lower panels show the corresponding y, log y and S distributions, with the same line types and colours.

Open in new tabDownload slide

3 DIRECT DERIVATION OF THE MRD

3.1 Likelihood derivation of the MRD

Let us assume that the MRD is characterized by a set of parameters |${\boldsymbol {c}}=\lbrace c_k\rbrace$|⁠. This could be, for example, a Gaussian distribution |$f_q\propto \exp [ (q-c_1)^2/2c_2^2 ]$|⁠, a power-law distribution |$f_q\propto q^{c_1}$|⁠, a flat distribution between c1 = qmin and c2 = qmax, and such like, or any combination of the above. Using equation (8), we transform the |$f_q(q;{\boldsymbol {c}})$| into |$f_S(S;{\boldsymbol {c}})$|⁠, which has the same set of parameters {ck}.

We wish to find the values of ck that best match the given sample of observed SB1s, with a set of {yi} that we transfer via equation (4) to the corresponding set of {Si}. The search is done by maximizing the log-likelihood of |${\boldsymbol {c}}$|⁠,
(10)
given {Si}. The core of the algorithm is the search for the best MRD in the S domain. Because the fitted ck values are shared by both fS and fq, it is clear that by fitting fS we can readily derive its underlying fq.
In practice, examples in this work were analysed with the emcee ensemble sampler (Goodman & Weare 2010; Foreman-Mackey et al. 2013). Each step in the generated chain yields a set of values for ck, from which the MRD, fq, is derived over a dense set of predetermined {qj}. The chain yields a posteriori distributions for each {fq(qj)}. We use these distributions to derive the median |$\lbrace \skew4\hat{f}_q(q_j)\rbrace$| and their 1σ confidence intervals |$\lbrace \hat{\delta }(q_j)\rbrace$| to finally yield
(11)

3.2 Expansion of MRD by a set of basis functions

Likelihood derivation of the MRD requires a predetermined functional model, whose parameters are searched to fit best the observed set of {Si}. It is therefore desirable to use models that can span a broad class of functions, thus avoiding a priori assumptions about the functional shape of the MRD. This can be achieved by approximating fq with a set of basis functions,
(12)
where ϕk(q) is the kth function and ck is its corresponding coefficient.
For each basis function, we derive its corresponding function in the S plane, denoted |$\tilde{\phi }_k(S)$|⁠, using equation (8),
(13)
The linearity of equation (8) links the expansion of fS to that of fq via the modified functions, namely
(14)
where ck is the kth coefficient of the fq series expansion. In this case, the parametrized probability density takes a simple form
(15)
where Mik is the value of |$\tilde{\phi }_k$| at Si. This can be written in a matrix form
(16)
where |$\boldsymbol {M}$| is the Mik design matrix.
The design matrix |$\boldsymbol {M}$| can be calculated given the sample {Si}, the basis {ϕk(q)} and its corresponding |$\lbrace \tilde{\phi }_k\rbrace$|⁠. The log-likelihood in terms of |$\boldsymbol {M}$| is
(17)
according to which the best |${\boldsymbol {c}}$| can be found.

In the following subsections, we suggest two possible sets of basis functions: the shifted Legendre polynomials and the boxcar functions. The two have complementary properties in terms of smoothness and locality. Obviously, other possibilities, such as harmonic functions or power series, can be considered and implemented in a similar manner.

3.2.1 Shifted Legendre polynomials

A possible basis is, for example, the shifted Legendre polynomials {Pk}
(18)
where the first four Pk(q), P0–P3, are plotted in Fig. 4, together with their corresponding modified functions |$\tilde{P}_k(S)$|⁠.
First four shifted Legendre polynomials (dashed black) and their corresponding modified functions (solid red).
Figure 4.

First four shifted Legendre polynomials (dashed black) and their corresponding modified functions (solid red).

Open in new tabDownload slide

The shifted Legendre polynomials have the property |$\int _0^1 P_k(x)\, {\rm d}x = \delta _{k0}$|⁠. This makes them suitable as a set of basis function for any PDF, as the integral of any combination of them over the range [0, 1] is unity, as long as c0 = 1.

3.2.2 Boxcar functions

Another basis is the set of boxcar functions, |$\lbrace \Pi _{_{N,k}}(x), k=1,\ldots ,N\rbrace$|⁠, that are simple unit pulses of the form
(19)
spanning the histogram-like models. The corresponding modified functions are derived through equation (8).

3.3 Starting point of the MCMC

In any MCMC search, it is important to start the chain not too far from the global maximum of the log-likelihood. Our starting point relies on the histogram of observed {Si}. For the boxcar basis, the starting point is taken as the normalized number of counts in the S histogram bins, whereas the starting point of the Legendre polynomial set was derived by a simple linear least-squares fit to the histogram bins (see Barlow 1989).

To choose the number of bins, Nbin, for the histogram we use the Rice rule (Terrell & Scott 1985),
(20)
where n is the size of the observed sample.

4 TESTING THE ALGORITHM

In order to test our algorithm, and the two bases presented above in particular, we ran numerous simulations, two of which are presented here. In each numerical experiment, we assumed an underlying MRD and prepared a simulated SB1 sample, drawing at random the values for the mass ratio and the inclination of each binary. We then derived the S value for each binary and applied our algorithm to the sample of modified mass functions twice, using at each time one of our two bases.

In all our simulations, we were able to retrieve the correct shape of the underlying MRD, with each of the two bases.

Here we present two simulations, one (Fig. 5) with an MRD composed of a fourth-degree polynomial, fq(q) ∝ 25(2q − 1)4 + 4, that peaks at q = 0 and q = 1, and the other (Fig. 6) composed of a Gaussian with a mean at q = 0.2 and a standard deviation of 0.15 (77 per cent of the population) together with a flat distribution in the range q = [0, 1] (23 per cent). Because, typically, the number of SBs in modern spectroscopic surveys is of the order of 100 – for example, Goldberg et al. (2003) analysed 129 SBs – we chose the size of the simulated sample to be 100 SB1 systems in both examples.

Derivation of MRDs from a simulation of a sample of 100 SB1s. The top panel shows the simulated sample of 100 SB1s, with random orientations, using as an MRD (dashed black line) a fourth-degree polynomial, fq(q) ∝ 25(2q − 1)4 + 4, that peaks at q = 0 and q = 1. The specific drawn sample is presented by a seven-bin histogram. The bottom panel shows two independently derived MRDs, one using the first five shifted Legendre polynomials as a basis (dashed black line) and the other using the boxcar basis of seven bins (thick black line). Each derived MRD is associated with an error for each value of q (see text).
Figure 5.

Derivation of MRDs from a simulation of a sample of 100 SB1s. The top panel shows the simulated sample of 100 SB1s, with random orientations, using as an MRD (dashed black line) a fourth-degree polynomial, fq(q) ∝ 25(2q − 1)4 + 4, that peaks at q = 0 and q = 1. The specific drawn sample is presented by a seven-bin histogram. The bottom panel shows two independently derived MRDs, one using the first five shifted Legendre polynomials as a basis (dashed black line) and the other using the boxcar basis of seven bins (thick black line). Each derived MRD is associated with an error for each value of q (see text).

Open in new tabDownload slide
Derivation of MRDs from a simulation of a sample of 100 SB1s. The simulated MRD is composed of a Gaussian with a mean at q = 0.2 and width of 0.15 (77 per cent of the population) and a flat part in the range q = [0, 1] (23 per cent). The simulated sample and the two derived MRDs are plotted as in Fig. 5.
Figure 6.

Derivation of MRDs from a simulation of a sample of 100 SB1s. The simulated MRD is composed of a Gaussian with a mean at q = 0.2 and width of 0.15 (77 per cent of the population) and a flat part in the range q = [0, 1] (23 per cent). The simulated sample and the two derived MRDs are plotted as in Fig. 5.

Open in new tabDownload slide

The best MRD model and its uncertainty were derived by calculating the median and scatter of the values obtained for each q along the MCMC run, as described in Section 3.1. The top panels of Figs 5 and 6 show the MRD used and the mass-ratio histogram of the simulated sample, while the bottom panels present the MRDs derived with a basis of seven boxcar functions and with the first five shifted Legendre polynomials.

An alternative method of deriving an explicit expression for the best-fitting model is to take the median value of each parameter obtained along the chain. For example, in terms of the shifted Legendre polynomials given in equation (18), for the two experiments presented above in this section the fitted MRDs are
\begin{eqnarray*} \skew4\hat{f}_1(q) &=& 1.22\, P_4(q) - 0.35\, P_3(q) + 1.54\, P_2(q)\\ && +\, 0.11\, P_1(q) + P_0(q) \end{eqnarray*}
and
\begin{eqnarray*} \hat{f}_2(q) &=& - 0.65\, P_4(q) + 0.99\, P_3(q) + 0.14\, P_2(q)\\ && -\, 1.22\, P_1(q) + P_0(q). \end{eqnarray*}
Here, |$\hat{f}_1$| and |$\hat{f}_2$| are the fitted models for the simulations presented in Figs 5 and 6, respectively. By gathering terms of the same power in q, the derived MRDs become
\begin{equation*} \hat{f}_1(q) = 85.6\, q^4 - 178.2\, q^3 + 129.8\, q^2 - 37.7\, q + 4.0 \end{equation*}
and
\begin{equation*} \hat{f}_2(q) = -45.7\, q^4 + 111.0\, q^3 - 87.4\, q^2 + 21.6\, q + 0.7. \end{equation*}
The differences between the MRDs derived by the two methods are found to be ≲σ/5.

The two examples demonstrate the power of our algorithm, as the derived MRDs are very close to the underlying functions, even though the algorithm was applied without any assumption on the shape of the MRD.

5 ACCOUNTING FOR AN OBSERVATIONAL DETECTION THRESHOLD

Samples of observed spectroscopic binaries are subjected to many observational biases. An obvious bias (e.g. Mazeh & Goldberg 1992; Tokovinin 1992) is the detection threshold – RV surveys can identify SB systems only if their RV amplitude is large enough. The impact of such a selection effect becomes increasingly significant for small values of q, causing the derived |$\hat{f_q}$| at small q values to be underestimated. In this section, we describe our way to account for this observational bias, following the approach of Mazeh & Goldberg (1992).

To model this effect, we assume that only (and all) binaries with RV amplitudes larger than some Kmin are detected. Therefore, for each q and m1, there exists the longest orbital period that can be detected:
(21)
Here, the orbits are assumed to be circular. For periods shorter than Pmax, the detectability depends on the inclination angle and therefore only a fraction of the population of binaries are detectable. We define the detection function D, which is the fraction of detected binaries out of all systems with identical P, m1 and q. The detection function is the probability for a system to have an inclination such that its observed RV semi-amplitude will be larger than the detection threshold,
(22)
Here, P is in d, m1 is in solar mass and Kmin is in km s−1.
Let us further assume that the primary stars in the sample are of nearly identical mass, |$\overline{m}$|⁠, and that the distribution of the orbital period, fP, is independent of q. The fraction of detected systems with some specific q is composed of the probability that both the period and the inclination allow a detection, namely
(23)
where P1 and P2 are the shortest and longest periods of the population, respectively.
The derived |$\hat{f_q}$| can now be corrected by the detection function,
(24)
where |$\hat{h}_{q}(q)$| is the unbiased distribution of q. This time, the corrected function has to be normalized in order to be used as a PDF. In the case of a boxcar fit, the correction factor of each bin is taken according to its value at the bin's centre. Notably, for very small mass ratios, the correction factor, |$1/\overline{D}(q)$|⁠, becomes very large and consequently uncertain. It is therefore advisable to cautiously address only a domain where the correction factor is a small number (i.e., |$1/\overline{D}(q) \lesssim 2$|⁠).

Again, to test the correcting part of the algorithm, we ran numerous simulations, one of which is presented in Fig. 7. Here, we used the same population as in Fig. 5, but now with a 1-M⊙ primary for each binary, orbital periods with log-uniform distribution between 1 and 103 d, and a detection threshold of Kmin = 3 km s−1, which caused 22 simulated binaries not to be detected. A histogram of the detected and missed binaries is plotted in the top panel of Fig. 7. As can be seen in the figure, most of the missed binaries had low mass ratio, as expected. The lower panel shows the uncorrected and corrected distributions. The uncorrected MRD suffers from serious suppression of its lower part, while the correction succeeded in producing the correct underlying MRD. As expected, for small mass ratios, q ≲ 0.05, the correction factor |$1/\overline{D}(q)$| became large, and therefore we refrained from obtaining the corrected function for this range of q.

MRD derivation from a simulated SB1 sample with a detection threshold of Kmin = 3 km s−1. The top panel shows a simulated sample of 100 SB1s with random orientations, using an MRD (dashed black line) that peaks at q = 0 and q = 1. The specific drawn sample is presented by a seven (grey) bin histogram (see Fig. 5). Because of the detection threshold, the upper bins of 22 systems (blue) were not detected. The bottom panel shows the derived MRDs (dotted line), using the first five shifted Legendre polynomials as a basis, together with the corrected MRD (dashed line) and its confidence 1σ range (see text).
Figure 7.

MRD derivation from a simulated SB1 sample with a detection threshold of Kmin = 3 km s−1. The top panel shows a simulated sample of 100 SB1s with random orientations, using an MRD (dashed black line) that peaks at q = 0 and q = 1. The specific drawn sample is presented by a seven (grey) bin histogram (see Fig. 5). Because of the detection threshold, the upper bins of 22 systems (blue) were not detected. The bottom panel shows the derived MRDs (dotted line), using the first five shifted Legendre polynomials as a basis, together with the corrected MRD (dashed line) and its confidence 1σ range (see text).

Open in new tabDownload slide

Another way to correct for the undetected binaries, not presented here, is to apply the derived |$\overline{D}(q)$| factor to the base functions, and to use these corrected functions along the MCMC fitting. One then constructs the true MRD by using the uncorrected base functions with the parameters obtained with the MCMC.

6 CONCLUSIONS

We have presented a novel direct algorithm to derive the MRD of short-period binaries from an observed SB1 sample. The algorithm considers a parametrized family of MRDs and finds the set of parameters that best fits the observed sample.

The algorithm consists of four parts. First, we define a new observable, the modified mass function S, derived for each SB1 in the sample. We show that the distribution of the modified mass function of an SB1 sample follows the shape of the underlying MRD, making the use of the modified mass function more advantageous than the previously used mass function, reduced mass function or the reduced mass function logarithm. Secondly, given an assumed MRD, we derive the likelihood of obtaining the observed sample of SB1s with the derived modified mass functions. Maximizing this likelihood by an MCMC search enables the algorithm to find the best parameters of the underlying MRD. Thirdly, we suggest that the unknown MRD is expressed by a basis of functions with some unknown coefficients that linearly span the space of possible MRDs. We suggest two such bases. Fourthly, we have shown how to account for the undetected systems that have an RV amplitude below a certain threshold. The correction is calculated per mass ratio and therefore can be applied to the derived MRD. Without the correction, this observational bias suppresses the derived MRD for low mass ratios. Numerous simulations show that the algorithm works with either of the two bases.

The algorithm is based on three simplifications. We consider only circular orbits, we ignore the double-lined binaries and we assume that there are no uncertainties associated with y and therefore with S. With the present layout, it is straightforward to generalize the algorithm to include eccentric orbits and uncertainties in S. However, it is an obvious drawback to ignore the extra information about the known mass ratio of the SB2s (e.g. Mazeh et al. 2003; Prato 2007; Fernandez et al. 2017). A further development of the algorithm to use the SB2 information is planned for a further publication. At present, the algorithm treats those systems as SB1s.

The detection threshold correction presented here depends on the orbital period distribution of the binary population and on the assumption that the MRD does not depend on the binary period; see the discussion by Moe & Di Stefano (2017), which put this assumption into question for O- and B-type stars. These two assumptions are inherent to any correction algorithm, as the RV amplitude does depend on the mass ratio and the orbital period. The simulations show that the correction succeeded in producing the correct MRD for low mass ratios.

The correction is based on a simplistic conception of the detection threshold. In reality, the observational bias does not act as a stiff threshold, but instead the detection probability of a binary is a continuous monotonic increasing function of its amplitude, which depends on the period, determined by the time stamp of the RV observations. However, it is easy to adopt the algorithm for any detection sensitivity through equations (22) and (23), by which one can derive a more sophisticated correction for any value of mass ratio. Needless to say, any derivation of the MRD can be based only on a sample that was obtained by a complete systematic survey that searches for spectroscopic binaries with known detection thresholds, so that the corrections can be derived and applied to the observed sample.

Obviously, the correction procedure introduces additional errors to the derived MRD, because of the inexact period distribution and inaccurate detection threshold used. Therefore, the correction becomes less valuable for low mass ratios, a range for which we have a small number of systems and for which the correction factor becomes large. In the simulated case presented above, for example, we refrained from plotting the corrected MRD for mass ratio smaller than 0.05. The exact limit depends on the specific SB1 sample.

In the next paper of this series (Shahaf et al., in preparation), we apply the algorithm to a few samples published in the literature, in particular those of Mazeh et al. (2003), Fisher et al. (2005), Prato (2007), Mermilliod et al. (2007) – see also North (2014) and Van der Swaelmen et al. (2017) – and Tal-Or, Faigler & Mazeh (2015).

Furthermore, in the near future we anticipate extremely large new samples of SBs coming from the APO Galactic Evolution Experiment (APOGEE)1 and the Gaia2 project. The release of the Gaia distances will enable us to better estimate the primary masses of these samples, which is a key element in the derivation of the reduced and modified mass functions. The new algorithm will be ready for these large samples to determine the MRDs of SBs. In addition, we anticipate two additional large samples, which will explore binaries with very short and very long period ranges: eclipsing binaries from large photometric data bases – see, for example Mazeh, Tamuz & North (2006) and Mowlavi et al. (2017) for the analysis of the OGLE LMC binaries – and astrometric binaries from Gaia. The new samples will finally give us the full picture of the different binary populations.

Acknowledgements

We are indebted to the referee for the thorough reading of the manuscript and very helpful comments. We thank Shay Zucker for the insightful discussions. We acknowledge support from the Israel Science Foundation (grant No. 1423/11) and the Israeli Centers of Research Excellence (I-CORE, grant No. 1829/12).

1

http://www.sdss3.org/surveys/apogee.php

2

http://sci.esa.int/gaia/

REFERENCES

Barlow
R.
,
1989
,
Least squares, in Statistics. A Guide to the Use of Statistical Methods in the Physical Sciences
.
Wiley
,
New York
, p.
105

Google Scholar

OpenURL Placeholder Text

 

Bate
M. R.
,
Bonnell
I. A.
,
1997
,
MNRAS
,
285
,
33

Crossref
Search ADS
 

Benacquista
M. J.
,
Downing
J. M. B.
,
2013
,
Living Reviews in Relativity
,
16
,
4

Crossref
Search ADS
PubMed
 

Boffin
H. M. J.
,
2010
,
A&A
,
524
,
A14

Crossref
Search ADS
 

Boffin
H. M. J.
,
2015
,
A&A
,
575
,
L13

Crossref
Search ADS
 

Curé
M.
,
Rial
D. F.
,
Cassetti
J.
,
Christen
A.
,
Boffin
H. M. J.
,
2015
,
A&A
,
573
,
A86

Crossref
Search ADS
 

Ducati
J. R.
,
Penteado
E. M.
,
Turcati
R.
,
2011
,
A&A
,
525
,
A26

Crossref
Search ADS
 

Duquennoy
A.
,
Mayor
M.
,
1991
,
A&A
,
248
,
485

 

Fernandez
M. A.
et al. ,
2017
,
PASP
,
129
,
084201

Crossref
Search ADS
 

Fisher
J.
,
Schröder
K.-P.
,
Smith
R. C.
,
2005
,
MNRAS
,
361
,
495

Crossref
Search ADS
 

Foreman-Mackey
D.
Hogg
D. W.
Lang
D.
Goodman
J.
,
2013
,
emcee: The MCMC Hammer
,
preprint (arXiv:1202.3665)

Google Scholar

OpenURL Placeholder Text

 

Goldberg
D.
,
Mazeh
T.
,
Latham
D. W.
,
2003
,
ApJ
,
591
,
397

Crossref
Search ADS
 

Goodman
J.
,
Weare
J.
,
2010
,
Communications in Applied Mathematics and Computational Science
,
5
,
65

Crossref
Search ADS
 

Grether
D.
,
Lineweaver
C. H.
,
2006
,
ApJ
,
640
,
1051

Crossref
Search ADS
 

Heacox
W. D.
,
1995
,
AJ
,
109
,
2670

Crossref
Search ADS
 

Hogeveen
S. J.
,
1992
,
Ap&SS
,
196
,
299

Crossref
Search ADS
 

Hut
P.
et al. ,
1992
,
PASP
,
104
,
981

Crossref
Search ADS
 

Kouwenhoven
M. B. N.
,
Brown
A. G. A.
,
Goodwin
S. P.
,
Portegies Zwart
S. F.
,
Kaper
L.
,
2009
,
A&A
,
493
,
979

Crossref
Search ADS
 

Lucy
L. B.
,
1974
,
AJ
,
79
,
745

Crossref
Search ADS
 

Lucy
L. B.
,
Ricco
E.
,
1979
,
AJ
,
84
,
401

Crossref
Search ADS
 

Mazeh
T.
,
Goldberg
D.
,
1992
,
ApJ
,
394
,
592

Crossref
Search ADS
 

Mazeh
T.
,
Simon
M.
,
Prato
L.
,
Markus
B.
,
Zucker
S.
,
2003
,
ApJ
,
599
,
1344

Crossref
Search ADS
 

Mazeh
T.
,
Tamuz
O.
,
North
P.
,
2006
,
MNRAS
,
367
,
1531

Crossref
Search ADS
 

Mermilliod
J.-C.
,
Andersen
J.
,
Latham
D. W.
,
Mayor
M.
,
2007
,
A&A
,
473
,
829

Crossref
Search ADS
 

Moe
M.
,
Di Stefano
R.
,
2017
,
ApJS
,
230
,
15

Crossref
Search ADS
 

Mowlavi
N.
et al. ,
2017
,
preprint (arXiv:1703.10597)

Nelsen
R.
,
2013
,
An Introduction to Copulas, Lecture Notes in Statistics
.
Springer
,
New York

Google Scholar

OpenURL Placeholder Text

 

North
P.
,
2014
, in
Mathys
G.
Griffin
E. R.
Kochukhov
O.
Monier
R.
Wahlgren
G. M.
, eds,
Putting A Stars into Context: Evolution, Environment, and Related Stars
.
Pero
,
Moscow
, p.
63
(arXiv:1309.7636)

Google Scholar

OpenURL Placeholder Text

 

Prato
L.
,
2007
,
ApJ
,
657
,
338

Crossref
Search ADS
 

Raghavan
D.
et al. ,
2010
,
ApJS
,
190
,
1

Crossref
Search ADS
 

Satsuka
T.
,
Tsuribe
T.
,
Tanaka
S.
,
Nagamine
K.
,
2017
,
MNRAS
,
465
,
986

Crossref
Search ADS
 

Scherrer
R. J.
,
Berlind
A. A.
,
Mao
Q.
,
McBride
C. K.
,
2010
,
ApJ
,
708
,
L9

Crossref
Search ADS
 

Tal-Or
L.
,
Faigler
S.
,
Mazeh
T.
,
2015
,
A&A
,
580
,
A21

Crossref
Search ADS
 

Terrell
G. R.
,
Scott
D. W.
,
1985
,
Journal of the American Statistical Association
,
80
,
209

Crossref
Search ADS
 

Tokovinin
A. A.
,
1992
,
A&A
,
256
,
121

 

Tokovinin
A. A.
,
2000
,
A&A
,
360
,
997

 

Toonen
S.
,
Nelemans
G.
,
Portegies Zwart
S.
,
2012
,
A&A
,
546
,
A70

Crossref
Search ADS
 

Trimble
V.
,
1990
,
MNRAS
,
242
,
79

Crossref
Search ADS
 

Van der Swaelmen
M.
,
Boffin
H. M. J.
,
Jorissen
A.
,
Van Eck
S.
,
2017
,
A&A
,
597
,
A68

Crossref
Search ADS
 

Yungelson
L. R.
,
Kuranov
A. G.
,
2017
,
MNRAS
,
464
,
1607

Crossref
Search ADS
 

APPENDIX A: DISTRIBUTION OF THE REDUCED MASS FUNCTION

The mass ratio of a binary system is q ≡ m2/m1, where m1 and m2 are the stellar masses of the primary and secondary, respectively. The reduced mass function y is
(A1)
where i is the inclination. Notably, each value of y is associated with a minimal possible q value, which can be determined by setting the inclination angle i to be 90°. This q minimum, denoted |$\mathcal {Q}_{y}$|⁠, is the only real root of the polynomial Py(q),
(A2)
The explicit expression for |$\mathcal {Q}_{y}$|⁠, as given by Heacox (1995), is
(A3)
where
(A4)

A rigorous development of the y PDF,  fy, for a sample of randomly oriented binaries has been previously presented by Heacox (1995). Nevertheless, an alternative geometrical derivation of this might be instructive in the context of this work.

We choose to work in the parameter plane of (1 − cos i, q), where 0 ≤ q ≤ 1 and 0 ≤ 1 − cos i ≤ 1, as the distribution of 1 − cos i is uniform for random orientation of the orbits. Equation (A1) implies that y values are uniquely associated with contours on that plane, as demonstrated in Fig. A1. A specific system with some given y ± δy/2 and q ± δq/2 occupies an area on the parameter plane,
(A5)
where by means of equation (A1), |∂cos i/∂y| is
(A6)
An example of δA assuming 0.0100 ≤ y ≤ 0.0101, at q = 0.3, is given in Fig. A1.
A contour of y = 0.01005 plotted in the (1 − cos i), q) plane. The red dot corresponds to the q minimum value of y = 0.01005. The grey circle indicates the point where q = 0.3. The enlarged window shows the area bounded by 0.0100 ≤ y ≤ 0.0101. The horizontal width of the parallelogram, δq, equals 0.0012 (dashed red). The vertical height of the parallelogram, δ(1 − cos i), equals $\mathbb {K}(y,q)\times \delta y = 0.0047$ (dashed red).
Figure A1.

A contour of y = 0.01005 plotted in the (1 − cos i), q) plane. The red dot corresponds to the q minimum value of y = 0.01005. The grey circle indicates the point where q = 0.3. The enlarged window shows the area bounded by 0.0100 ≤ y ≤ 0.0101. The horizontal width of the parallelogram, δq, equals 0.0012 (dashed red). The vertical height of the parallelogram, δ(1 − cos i), equals |$\mathbb {K}(y,q)\times \delta y = 0.0047$| (dashed red).

Open in new tabDownload slide
Because 1 − cos i is uniformly distributed, the probability of drawing a system with specific y and q values from a sample of randomly oriented binaries is ∼fq(q)δA, where fq is the sample's underlying MRD. Considering all possible values of q, taking δA to be infinitesimal and assuming 0 < q ≤ 1, fy becomes
(A7)
Some attempts have been made to use the PDF of log (y), flog y, as a more informative representation of the data (e.g. Boffin 2010, 2015). In terms of equation (A7), flog y is
(A8)

APPENDIX B: DERIVATION OF THE MODIFIED MASS FUNCTION

The modified mass function S is required to be a strictly increasing continuous transformation of y, from [0, 0.25] on to [0, 1]. Additionally, if the underlying MRD fq is uniform, then its resulting S distribution fS is also required to be uniform.

According to equation (2) cos i can be expressed in terms of y and q,
(B1)
The probability of observing a system at some y΄ < y is provided by integrating over the surface bounded by the axis, a contour of 1 − cos i(y, q) within the (1 − cos i, q) plane, namely
(B2)
Because 1 − cos i is uniformly distributed,
(B3)
The modified mass function is defined by taking equation (B3) with uniform MRD:
(B4)
The transformation |$\mathbb {S}$| is, by definition, the CDF of y assuming a uniform MRD, and therefore it obeys the requirements given at the beginning of this subsection.

|$\mathbb {S}$| is unique. Let |$\mathbb {T}$| and |$\mathbb {S}$| uphold the stated requirements. Because both are transformations of y, the PDFs obey fS|dS/dy| = fT|dT/dy|. Specifically for uniform fq, this relation becomes |dS/dy| = |dT/dy|. Because both are strictly increasing and continuous from [0, 0.25] on to [0, 1], |$\mathbb {T} \equiv \mathbb {S}$|⁠.

© 2017 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society
Issue Section:
Article
Download all slides