The Time-Domain Spectroscopic Survey: Radial Velocity Variability in Dwarf Carbon Stars

TL;DR

This study investigates the binary nature of dwarf carbon stars using radial velocity data, revealing a high binary fraction and supporting their formation through mass transfer in close binary systems.

Contribution

It provides the first comprehensive analysis of binary frequency and orbital properties of dwarf carbon stars using time-domain spectroscopy.

Findings

Dwarf carbon stars have a binary fraction of approximately 95%.

Most dC binaries have separations less than 1 AU, with periods around 1 year.

Radial velocity variability confirms mass transfer as a key formation mechanism.

Abstract

Dwarf carbon (dC) stars, main sequence stars showing carbon molecular bands, were initially thought to be an oxymoron since only AGB stars dredge carbon into their atmospheres. Mass transfer from a former AGB companion that has since faded to a white dwarf seems the most likely explanation. Indeed, a few types of giants known to show anomalous abundances --- notably, the CH, Ba and CEMP-s stars --- are known to have a high binary frequency. The dC stars may be the enhanced-abundance progenitors of most, if not all, of these systems, but this requires demonstrating a high binary frequency for dCs. Here, for a sample of 240 dC stars targeted for repeat spectroscopy by the SDSS-IV's Time Domain Spectroscopic Survey, we analyze radial velocity variability to constrain the binary frequency and orbital properties. A handful of dC systems show large velocity variability ($>$100 km s$^{-1}$). We compare the dCs to a control sample with a similar distribution of magnitude, color, proper motion, and parallax. Using MCMC methods, we use the measured $\Delta$RV distribution to estimate the binary fraction and the separation distribution assuming both a unimodal and bimodal distribution. We find the dC stars have an enhanced binary fraction of 95\%, consistent with them being products of mass transfer. These models result in mean separations of less than 1 AU corresponding to periods on the order of 1 year. Our results support the conclusion that dC stars form from close binary systems via mass transfer.

Tables

Table 1. Table 1: dC Sample Properties

Index	${\alpha }_{J2000}$		${\delta }_{J2000}$		$r$		$g-r$		Plate1	MJD1	FiberID1	Plate2	MJD2	FiberID2	$\Delta$RV	$\Delta$RVerror
	[Deg]		[Deg]												[km s-1]	[km s-1]
1	0\@alignment@align.0626	28.1693	17\@alignment@align.97	0.74	2824\@alignment@align	54452	253\@alignment@align	7696	57655	83	5	13
2	0\@alignment@align.1483	-0.1875	18\@alignment@align.74	1.42	1489\@alignment@align	52991	156\@alignment@align	7850	56956	704	-18	22
3	0\@alignment@align.9212	23.9270	16\@alignment@align.94	1.74	2801\@alignment@align	54331	201\@alignment@align	7665	57328	734	25	14
4	1\@alignment@align.2380	1.1606	18\@alignment@align.33	1.14	1490\@alignment@align	52994	364\@alignment@align	7862	56984	626	6	11
5	3\@alignment@align.1909	-1.0895	19\@alignment@align.47	1.59	687\@alignment@align	52518	297\@alignment@align	7863	56975	804	44	33
6	6\@alignment@align.8651	6.6597	17\@alignment@align.47	0.41	3106\@alignment@align	54714	139\@alignment@align	3106	54738	131	14	25
7	7\@alignment@align.3310	0.7206	18\@alignment@align.84	0.7	1087\@alignment@align	52930	504\@alignment@align	7855	57011	424	26	31
8	9\@alignment@align.3389	0.2069	19\@alignment@align.01	0.88	1495\@alignment@align	52944	353\@alignment@align	7868	57006	812	-252	15
9	9\@alignment@align.9056	15.4863	18\@alignment@align.48	0.95	419\@alignment@align	51812	346\@alignment@align	419	51868	346	11	26
10	10\@alignment@align.8422	0.4788	16\@alignment@align.72	0.55	1904\@alignment@align	53682	495\@alignment@align	7870	57016	562	13	8

Table 2. Table 2: Control Sample Properties

Index	${\alpha }_{J2000}$		${\delta }_{J2000}$		$r$		$g-r$		Plate1	MJD1	FiberID1	Plate2	MJD2	FiberID2	$\Delta$RV	$\Delta$RVerror
	[Deg]		[Deg]												[km s-1]	[km s-1]
1	47\@alignment@align.6415	-0.2745	17\@alignment@align.91	0.73	2068\@alignment@align	53386	74\@alignment@align	7255	56597	160	1	9
2	172\@alignment@align.5454	20.1782	18\@alignment@align.7	1.43	3170\@alignment@align	54859	640\@alignment@align	3170	54907	582	-54	26
3	0\@alignment@align.9212	23.9269	16\@alignment@align.94	1.74	2801\@alignment@align	54331	201\@alignment@align	7665	57328	734	-18	19
4	35\@alignment@align.3647	-0.2364	18\@alignment@align.22	1.18	704\@alignment@align	52205	234\@alignment@align	703	52209	30	-53	19
5	170\@alignment@align.6615	45.6529	19\@alignment@align.65	1.57	3216\@alignment@align	54853	143\@alignment@align	3216	54908	156	66	55
6	6\@alignment@align.8651	6.6597	17\@alignment@align.47	0.41	3106\@alignment@align	54714	139\@alignment@align	3106	54738	131	22	11
7	10\@alignment@align.6953	-0.6110	18\@alignment@align.83	0.7	1905\@alignment@align	53613	213\@alignment@align	1905	53706	219	32	27
8	113\@alignment@align.8683	41.4378	19\@alignment@align.07	0.92	3658\@alignment@align	55205	890\@alignment@align	5941	56193	884	12	38
9	44\@alignment@align.3177	0.9009	18\@alignment@align.39	0.95	1512\@alignment@align	53035	590\@alignment@align	1512	53742	579	-2	26
10	108\@alignment@align.2489	38.7804	16\@alignment@align.77	0.54	2938\@alignment@align	54503	6\@alignment@align	2938	54526	18	3	8

Table 3. Table 3: XD-GMM Component Fits

Parameter	dC		Control
${\alpha }_1$	0\@alignment@align.688	1.00
${\mu }_1$ [km s-1]	2\@alignment@align.02	9.69
${\sigma }_1$ [km s-1]	251\@alignment@align.53	1035.82
${\alpha }_2$	0\@alignment@align.312	.
${\mu }_2$ [km s-1]	2\@alignment@align.36	.
${\sigma }_2$ [km s-1]	12365\@alignment@align.92	.

Equations

$$f(a)=\frac{1}{a\sigma \sqrt{2\pi }}\exp \left(-\frac{(\ln a-\mu {)}^2}{2{\sigma }^2}\right)$$

$$\begin{gathered} f(a)=\frac{\alpha }{a{\sigma }_1\sqrt{2\pi }}\exp \left(-\frac{(\ln a-{\mu }_1{)}^2}{2{\sigma }_1^2}\right) \\ +\frac{1-\alpha }{a{\sigma }_2\sqrt{2\pi }}\exp \left(-\frac{(\ln a-{\mu }_2{)}^2}{2{\sigma }_2^2}\right) \end{gathered}$$

Full text

The Time-Domain Spectroscopic Survey: Radial Velocity Variability in Dwarf Carbon Stars

Benjamin R. Roulston

SAO Predoctoral Fellow

Department of Astronomy, Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA

Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA

Paul J. Green

Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA

John J. Ruan

McGill University, 3550 University Street Montreal, QC H3A 2A7, Canada

Chelsea L. MacLeod

Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA

Scott F. Anderson

Department of Astronomy, University of Washington, Box 351580, Seattle, WA 98195, USA

Carles Badenes

Department of Physics and Astronomy and Pittsburgh Particle Physics, Astrophysics and Cosmology Center (PITT PACC), University of Pittsburgh, 3941 O’Hara St, Pittsburgh, PA 15260, USA

Joel R. Brownstein

Department of Physics and Astronomy, University of Utah, 115 S. 1400 E., Salt Lake City, UT 84112, USA

Donald P. Schneider

Department of Astronomy and Astrophysics, The Pennsylvania State University, 525 Davey Laboratory, University Park, PA 16802, USA

Keivan G. Stassun

Department of Physics & Astronomy, Vanderbilt University, 6301 Stevenson Center Lane, Nashville, TN 37235, USA

(Received December 19, 2018; Revised April 11, 2019; Accepted April 15, 2019)

stars: binaries — stars: carbon — stars: chemically peculiar

††journal: The Astrophysical Journal††software: Astropy (The Astropy Collaboration et al., 2018), corner (Foreman-Mackey, 2017), IRAF (Tody, 1993), matplotlib (Hunter, 2007), Numpy (Oliphant, 2006), Scipy (Jones et al., 2001–), Scikit-learn (Pedregosa et al., 2011)

\turnoffediting

1 Introduction

The first carbon star was discovered by Secchi (1869) showing strong bandheads of C2 in their optical spectra. The phrase “dwarf carbon star” would have long been considered an oxymoron since it was assumed that all carbon stars were giants AGB stars that had dredged up carbon produced in their cores which polluted their atmospheres, producing an atmospheric C/O ratio above unity. Carbon preferentially binds with oxygen to form CO. Remaining carbon can combine to form the carbon compounds C2, CN, and CH giving rise to strong carbon molecular bandheads dominating their optical spectra.

This assumption was shown to be invalid when Dahn et al. (1977) discovered the first main-sequence carbon star, the dwarf carbon (dC) star G77-61. How could a main-sequence star have enough carbon in its atmosphere to push the C/O ratio beyond unity and create the strong carbon bandheads seen? The currently favored hypothesis is that dC stars are the products of binary mass transfer where the former AGB companion has become a white dwarf, leaving a carbon-enhanced dC primary. Indeed, about a dozen “smoking gun” systems, having composite spectra with a hot DA white dwarf component (Heber et al., 1993; Liebert et al., 1994; Green, 2013; Si et al., 2014), bolster this hypothesis.

As main sequence stars with carbon-enriched atmospheres, dC stars are probably the progenitors of the typically more luminous carbon-enhanced metal-poor (CEMP), sgCH, CH, and, perhaps, barium (Ba II) stars, which all have carbon and $s$-process enhancements (see discussion and references in Jorissen et al. 2016; De Marco & Izzard 2017). Samples of such stars have all been targets of spectroscopic radial velocity (RV) monitoring campaigns to test for binarity and characterize orbits. Sperauskas et al. (2016) recently studied the CH-like stars — giants showing CH and $s$-process abundances like CH stars, but without their halo kinematics — and showed $\sim 6$ times higher RV variability for C-rich than C-normal giants. Leiner et al. (2017) and others reported that binary interactions occur with surprising frequency, modifying the evolution of a considerable fraction of stars. Binary systems with mass transfer follow an impressive variety of evolutionary channels and may result in important and spectacular systems such as luminous red novae, Type Ia supernovae, planetary nebulae, and more. Here, we demonstrate that dCs definitively belong to the family of mass transfer binary systems.

The first dC star confirmed to be in a binary system, with a measured period of 245d, is the dC prototype G77-61 (Dearborn et al., 1986). Margon et al. (2018) recently mapped the orbit of a second system with a much shorter (3d) period. Whitehouse et al. (2018) detected radial velocity variability in 21 of 28 dC stars monitored spectroscopically with the William Herschel Telescope between 2013 and 2017.

We seek here, using multi-epoch spectroscopy of many dC stars discovered by the Sloan Digital Sky Survey (SDSS; Blanton et al., 2017), to measure the dC binary frequency, which should be near unity in this mass transfer scenario. As was the case in Whitehouse et al. (2018), our SDSS sampling lacks enough epochs to determine individual orbital parameters. However, with a significantly larger sample of 240 dC stars, we can use the distribution of radial velocity variations and Markov chain Monte Carlo methods to characterize the dC population’s separation distribution as was done by Maoz et al. (2012) and Badenes & Maoz (2012).

2 Dwarf Carbon Star Sample Selection

Dwarf carbon stars for this study were selected from the Green (2013) and Si et al. (2014) carbon star samples. Green (2013) identified carbon stars by visual inspection of single-epoch SDSS spectra compiled from the union of (1) SDSS DR7 spectra (Abazajian et al., 2009) having strong cross-correlation coefficients with the SDSS carbon star templates with (2) SDSS spectra with a DR8 pipeline class of STAR and subclass including the word carbon (Aihara et al., 2011). From within this carbon star parent sample, definitive main sequence dwarfs were selected by having significant proper motions ($\geq 3\sigma$ and 11 mas yr*-1*) in the catalog of Munn et al. (2004) and/or having SDSS spectra visibly identifiable as composite DA/dC spectroscopic binaries . Si et al. (2014) selected dCs using

For our current work, which aims to measure radial velocity variability, the primary additional selection criterion was that the selected dC stars have more than one epoch of spectroscopy in the SDSS as of November 2017. The majority of such objects were intentionally targeted for a second epoch of spectroscopy with the Time Domain Spectroscopic Survey (TDSS; Morganson et al., 2015a), a subprogram of the SDSS-IV extended Baryon Acoustic Oscillation Sky Survey (eBOSS; Dawson et al., 2016) project. Within TDSS, the main single-epoch-spectroscopy (SES) program (Morganson et al., 2015b) — along with its pilot survey, dubbed SEQUELS within SDSS-III (Ruan et al., 2016) — primarily targets optical point sources (unconfirmed quasars and stars) for a first epoch of spectroscopy based on variability. However, within several “few-epoch spectroscopy” (FES) subprograms, TDSS also acquires repeat spectroscopic observations for subsets of known stars and quasars that are astrophysically interesting. The FES programs are described by MacLeod et al. (2018) and include several classes of quasars and stars, including dC stars, re-targeted to study their spectroscopic variability. For the dC FES program, we selected all 730 SDSS dC stars from Green (2013) as well as another 99 dC stars found by Si et al. (2014), totaling 829 unique dC stars

We then searched the SDSS database (using the CasJobs query tool) for spectroscopy from DR14 (Abolfathi et al., 2018) for the dC stars that have been observed in the dC FES program. We also checked the DR14 database for all dC stars in the Green (2013) sample in search of any dC stars that may have been observed more than once, but not as part of the TDSS FES program. Our final sample contains FES spectra obtained up until October 31, 2017 and spectra from DR14.

We visually inspected all spectral epochs and removed any spectra that had strong broad artifacts. The final sample for this study was 240 dCs with a total of 540 spectra within the SDSS.

Figure 1 shows the correlation between $r$ mag (Fukugita et al., 1996) and the median spectroscopic signal-to-noise ratio (S/N). The color axis is the Modified Julian Date (MJD) of each epoch, and a clear distinction can be seen between the early epochs and later epochs in regards to S/N due to the enhanced capabilities of the BOSS spectrograph.

3 Control Sample Selection

3.1 Selection Criteria

The control sample was selected from the SDSS DR14 catalog using the properties of the dC sample as a selection criteria. The control sample criteria were as follows: (1) objects must have CLASS=STAR from the SDSS spectroscopic pipeline (2) a significant proper motion detection following the criteria of Green (2013)111Proper motion in at least one coordinate larger than $3\sigma$ where $\sigma$ is the proper motion uncertainty in that coordinate, and total proper motion larger than 11 mas yr*-1*. (3) select only stars within the 2 – 98% parameter ranges of the dC sample (i.e., total proper motion between 11 and 143 milliarcseconds yr*-1*, SDSS $r$ mag between 15.9 and 20.3, and a $g-r$ color between 0.375 and 1.908 using extinction-corrected magnitudes and colors). All carbon stars (including dCs) were removed from the sample by SDSS CLASS and SUBCLASS keywords and by matching to all known dCs. Since their binary fraction is likely to be highly biased, we further removed stars originally targeted for reasons of X-ray emission or variability.222We removed from the control sample any eBOSS_TARGET0 stars that are selected by variability as TDSS target (8). Most of these variables are RR Lyr or close eclipsing binaries and some are dC stars. We further removed stars where LEGACY PrimTarget keyword contained ROSAT or where BOSS ANCILLARY_TARGET1 = QSO_VAR, QSO_VAR_LF or QSO_VAR_SDSS. Finally, common proper motion binaries were removed by eliminating control stars with BOSS ANCILLARY_TARGET2 = SPOKE2. Finally, all control stars were required to have a match in the Gaia DR2 data release (Gaia Collaboration et al., 2018). These criteria returned 9,822 stars that had more than one SDSS spectrum for a total of 21,820 spectra.

3.2 Property Matching

To reduce the effects of differing properties between the dC sample and the control sample, we matched the control stars to each dC by finding the normalized “distance” in a “four-property space”: $r$ mag, $g-r$ color, Gaia DR2 total proper motion, and Gaia DR2 parallax.333We use parallax rather than distance due to the subtleties of converting Gaia DR2 parallaxes to distance as noted in Luri et al. (2018)

This distance matching was performed by creating a “normalized coordinate” out of each of the four properties. This coordinate was constructed by subtracting the minimum property value, then dividing by the maximum value for the property. This approach scales all of the values for each property into the range of $[0,1]$ based on the dC sample so that all of the properties are similarly weighted.

These coordinates were used to find the distance from each dC to all of the control stars. These distances, once sorted, give the closest matching control stars for each individual dC based on the chosen four properties. With the control sample sorted for closest matching properties to the dC sample, we drew the closest matches for each dC to create the final, property-matched control sample to analyze along the dC sample. Figure 2, compares histograms of these four properties for the dC and control samples.

3.3 Control Sample Issues

The control sample, even given the matching process we used, is not perfect for several reasons.

(1) The SDSS stellar sample produced by a hodgepodge of different targeting programs, some of which may skew the $\Delta$RV distribution.

(2) It could be more difficult to detect binarity in the control sample because the single spectrum of an unresolved binary contains (by definition) both components. If the two components have significantly different main sequence spectral types or evolutionary stage (e.g., giant + dwarf), then one component is much more luminous than the other — similar to the dC systems we expect, which likely contain a white dwarf too cool to detect in most spectra. However, if the two components have close spectral types (e.g., a K7+M2 binary), they contribute similarly to the spectral flux. Thus, the observed velocity changes are muted because if one component is approaching, the other is receding at any epoch. There are techniques that could mitigate this issue such as attempting to fit the sum of two spectral templates to each spectrum (e.g., as proposed by El-Badry et al. 2018, but this approach would be effective only for some combinations of mass ratios and S/N.

(3) The control sample has a significantly different MJD sampling than that of the dC sample. A majority of the control sample was observed in the earliest versions of the SDSS and have $\Delta$MJDs between spectroscopic epochs on average of 100 days. Most of the dC stars have been specifically targeted by TDSS for repeat spectroscopy during SDSS-IV; so, they have a $\Delta$MJD distribution of typically 1000s of days. While this sampling does affect the range of periods our methods are sensitive to, it should not impact our results. Since we are searching for close binary systems which have large $\Delta$RVs and, therefore, short periods, the control sample’s accessible $\Delta$MJD distribution would only impede our ability to detect wide binary systems for which our sensitivity is already severely limited by the RV errors .

The first two items are observational and may diminish the discriminating power of our tests. Other intrinsic differences may complicate our analysis and interpretation of the results. For instance, we expect dC stars to have a 100% binary fraction, but a very narrow range of companion masses (all white dwarfs, therefore, strongly peaked near 0.5$\,M_{\odot}$). By contrast, the control sample has a certain binary fraction, but the distribution of companion masses in those binaries will have a wider range. The orbital properties of binaries in the control sample may also have a wider range. We expect that the dC has interacted with its (former AGB) companion (e.g., either by wind accretion or Roche lobe overflow) which sets upper (and perhaps even lower) limits on the orbital separation. The only effective limit on orbital separation in the control sample is that the pair be spatially unresolved (\sim$$< $\,2\arcsec$).

4 Radial Velocity Analysis

4.1 Cross-Correlation Method

We measured radial velocity variations ($\Delta$RV) using the IRAF444IRAF is distributed by the National Optical Astronomy Observatories,which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. (Tody, 1993) package FXCOR that cross-correlates between a template and object spectrum following the methods of Tonry & Davis (1979).

Each spectrum was visually inspected to insure the S/N was sufficiently high for cross-correlation as well as to identify wavelength regions with corrupt data. We also searched for any problematic features that could affect the cross-correlation. Those objects that had corrupted regions were marked and individually run through the cross-correlation, ignoring those corrupted regions. The rest of the sample was cross-correlated in a batch, all using the same constraints and regions.

Each epoch combination’s cross-correlation function was manually inspected to check the quality of the cross-correlation. In a small ($\sim$10%) fraction of cases, the cross-correlation function is best fit manually. If the cross-correlation function could not be fit (e.g., no peak in the CCF is preferred), that epoch combination was removed from the sample.

4.2 Cross-Correlation Errors

To validate the cross-correlation process, we ran a variety of tests. The first was to verify and, if possible, minimize the reported errors from FXCOR.

To minimize uncertainties in the $\Delta$RV measurements given by the cross-correlation, we two techniques: (1) direct cross-correlation of one object against itself across different epochs and (2) cross-correlation of each epoch for one object against a SDSS C star template spectrum. For each method we also experimented with changing the regions sampled (e.g., only narrow atomic lines, excluding the carbon bands, or only including carbon bands).

We use one epoch (the early MJD) as the “template” and the other epoch (the later MJD) as the “object”. This method produces some benefits over using the usual template method: (1) This cross-correlation directly provides the $\Delta$RV shift. (3) The SDSS C star template spectrum is for AGB C stars; there are no templates for dC stars.

The second test performed was to determine if the reported values and errors from FXCOR are believable for both dC and control spectra.

FXCOR was generally able to recover the applied shift in both the dC and control samples. However, the reported errors from FXCOR generally are overestimated. By comparing each object’s FXCOR-measured shift for each of the 30 different applied shifts, we determined “multi-shift” errors for each sample as the RMS of the measured $-$ applied shift (see Figure 3).

Figure 3 presents histograms for the reported FXCOR errors and the measured “multi-shift” errors for both the dC and control samples. The top panel (a) shows how across the sample, the errors are smaller for the “multi-shift” errors than those reported by FXCOR. The bottom panel (b) displays that as FXCOR error increases, so do the “multi-shift” errors (a plausible result as spectral S/N is a key factor in the error).

Figure 4 shows the relationship between the brightness ($r$ mag) and the $\Delta$RV errors from the cross-correlation. Brightness (and by proxy S/N) determines the $\Delta$RV errors, and bluer objects tend to have smaller errors (again because these stars are usually brighter and have better S/N).

4.3 dC and Control $\Delta$RVs

The dC and control samples were both cross-correlated using the same method. For every object, each possible combination of epochs was cross-correlated (with the earlier epoch as the template and later epoch as the object). From all possible combinations for an object, we selected the maximum $\Delta$RV for our statistical analysis.

Extremely large $\Delta$RV values (e.g., $>$600 km s*-1*) in a binary with a main sequence component are suspect as in such cases we would expect extremely close orbits and strong signs of interactions and mass transfer. Therefore, any object whose $\Delta$RV was measured to be larger than this value was manually cross-correlated again and had its cross-correlation function manually fit to try and obtain a better $\Delta$RV. If the cross-correlation is unsuitable for fitting, the object was removed from the sample (this only resulted in the removal of two dCs and three control stars).

Figure 5 is a normalized (note the log scale) histogram showing the $\Delta$RV measurements for both the dC (blue) and control (red) samples. The bins are used for each of the samples. This figure demonstrates that both samples have a central core whose width is dominated by the errors. The dC sample, however, has a tail of high $\Delta$RV systems that extends beyond this core. These systems likely represent close binary systems.

In the dC sample, these high $\Delta$RV systems those objects which display $\Delta$RV values $\geq$ 100 km s*-1*. Stars that display such high $\Delta$RV values are indicative of close binary orbits. these large $\Delta$RVs were inspected by shifting the later epoch by the measured $\Delta$RV amount and visually checking to determine if features in the spectrum align.

Given that these systems display no strong signs of interaction (such as explosive variability or an accretion disk continuum emission component), few if any of the dCs are likely to have filled their Roche lobes and be transferring mass to the presumed white dwarf companion. Figure 6 shows the largest possible $\Delta$RV we would expect if the dC was filling its Roche lobe (for varying dC masses and a 0.5M*⊙* WD). Using equation 2 of Eggleton (1983), we calculate the separation for a main-sequence star of each mass to fill its Roche lobe and calculate the corresponding critical $\Delta$RV and period ($\Delta RV_{crit}$, $P_{crit}$ ) that corresponds to the Roche lobe limit. The figure suggests that while we detect dC systems that have large $\Delta$RVs; we have not detected any dCs near the Roche lobe limit .

Figure 7 shows spectra for the dC with one of the largest measured $\Delta$RVs. Both epochs are plotted with the early epoch in black and the late epoch in red. The top panel is of the original spectrum as observed by the SDSS. The bottom panel presents the same epochs, but the later epoch (red) has been shifted by the measured $\Delta$RV = -252 $\pm$ 15 km s*-1* amount. After this shift, the absorption features in the spectrum align confirming this measured $\Delta$RV. All pre-BOSS spectra in this figure have been smoothed by a box-car of 20 pixels, and all later spectra have been smoothed by a box-car of 15 pixels so the SDSS legacy spectra match the resolution of the new BOSS spectrograph; because otherwise, there is a spurious appearance of variability.

Table 1 and Table 2 list the properties for the dC and control sample respectively. Only the first 10 rows for each are shown, the full machine readable tables are available online.

5 Statistical Comparison of $\Delta$RV Distributions

5.1 Anderson-Darling Test

We used a standard, two-sample Anderson-Darling (AD; Scholz & Stephens, 1987) test to determine the similarity between the dC and control $\Delta$RV distributions. From the measured dC and control $\Delta$RVs, the null hypothesis that the dC and control $\Delta$RVs are drawn from the same distribution can be rejected at the 99.95% level ($\log{p}=-3.31$).

5.2 Extreme Deconvolution

A drawback of the AD test is that it does not take measurement uncertainties into account when comparing two distributions. For example, two distributions can look dissimilar if their uncertainties are different even if their true underlying distributions are identical. To ensure that the wider $\Delta$RV observed in our dC sample in Figure 5 is not simply due to differences in the measurement uncertainties , we use the extreme deconvolution (XD) method of Bovy et al. (2011) to deconvolve the underlying distribution of our $\Delta$RV measurements. This XD method employs a Gaussian Mixture Model (GMM) to infer the underlying distribution from a set of heterogeneous, noisy observations or samples while incorporating the errors.

We tested the number of components for the XD-GMM for both the dC and control samples using the Bayesian Information Criterion (BIC) of each model. The BIC approach suggests that the dC sample is best modeled by a mixture of three Gaussians. However, the third Gaussian component for the dC population converges to a small normalization and an unphysically large width; so, we constrain the dC sample to be fit with two components. This decision allows for a central core and for a possible large $\Delta$RV tail that contains close binary systems. Table 3 lists the parameters for these fit components for both the dC and control XD-GMMs.

Figure 8 shows the results of the XD analysis, displaying both XD-GMMs for the two samples (smooth curves) and the histogram of the measured $\Delta$RVs (both the smooth PDFs and histograms have been normalized to an integral of one). This figure demonstrates that both the dC and control samples have a core in their $\Delta$RV distribution, but the dC distribution has a much broader wide component that flares out from the core.

The width of the single component as fit to the control sample is wider than that narrow component of the dC sample. At the risk of over-interpreting this difference, we mention several effects that could contribute to this difference. First, we have used the FXCOR reported errors, which in Section 4.2 were shown to be overestimated. Since the control sample is primarily from legacy SDSS spectra with lower S/N (therefore larger errors), this overestimation is larger and may inflate the error-deconvolved core of the control distribution. Second, the single control sample fit component must accommodate the full range of single and multiple systems. Third, the narrow core of the dC sample could be real; perhaps, some fraction of dC binary orbits have actually widened due to processes related to mass transfer. Chen et al. (2018) report that some wider binaries may undergo Bondi-Hoyle-Littleton mass transfer (Edgar, 2004) and further separate since orbit-synchronized rotation of the giant star could serve as an angular momentum reservoir. However, we would still expect a narrower core for the control sample since a substantial fraction should be single stars.

It should be noted that this XD method simply uses the GMM method with errors to determine the underlying PDF as a mixture of Gaussians, but does not imply or impose any physical meaning or model on our data. However, it clearly allows us to determine that the dC sample has a tail that extends far beyond that of the control sample and is indicative of close binary systems included in the dC population.

6 Binary Orbit Simulations

While our sample lacks sufficient epochs per dC to fit individual orbits, we can use the $\Delta$RV distribution to model the separation distribution. Assuming a primitive dC mass distribution further allows us to characterize the expected period distribution of the dC sample.

Since little is known of dC orbital properties outside of G77-61, we make some physically-driven assumptions. First, we assume that the dC orbits have been circularized ($e=0$) since we expect all of them to have undergone mass transfer. Third, we assume that the WD mass distribution follows that found by Kepler et al. (2007) (i.e., a combination of four Gaussian components. The dominant component is centered on $0.58$ M*⊙* with a width of $0.047$ M*⊙*). We use the distribution for the hot WD sample in Kepler et al. (2007) since they state the distribution for the cooler WDs is not reliable. We also assume a probability density function (PDF) that is uniform over $\cos{i}$ in order to determine the PDF for $\sin{i}$. Finally, since there are no known constraints on the dC mass distribution, we assume a uniform PDF over the range of 0.2$M_{\odot}$ and 1.0$M_{\odot}$, simply assuming that dCs span the same range of masses as normal main sequence stars of the same $g-r$ color distribution.

With these assumed PDFs we simulate a population of stars and sample those orbits to obtain a simulated $\Delta$RV distribution. Comparing the simulated $\Delta$RV distribution to the measured one allows the MCMC to map the separation distribution parameter space.

For the first simulation, we assumed that separations that follow a log-normal distribution with unknown mean $\mu$ and standard deviation $\sigma$, as shown in Equation 1. We placed no constraints on the model parameters, aside from those required by the log-normal PDF (i.e., $\sigma\geq 0.0$ km s*-1* ), and allowed the MCMC walkers to explore the parameter space freely.

[TABLE]

We ran a MCMC for steps with 100 walkers using the Goodman & Weare (2010) algorithm. This approach allowed our walkers to explore all of the parameter space and sample the posterior of our model, which we checked for with the convergence of the chains. Figure 9 shows the resulting MCMC posterior distributions for our model parameters .

The resulting separation distribution from the MCMC has a mean of 0.39 AU, a variance of 0.28 AU, and a median of 0.36 AU. These distances correspond to mean periods of 79-100 days depending on dC mass (G77-61 has a period of 245 days) and a minimum period for this distribution is on order 2.5 days (consistent with Margon et al. (2018), who found a dC with a period of 2.9 days using photometry from the Palomar Transient Factory). The separation distribution generated by our MCMC results in periods that are consistent with the periods known of individual dC systems.

However, de Kool & Green (1995) predicted that dC stars should follow a bimodal period distribution with one peak between $10^{2}-10^{3}$ days and another at $10^{3}-10^{5}$ days. Therefore, we also use our MCMC to model a bimodal mixture model (made of two log-normal separation distributions) of the form in Equation 2.

[TABLE]

In Equation 2, $\mu_{i}$ and $\sigma_{i}$ are the same parameters as in the unimodal distribution, and $\alpha$ is the mixing parameter in this mixture model that controls how much of each distribution contributes to the total PDF. As before, we place no constraints outside of those required by the log-normal PDFs and required by the mixing parameter (i.e. $0.0\leq{}\alpha\leq{1.0}$).

This bimodal distribution MCMC simulation was run for 100,000 steps with 100 walkers. The reduction in steps is required by increased computational load when drawing from this bimodal PDF distribution. While this change does reduce the number of points in the parameter space, the MCMC walkers still mapped the posterior quite well, which we checked for with the convergence of the chains.

Figure 11 shows the MCMC posterior distributions for the bimodal mixture model for all five of our model parameters . In this bimodal mixture model, the total separation distribution has a mean of 0.71 AU and a variance of 1.45 AU. This distribution gives (for the previously stated uniform dC mass range) a range of the mean period of 298-413 days and a minimum period of 1.6 days. Although the number of measured dC periods is quite sparse, the period distribution (calculated from the separation distribution) is in agreement with those few periods in the literature.

One significant improvement can be achieved by measuring dC masses via orbital fits from a follow up spectroscopy campaign. Fitting an orbit to even a few dCs will place initial physical constraints on the dC mass distribution. With a more physical and realistic dC mass distribution the models and MCMC simulations can fit a more accurate separation distribution than can be done with the currently used uniform mass distribution.

7 Balmer Emission Lines

The multi-epoch spectra present an opportunity to survey the dC sample for H$\alpha$ emission line strength and variability. Balmer line emission has been observed in dCs, and Green (2013) found that about 2.6% of dCs showed H$\alpha$ emission.

There are 10 objects with H$\alpha$ emission. Balmer line emission might be expected among dCs for several reasons: (1) coronal emission that may be a result of increased activity from spin-up during the accretion phase of the dC evolution , (2) irradiation of the dC by a hot white dwarf companion, or (3) spin-orbit coupling in a close WD/dC binary.

To explore case 1, in a related effort, we are currently analyzing Chandra observations of a small sample of dC stars to test whether their X-ray emission is consistent with dynamo rejuvenation by accretion spin-up (Green, P.J. et al. 2019, in preparation).

If the emission is from case 2, we expect to detect the WD component in the dC spectra. Indeed, all four of our DA/dCs show emission in their spectra. The remaining six of the H$\alpha$ emission line dCs are of “normal” type (i.e., no visible WD in the spectrum). H$\alpha$ emission is variable in only one normal dC and in none of the DA/dCs. Since close orbits should be involved for cases 2 and 3, we will pursue further multi-epoch spectroscopy for emission line systems.

8 Discussion

Using multi-epoch spectroscopy we have measured the radial velocity variations of a SDSS sample of dC stars. Through MCMC methods, our modeling was able to construct the separation distribution of this dC population that best recreates the observed $\Delta$RV.

We presented the best parameters for two separation models: a unimodal log-normal distribution and a bimodal mixture model of log-normal distributions. Both models result in close binary separation distributions with means less than 1 AU, corresponding to mean periods on the order of 1 year (varying depending on dC mass).

Our sample contains a handful of objects with large ($\geq$ 100 km s*-1*) $\Delta$RV measurements that are indicative of close binary systems. These objects will be targeted for future spectroscopy to constrain orbital parameters thereby better characterizing the separation distribution. In addition, orbital fits will also allow us to determine the masses of the dCs assuming a WD component.

Badenes et al. (2018) analyze the RV variability of main sequence stars and report that the binary fraction is likely higher for more metal-poor stars. Carbon stars are suspected to form more easily at lower metallicity; indeed, about 20% of stars with show carbon-enhancement (e.g., Christlieb et al. 2001; Lee et al. 2013), but that frequency is increasing rapidly as metallicity decreases (Placco et al., 2014). Close binaries ($<10$ AU) also show increases in lower metallicity populations (Moe et al., 2018). The dC in G77-61 is thought to be extremely metal-poor (Gass, 1988). The measured dC $\Delta$RV distribution being wider than the control sample could in part be due both to low metallicity and to evolutionary effects since dC stars are carbon-enhanced by binary mass transfer. If the mass transfer results in inward evolution of the binary, then that should further widen the $\Delta$RV distribution for dCs. Binaries with an AGB primary can be at large separations and still, via wind-Roche lobe overflow, lose orbital angular momentum, evolve towards direct Roche lobe overflow and/or tidal friction towards a common envelope (Chen et al., 2018). Therefore, the fraction of binaries that result in mass transfer in a common envelope and a tight binary configuration may be quite large. The dC stars present a population of post mass transfer binaries that are unusually easy to identify, but may represent just a tiny fraction of such stars — those sufficiently cool and with large enough C/O to produce C2 and/or CN bands. In some cases, the AGB evolution may have been truncated during the common envelope phase before significant carbon dredge up. A much larger space density of post mass transfer M dwarfs may remain unidentified until massive multi-epoch RV surveys become available.

Sloan 2.5-meter

We thank John Ruan for longstanding TDSS contributions, important aspects for this research. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.