© ESO, 2015

1. Introduction

Extremely low-mass white dwarfs (ELMs) are thought to be the end products of binary star interactions and generally have a companion (e.g., Marsh et al., 1995; Rebassa-Mansergas et al., 2011). Although still far from understood, binary models predict that the companions to these ELMs are white dwarfs because the system has undergone one or two common-envelope phases (Nelemans et al., 2001; van der Sluys et al., 2006; Woods et al., 2012; Toonen et al., 2014). Knowing the distribution of the companion’s mass to these ELMs could provide useful constraints on the various parameters that enter binary star evolution and the common-envelope phase, as well as to predict whether and when they will merge and, more generally, how they will evolve (Dan et al., 2012; Kilic et al., 2012). Thanks to the ELM Survey, which has identified 61 ELMs and provided orbital parameters for 54 of them¹ (see Brown et al., 2013, and references therein), it is now possible to begin envisaging conducting a statistical analysis as was done by Andrews et al. (2014). The problem is that ELMs are single-lined spectroscopic binaries and as such, it is not possible to obtain the mass ratio directly (see, e.g., Boffin, 2010, 2012; Curé et al., 2014), but one needs to apply statistical methods to derive the mass-ratio (or companion mass) distribution. Andrews et al. (2014) developed a Bayesian probabilistic model to infer the companion mass distribution for the above-mentioned sample, assuming a functional form – a two-component Gaussian, with one component representing white dwarfs with masses between 0.2 and 1.44  M_⊙, and the other neutron stars with masses centred around 1.4  M_⊙ and a standard deviation of 0.05  M_⊙. Using a Markov chain Monte Carlo algorithm, they found that their best fit is given by a population of white dwarfs centred around 0.74  M_⊙ and a standard deviation of 0.24  M_⊙, without a neutron star. This is quite an interesting result that also indicates that in contrast to population synthesis models, the majority of companions to ELMs are CO-core WDs, and not another He WD. As such, it is important to examine whether this result holds when using different methods, including when the functional form is not fixed a priori. This is what I present here.

2. Methods

The derivation of the orbital elements for a single-lined spectroscopic binary (period, radial velocity amplitude, and eccentricity) allows obtaining the spectroscopic mass function, f(m), which is a combination of the masses of the two components and the (unknown) inclination of the orbit on the line of sight, i:

$\mathit{f}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{2}}^{\mathrm{3}}}{\mathrm{(}{\mathit{M}}_{\mathrm{1}}\mathrm{+}{\mathit{M}}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}{\sin }^{\mathrm{3}}\mathit{i,}$where M₁ is the mass of the primary (in this case, the ELM), and M₂, the mass of its companion. If M₁ is known, as it is the case here, then one can rewrite this as a function of the mass ratio, q = M₂/M₁:

$\mathit{Y}\mathrm{=}\frac{\mathit{f}\mathrm{\left(}\mathit{m}\mathrm{\right)}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\frac{{\mathit{q}}^{\mathrm{3}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{q}{\mathrm{)}}^{\mathrm{2}}}\mathrm{·}$The distribution of the logarithm of f(m) – or Y – can be used to determine the distribution of M₂ or q. This was done here using two different methods, in which I always assumed that the inclination i is randomly distributed on the sky, that is, P(i) = sini. I used the sample of Brown et al. (2013, similar to that used by Andrews et al. 2014, which provides a list of 54 systems with known f(m) and M₁.

In the first method, I assumed a functional form for the distribution of M₂: this is either a Gaussian with mean μ and standard deviation σ, or a uniform distribution, defined between a lower, M_2,l, and an upper, M_2,u value of the companion mass. I then applied a Kolmogorov-Smirnov statistical test. I assumed a single population of companions, without distinguishing between a neutron star (NS) and a white dwarf companion population, for instance. This is based on the fact that Andrews et al. (2014) found the NS fraction to be very low, as I confirm here as well.

In each case, I ran 10 000 Monte Carlo simulations, where M₁ is distributed according to the observed distribution (see Fig. A.1), M₂ according to the chosen distribution with one set of parameters, and i is assumed to be randomly distributed on the sky. For each sample of simulations, I calculated the cumulative distribution of log f(m), which I compared with the observed one. To do this, I calculated the largest deviation between the two distributions, D^∗ = D_{54,10 000}. Running 10 000 simulations provides very good precision, and there would be no gain in running more, as the estimator D^∗ saturates for large numbers².

This estimator allows determining the probability with which the simulated and the observed distribution are extracted from the same population. Thus, if D^∗> 0.2628, the chance is only 0.1% that the two populations are extracted from the same population, and we may most likely ignore such solutions. For D^∗> 0.1274 and D^∗> 0.0506, these probabilities become 33% and 99.9%, respectively. The first value provides a 1σ estimate of the parameters that are allowed, while the latter number can be used to estimate the 0.1% confidence interval of parameters that provide a very good match to the observed distribution.

In the second method, I took a more direct approach and inverted the distribution of Y to derive the mass-ratio distribution.

3. Results

3.1. Functional form

3.1.1. Gaussian

Fig. 1 Standard deviation, σ, versus the mean mass of the companion, μ, for all simulated samples where the companion mass is distributed according to a Gaussian, and which have D∗ < 0.1274 (in black) and those that have D∗ < 0.0506 (in green). The heavy, red dot shows the location μ = 0.74, σ = 0.24. (A colour version of the figure is available online.)

As mentioned above, I here assumed that the distribution of the companion mass, Φ(M₂), is given by a Gaussian:

$\mathrm{\Phi }\mathrm{\left(}{\mathit{M}}_{\mathrm{2}}\mathrm{\right)}\mathrm{\propto }\exp \left(\mathrm{-}\frac{{\mathit{M}}_{\mathrm{2}}\mathrm{-}\mathit{\mu }}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right)\mathit{,}$and I determined for a plane 0 <μ < 2,0 <σ < 1 the value of D^∗. The relation of D^∗ with μ is shown in Fig. A.2, where I also show the lines corresponding to D^∗ = 0.0506, 0.1274, and 0.2628. This figure shows that one can find simulations spanning the whole range 0 <μ < 1.22 M_⊙ that would lead to D^∗ < 0.2628, while the 1σ range covers the values 0.3 <μ < 0.92 M_⊙. This indicates that it is difficult to constrain the parameters of such a functional form based on the small observed sample. It is clear, however, that the simulations with μ between 0.7 and 0.8  M_⊙ correspond to the lowest values of D^∗, with the minimum being around μ = 0.76 M_⊙ and σ = 0.27 M_⊙, leading to D^∗ = 0.027. This value indicates that this simulation and the observed distribution are consistent at a very high level because the null hypothesis that the two populations are drawn from the same population can only be rejected at a level of 1.3 × 10^-12%.

Of course, the standard deviation σ is correlated with μ, and Fig. 1 shows the solutions that lead to D^∗ < 0.0506 and D^∗ < 0.1274, illustrating the allowed range. For D^∗ < 0.1274, I found that the mean value of the companion mass is μ_m = 0.70 ± 0.12 M_⊙, while the mean value of the standard deviation is σ_m = 0.45 ± 0.22 M_⊙. For D^∗ < 0.0506, these values become μ_m = 0.76 ± 0.03 M_⊙ and σ_m = 0.27 ± 0.06 M_⊙.

A χ² analysis that computes the deviation of the computed and observed distribution of log f(m) provides a similar result, with the best fit being given by μ = 0.76 ± 0.02 M_⊙ and σ = 0.27 ± 0.02 M_⊙ (see Fig. A.4).

Fig. 2 Highest mass of the companion as a function of the lowest mass of the companions for all simulated uniformly distributed samples that have D∗ < 0.1274 (in black) and those that have D∗ < 0.0506 (in green). The two masses are clearly correlated, and I show a linear fit to the green dots as the heavy, solid red line. (A colour version of the figure is available online.)

3.1.2. Uniform distribution

I repeated the same analysis but using a uniform distribution of the companion mass, between a lowest (M_2,l) and a highest value (M_2,u), which were assumed to be in the range 0 <M_2,l < 0.9 M_⊙ and M_2,l<M_2,u<M_2,l + 1.5 M_⊙. The results are shown in Figs. A.3 and 2. Again, the range of allowed values is very wide. If this is restricted to D^∗ < 0.2628, the whole range of M_2,l is allowed, while for M_2,u, it is restricted to values between 0.5 and 1.95  M_⊙. The range becomes narrower for lower values of D^∗, as shown in Fig. 2, which also shows that the acceptable values of M_2,l and M_2,u are correlated. A linear fit gives

${\mathit{M}}_{\mathrm{2}\mathit{,}\mathrm{u}}\mathrm{=}\mathrm{-}\mathrm{1.25}{\mathit{M}}_{\mathrm{2}\mathit{,}\mathrm{l}}\mathrm{+}\mathrm{1.60.}$The mean values of M_2,l and M_2,u are 0.36 ± 0.18 M_⊙ and 1.20 ± 0.23 M_⊙ for D^∗ < 0.1274, and 0.33 ± 0.07 M_⊙ and 1.19 ± 0.10 M_⊙ for D^∗ < 0.0506. The lowest value of D^∗ (0.02793) is reached for M_2,l = 0.25 M_⊙ and M_2,u = 1.28 M_⊙. This value corresponds to a probability that the two distributions are not extracted from the same population of 4 × 10^-10%. Perhaps most importantly, this shows that a uniform distribution of the companion mass fits the data very well and cannot be discarded. With such a functional form, a Gaussian or a uniform distribution cannot be preferred; they do, of course, quite overlap (see Fig. A.5).

Fig. 3 Comparison between the observed distribution of the logarithm of the spectroscopic mass function (solid black line) and the best fits for the uniformly distributed (red dotted line connected by heavy dots) and Gaussian-distributed (green dashed line connected by open squares) companion masses. For the former, I used M2,l = 0.25 M⊙ and M2,u = 1.28 M⊙, while for the latter, I used μ = 0.76 M⊙ and σ = 0.27 M⊙. The top panel shows the fraction of systems, while the bottom panel is the cumulative fraction of systems. It is clear that both samples are good fits, given the intrinsic errors of the observed distribution. (A colour version of the figure is available online.)

This is further illustrated in Fig. 3, which shows the distribution of log f(m) for the observed sample as well as for the two best fits of the functional form, a Gaussian, and a uniform distribution. It is clear that both simulated distributions are good fits to the observed distribution, while it is very hard to distinguish the results from the two different functional forms.

3.1.3. Inversion method

Using a functional form for the mass distribution allows a better control of the systematics of a method, but at the cost of risking missing some interesting deviations. This is for example the case in the paper by Andrews et al. (2014) in their test 4, where they compare the distribution they obtain for a sample of post-common envelope systems to that obtained from spectroscopy. Although they clearly reproduce the bulk distribution, they miss the tail and other details of the distribution (see their Fig. 3). Moreover, Andrews et al. (2014) also pointed out that for the sample of ELMs, their result “could indicate that the true WD distribution may not be exactly Gaussian” – and indeed the previous section has shown that a uniform distribution also provides a good fit to the observations. It is therefore useful to consider exploring methods that do not require the a priori input of a functional form. This is the case of the Richardson-Lucy (R-L) inversion method, as used by Boffin et al. (1993); Boffin (2010, 2012), and Curé et al. (2014).

Fig. 4 Distribution of the mass ratios as determined by the Richardson-Lucy algorithm (solid black line) and as obtained for the best fits with a Gaussian (red dotted line) and a uniform (blue dashed line) distribution of companions. (A colour version of the figure is available online.)

I refer to these papers and references therein for a full discussion of the method, and in particular to Cerf & Boffin (1994) for a more formal presentation. Here, I just mention that the Richardson-Lucy method relies on the Bayes theorem on conditional probabilities and solves the Fredholm integral equation that links Y with q by an iterative scheme. It is important to note that we do not directly have the distribution of the companion mass, but only of the mass ratio. This is, however, a very important parameter for binary evolution models, and given that the mass of the primary is very peaked at 0.17  M_⊙, the distribution of the mass ratio can give an idea of the distribution of the companion mass. The outcome of this method is shown in Fig. 4, where I also compare it with the distribution derived from the best fits of the functional form method. For the latter, I derived the mass-ratio distribution by using the functional form for the companion mass and the primary mass distribution as determined by Brown et al. (2013).

All three methods provide a rather similar mass-ratio distribution, with some small differences. The outcome of the R-L method gives a broad Gaussian distribution, but with two more pronounced peaks, around q ~ 2−2.5 and around q ~ 4−5. This may indicate that a simple functional form may be missing on some small structure in the data, although a much larger sample is needed to be able to confirm these (Boffin, 2012). A K-S test indicates that the hypothesis that the functional forms and the outcome of the R-L method are drawn from the same population can be rejected at the level of 92% – a rather high number, but perhaps not convincing enough. Indeed, given the small sample, this is at most a two-sigma result. However, it shows that the true mass-ratio distribution (and companion mass distribution) may have a more complicated structure than any simple functional form we can think of. Only with much larger samples will be able to know this.

At the suggestion of the referee, I have also examined the mass-ratio distribution derived with the R-L method when limiting to the least massive primary stars, that is, the 33 systems with masses lower than 0.2  M_⊙. The resulting distribution is shown in Fig. A.6. It shows a single-peaked distribution centred around q ~ 4−5, that is (given that we now examine systems with a primary mass M₁ = 0.17 M_⊙), M₂ = 0.68−0.85 M_⊙. In the figure, the corresponding mass-ratio distributions for the functional forms were computed with a single value of the primary mass. The outcome of the R-L method apparently agrees better with the Gaussian distribution, but given that we have now an even smaller sample, the data should not be overinterpreted because the three distributions are again compatible within 2σ. This peak may correspond to the similar peak seen in the mass-ratio distribution seen for the entire sample, while the fact that there is a possible second peak in Fig. 4 at a lower mass ratio is most probably linked to the more massive primaries.

The very small excess of high mass ratios seen in Figs. 4 and A.6 should not be given too much importance: firstly, it is a well-known effect of the R-L method to smooth the distribution, the effect becoming weaker depending on the number of iterations (Cerf & Boffin, 1994), and secondly, their value is compatible with zero, given the size of the sample.

4. Discussion and conclusions

The results of Andrews et al. (2014) in the study of ELM WDs are very important, therefore I have reanalysed the same sample of spectroscopic binaries they used with two different techniques. In the first one, I assumed functional forms and applied a K-S statistics. The obtained results were confirmed by a χ² statistical test. The parameters for the Gaussian functional form are very similar to those found using a different method by Andrews et al. (2014), but I showed that a uniform distribution of the companion mass can provide as good a fit to the observed data and that the range of

parameters allowed is rather large, making it hard to provide a definitive answer as to the real distribution. If the uniform distribution illustrated in Fig. A.5 is more representative of the companion mass distribution, then more double He WD binary systems may be expected, such as the eclipsing double white dwarf binary CSS 41177 (Bours et al., 2014): the uniform distribution shown in this figure leads to a 24% probability (this is the fraction of systems that have a secondary mass lower than 0.5  M_⊙), compared to 16% for the Gaussian fit shown in the same figure or as derived by Andrews et al. (2014). If I take into account all the possible values at the 1σ level (i.e. those with D^∗ < 0.1274), I also derive a 26% probability for both my Gaussian and uniform distributions.

In addition, I applied an inversion method to derive the mass-ratio distribution, without the need to assume any functional form. The results are compatible with those derived from the functional form, although they seem to indicate some additional fine-grained structures.

An important question that Andrews et al. (2014) addressed is the possibility to have a neutron star as a companion to the ELM WD. Andrews et al. (2014) found that the probability of this is very low. This is also confirmed by my results. For example, the uniform distribution of companion masses indicates that the range allowed for D^∗ < 0.0506 ends at 1.4  M_⊙. If this is relaxed to D^∗ < 0.1274, however, the highest companion mass can be up to 1.7  M_⊙, but all in all, the probability to have one system with such a high mass (>1.44  M_⊙) is very low: for the uniform distribution, the 1σ probability is zero, while for the Gaussian distribution it is 6.5%.

Online material

Appendix A: Additional figures

Fig. A.1 Distribution of primary mass in the sample of Brown et al. (2013). The sample contains many objects with primary mass equal to 0.17  M⊙.

Fig. A.2 Maximum deviation between the observed sample with N = 54 and a synthetic sample containing N′ = 10 000 objects, with the companion mass being distributed as a Gaussian with a mean mass as given by the abscissa. The blue dotted line is drawn at D∗ = 0.2628: all simulations with a higher value of the estimator can be rejected at the 99.9% level as being drawn from the same population as the observed distribution. The magenta dashed line is drawn at D∗ = 0.1274 and provide a 1σ estimate of the allowed parameter range, while the red solid line is drawn at D∗ = 0.0506.

Fig. A.3 Top: maximum deviation as a function of the highest mass of the companion if the latter is uniformly distributed between a lowest and highest mass. Bottom: the same, but as a function of the lowest mass. The lines are drawn at the same value of D∗ as in Fig. A.2.

Fig. A.4 Reduced χ2 map in the mean companion mass vs. standard deviation for an assumed Gaussian functional form of the companion mass distribution – contours are shown at levels of 1, 1.5, 2, and 2.5. The red dot shows the position of the best fit of Andrews et al. (2014). The χ2 map confirms the result that was obtained with the K-S statistics.

Fig. A.5 Comparison between the two best functional forms that lead to the lowest values of D∗: a Gaussian distribution with μ = 0.76 M⊙ and σ = 0.27 M⊙ (solid black line) and a uniform distribution with M2,l = 0.25 M⊙ and M2,u = 1.28 M⊙ (dotted red line).

Fig. A.6 Distribution of the mass ratios as determined by the Richardson-Lucy algorithm (solid black line) for all systems with a primary mass below 0.2  M⊙, and as obtained for the best fits with a Gaussian (red dotted line) and a uniform (blue dashed line) distribution of companions.

Acknowledgments

I would like to thank the referee, J.J. Andrews, for a careful reading of the manuscript and for providing suggestions to improve the paper.

References

All Figures

	Fig. 1 Standard deviation, σ, versus the mean mass of the companion, μ, for all simulated samples where the companion mass is distributed according to a Gaussian, and which have D∗ < 0.1274 (in black) and those that have D∗ < 0.0506 (in green). The heavy, red dot shows the location μ = 0.74, σ = 0.24. (A colour version of the figure is available online.)
In the text

	Fig. 2 Highest mass of the companion as a function of the lowest mass of the companions for all simulated uniformly distributed samples that have D∗ < 0.1274 (in black) and those that have D∗ < 0.0506 (in green). The two masses are clearly correlated, and I show a linear fit to the green dots as the heavy, solid red line. (A colour version of the figure is available online.)
In the text

Fig. 3 Comparison between the observed distribution of the logarithm of the spectroscopic mass function (solid black line) and the best fits for the uniformly distributed (red dotted line connected by heavy dots) and Gaussian-distributed (green dashed line connected by open squares) companion masses. For the former, I used M2,l = 0.25 M⊙ and M2,u = 1.28 M⊙, while for the latter, I used μ = 0.76 M⊙ and σ = 0.27 M⊙. The top panel shows the fraction of systems, while the bottom panel is the cumulative fraction of systems. It is clear that both samples are good fits, given the intrinsic errors of the observed distribution. (A colour version of the figure is available online.)

In the text

	Fig. 4 Distribution of the mass ratios as determined by the Richardson-Lucy algorithm (solid black line) and as obtained for the best fits with a Gaussian (red dotted line) and a uniform (blue dashed line) distribution of companions. (A colour version of the figure is available online.)
In the text

	Fig. A.1 Distribution of primary mass in the sample of Brown et al. (2013). The sample contains many objects with primary mass equal to 0.17  M⊙.
In the text

Fig. A.2 Maximum deviation between the observed sample with N = 54 and a synthetic sample containing N′ = 10 000 objects, with the companion mass being distributed as a Gaussian with a mean mass as given by the abscissa. The blue dotted line is drawn at D∗ = 0.2628: all simulations with a higher value of the estimator can be rejected at the 99.9% level as being drawn from the same population as the observed distribution. The magenta dashed line is drawn at D∗ = 0.1274 and provide a 1σ estimate of the allowed parameter range, while the red solid line is drawn at D∗ = 0.0506.

In the text

	Fig. A.3 Top: maximum deviation as a function of the highest mass of the companion if the latter is uniformly distributed between a lowest and highest mass. Bottom: the same, but as a function of the lowest mass. The lines are drawn at the same value of D∗ as in Fig. A.2.
In the text

	Fig. A.4 Reduced χ2 map in the mean companion mass vs. standard deviation for an assumed Gaussian functional form of the companion mass distribution – contours are shown at levels of 1, 1.5, 2, and 2.5. The red dot shows the position of the best fit of Andrews et al. (2014). The χ2 map confirms the result that was obtained with the K-S statistics.
In the text

	Fig. A.5 Comparison between the two best functional forms that lead to the lowest values of D∗: a Gaussian distribution with μ = 0.76 M⊙ and σ = 0.27 M⊙ (solid black line) and a uniform distribution with M2,l = 0.25 M⊙ and M2,u = 1.28 M⊙ (dotted red line).
In the text

	Fig. A.6 Distribution of the mass ratios as determined by the Richardson-Lucy algorithm (solid black line) for all systems with a primary mass below 0.2  M⊙, and as obtained for the best fits with a Gaussian (red dotted line) and a uniform (blue dashed line) distribution of companions.
In the text