© ESO 2021

1. Introduction

The star SDSS J160429.12+100002.2 (J1604+1000 for short) is a rather faint (G = 17.13 mag) blue (G_BP − G_RP = −0.16 mag) star of relatively high Galactic latitude (b = +41.5 deg) that was discovered in a search for white dwarf (WD) stars in the tenth data release (DR10) of the Sloan Digital Sky Survey (SDSS; Ahn et al. 2014) by Gentile Fusillo et al. (2015) and Kepler et al. (2015). While the former classify it as a “narrow-line hydrogen star”, that is, a star with a hydrogen-dominated atmosphere but a surface gravity that is lower than those of typical WDs, the latter report a WD spectral type of a class DAB star that is flagged as uncertain, with effective temperature T_eff = 23 819 ± 270 K and surface gravity log(g) = 5.68 ± 0.03, which is indeed too low for a WD. Because these atmospheric parameters are typical of hot subdwarf stars (sdO/B), Geier et al. (2017) assigned the star to the class of B-type subdwarfs (sdB).

J1604+1000 caught our attention because it exhibits a high negative radial velocity (−345 ± 3 km s⁻¹; Geier et al. 2017). Consequently, the object was selected for spectroscopic follow-up observations in an ongoing project that is looking for hot subdwarfs with extreme radial velocities in order to find high-amplitude radial-velocity variables, that is, short-period close binaries (MUCHFUSS; Geier et al. 2015b) as well as high-speed stars unbound to the Galaxy, such as the hypervelocity hot subdwarf star US 708 (Hirsch et al. 2005, Geier et al. 2015a). The initial survey, the Hyper-MUCHFUSS project, was restricted to sdO/B stars with radial velocities exceeding ±100 km s⁻¹ (Tillich et al. 2011, Németh et al. 2016, Ziegerer et al. 2017).

Based on the SDSS spectrum, which was the only one available at the time of the aforementioned study, a preliminary spectroscopic analysis was carried out using a grid of synthetic spectra that is based on model atmospheres in local thermodynamic equilibrium (LTE) that account for metal-line blanketing (Heber et al. 2000). By means of a χ² minimization technique (Napiwotzki et al. 1999), best fitting atmospheric parameters (T_eff = 16 800 ± 500 K, log(g) = 4.72 ± 0.09, and log(n(He)/n(H)) = −2.2) were found that are, on the one hand, considerably different from the ones derived by Kepler et al. (2015) and, on the other hand, consistent with those of blue horizontal branch (BHB) stars in globular clusters (see, e.g., Moni Bidin et al. 2007; Moehler et al. 2011).

Assuming a typical BHB mass of 0.5 ± 0.1 M_⊙, those preliminary atmospheric parameters gave a spectrophotometric distance d = 6.0 ± 0.9 kpc that, when combined with the ground-based proper motions available back then, resulted in such a high Galactic rest-frame velocity, v_Grf, that the star seemed to be a hypervelocity star, that is, an object that is gravitationally unbound to our Galaxy. This picture did not change when improved proper motions from Gaia Early Data Release 3 (EDR3; Gaia Collaboration 2021; Lindegren et al. 2021b) became available later on, which yielded v_Grf = 740 ± 140 km s⁻¹. Contrary to almost all known hypervelocity stars, the preliminary Galactic orbit of J1604+1000 pointed in our direction, which would hint at the origins of the star lying outside of our Galaxy if, indeed, it were to be a BHB star.

In the meantime, high-quality follow-up spectra were taken during the Hyper-MUCHFUSS project. The implications of the analysis of those spectra are presented here: We demonstrate that J1604+1000 is not a hypervelocity B-type star but instead a short-period single-lined spectroscopic binary system on a highly retrograde Galactic orbit that hosts the progenitor of a low-mass helium white dwarf (He WD) and an unseen companion that is most likely a brown dwarf (BD).

2. Analysis

This section presents the technical details of the spectroscopic analysis (Sect. 2.1), the investigation of the spectral energy distribution (SED; Sect. 2.2), the light-curve analysis (Sect. 2.3), the modeling of the radial-velocity curve (Sect. 2.4), the Bayesian inference of stellar parameters (Sect. 2.5), and the kinematic evaluation (Sect. 2.6).

2.1. Quantitative spectroscopic analysis

Our spectral investigation is based on data from four different instruments. The very first spectrum was taken with the spectrograph of the Baryon Oscillation Spectroscopic Survey (BOSS; Smee et al. 2013) attached to the 2.5 m telescope at Apache Point Observatory. Another five spectra were taken with the Echellette Spectrograph and Imager (ESI; Sheinis et al. 2002) mounted at the Keck II telescope, four with the blue arm of the Intermediate-dispersion Spectrograph and Imaging System¹ (ISIS) at the William Herschel Telescope, and six with the X-shooter (Vernet et al. 2011) at the ESO Very Large Telescope. Table 1 lists the date, exposure time, average signal-to-noise ratio (S/N), and measured radial velocity, v_rad, of the individual observations.

Following the analysis strategy outlined in Irrgang et al. (2014), we simultaneously fit all those spectra over their entire spectral range to determine the atmospheric parameters and elemental abundances. The underlying synthetic models were computed using a series of three codes. The structure of the atmosphere was computed in LTE with ATLAS12 (Kurucz 1996). Based on this atmosphere, population numbers in non-LTE were calculated with the DETAIL code (Giddings 1981; Butler & Giddings 1985), which numerically solves the coupled radiative transfer and statistical equilibrium equations. Using the non-LTE occupation numbers from DETAIL and more detailed line-broadening data as input, the emerging spectrum was eventually computed with the SURFACE code (Giddings 1981; Butler & Giddings 1985). Recent updates to all three codes (see Irrgang et al. 2018b for details) with respect to (i) non-LTE effects on the atmospheric structure, (ii) the implementation of the occupation probability formalism (Hubeny et al. 1994) for hydrogen, and (iii) new Stark broadening tables for hydrogen (Tremblay & Bergeron 2009) and neutral helium (Beauchamp et al. 1997) were also considered.

The results of the spectroscopic analysis are summarized in Table 2. The derived effective temperature, T_eff = 15 840 ± 160 K, shows that the star is of spectral type B. As exemplified in Fig. A.1, all available spectra only exhibit spectral lines of H, He, Mg, Ca, and Fe, which is quite uncommon for “normal” B-type stars. For reference, main sequence stars of similar temperature and solar chemical composition also show lines of C, N, O, Ne, Al, Si, S, and Ar. Owing to the absence of spectral features of those elements, we are only able to provide upper limits for the respective abundances. The measured abundances of He, Mg, Ca², and Fe significantly differ from those of the Sun. In terms of number fractions, He and Mg are under-abundant by more than 1 dex and Ca by 0.7 dex, while Fe is enriched by about 0.8 dex. Upper limits for the abundances of C, N, O, Ne, Al, Si, and S are also significantly below solar (from 0.5 dex for N to 2.5 dex for Si). The only exception from this trend is Ar, for which the upper limit is close to solar. This peculiar abundance pattern clearly hints at ongoing diffusion processes. The inferred value for the microturbulence, $\xi =2.0_{-1.5}^{+1.3}$ km s⁻¹, is rather inconspicuous. In contrast, the radial velocity, v_rad, turns out to be variable on a timescale of hours (see Table 1). Consequently, and as discussed in Sect. 2.4, our value for the projected rotational velocity, $\mathit{v}sin(i)=18.4_{-2.0}^{+2.2}$ km s⁻¹, is overestimated due to the effect of orbital smearing.

2.2. Analysis of the spectral energy distribution

The SED provides an important constraint to cross-check spectroscopic results and to obtain a more comprehensive picture of the star. For this particular object, photometric measurements covering the ultraviolet, optical, and infrared are available. In order to validate our spectroscopic results, we fit the observed SED with synthetic SEDs computed with ATLAS12. Because the surface gravity, microturbulence, He abundance, and metallicity are, if at all, only poorly constrained by photometry, they were set to the values determined from spectroscopy (see Table 2) using Fe as a proxy for metallicity. Consequently, the three quantities T_eff, Θ, and E(44 − 55) remained as free parameters to match the observed SED. The angular diameter, Θ = 2R₁/d (R₁ is the radius and d the distance of the star; see Heber et al. 2018 for more details), was used as a distance scaling factor, and the color excess E(44 − 55) was introduced to account for interstellar reddening, the effect of which was modeled here using the extinction law by Fitzpatrick et al. (2019) with a standard extinction coefficient of R(55) = 3.02. These extinction parameters are analogs of the more widely used color excess E(B − V) and extinction R(V) parameters, but with measurements in the Johnson B and V filters substituted by monochromatic ones at 4400 Å and 5500 Å (see, e.g., Fitzpatrick et al. 2019 for details). In order to empirically account for an apparent infrared excess appearing in the WISE data, we also modeled a blackbody component. This introduced two additional free yet highly anticorrelated parameters, namely a temperature and a flux weighting factor that was parameterized as a surface ratio relative to the stellar component. The results of the fitting procedure are summarized in Table 3 and illustrated in Fig. 1. Effective temperatures from spectroscopy and photometry are basically identical, which corroborates our spectroscopic results. Moreover, the inferred interstellar reddening of E(44 − 55) = 0.058 ± 0.021 mag is consistent with upper limits from reddening maps (Schlegel et al. 1998: E(B − V) = 0.060 ± 0.002 mag; Schlafly & Finkbeiner 2011: E(B − V) = 0.051 ± 0.002 mag). The best fitting parameters of the blackbody component are able to reproduce the observed infrared excess and hint at the presence of a very cool but relatively extended thermal source. Its effective radiation area is ${17}_{-5}^{+16}$ times the projected surface area of the star itself.

Fig. 1. Comparison of synthetic and observed photometry. The top panel shows the SED. The colored data points are filter-averaged fluxes that were converted from observed magnitudes (the respective full widths at tenth maximum of the filters are indicated by the dashed horizontal lines), while the solid gray line represents the best fitting model, i.e., it is based on the parameters from Table 3, degraded to a spectral resolution of 6 Å. The flux is multiplied by the wavelength to the power of three to reduce the steep slope of the SED on such a wide wavelength range. The individual contributions of the stellar (light blue) and blackbody (light red) components are shown as well. The bottom panel shows the residuals, χ, that is, the difference between synthetic and observed magnitudes divided by the corresponding uncertainties. The photometric systems have the following color code: GALEX (violet; Bianchi et al. 2017), SDSS (golden; Alam et al. 2015), Pan-STARRS1 (red; Chambers et al. 2017), Gaia (blue; Riello et al. 2021), and WISE (magenta; Schlafly et al. 2019).

2.3. Analysis of the light curve

The Catalina Sky Survey (CSS; Drake et al. 2009) took a light curve of our target that consists of 379 V-band measurements spread over 3101 days. The mean magnitude of this data set is V_CSS = 17.11 ± 0.01 mag, and the respective standard deviation is 0.08 mag. The latter is very close to the mean of the stated uncertainties on the individual measurements of 0.09 mag, which indicates that this light curve is heavily dominated by noise. Consequently, the available data set only allows an upper limit on the variability of the source to be determined. To this end, we employed a χ² testing scheme to compute p-values for the null hypothesis that the source is variable with a certain amplitude, A. To be independent of any assumptions regarding the form of the underlying light curve, we created one million arbitrarily shaped mock light curves of the form

$$\begin{array}{c}\hfill V_{\mathrm{CSS},i}=V_{\mathrm{CSS}}+A{X}_i=17.11\mathrm{mag}+A{X}_i,\end{array}$$(1)

where 1 ≤ i ≤ 379 is the index of a CSS data point and X_i is a uniformly distributed random number between −1 and 1. Starting from a constant source with A = 0 mag, this Monte Carlo procedure was repeated for steadily increasing values of A until a p-value consistent with 0 was reached, that is, until it became completely unlikely to erroneously reject the null hypothesis despite it being true. The result of this exercise, which is illustrated in Fig. 2, indicates that the available CSS light curve is consistent with the source being photometrically variable with a semi-amplitude of ≲110 mmag.

Fig. 2. P-value inferred from χ2 testing the null hypothesis that the observed CSS light curve exhibits a certain amount of variability, A. The lower the p-value, the more unlikely this null hypothesis is. The width of the line represents 1σ uncertainties derived from modeling the unknown shape of the variability via a Monte Carlo simulation.

2.4. Analysis of the radial-velocity curve

Already in the first follow-up observation run in July 2016, variations in the radial velocity became obvious. Within about an hour, v_rad changed by more than 40 km s⁻¹. In a second run in April 2019, enough follow-up observations were collected to construct the radial-velocity curve. For the modeling of this curve, we did not consider the BOSS spectrum because its exposure time turned out to be roughly half the orbital period, which is too long to infer a representative radial velocity.

Given that it is much smaller than the gap between the observation blocks, the determination of the orbital period, P, was not trivial. To estimate an upper limit for P, we looked at the radial-velocity variations in the night on heliocentric Julian date (HJD) 2 458 565 (see Table 1). Within 2.5 h, v_rad increased from about −376 km s⁻¹ to roughly −317 km s⁻¹. This change in v_rad coincides fairly well with the maximum velocity amplitude inferred from all available observations. For a sinusoidal curve, it would take 1/2 + N times the orbital period to accomplish this maximum change, where N ≥ 0 is an integer. Consequently, P ≈ 2.5 h/(1/2 + N)≲5 h. To account for the various approximations in this reasoning, we adopted a conservative upper limit of 18 h. A lower limit for P can be estimated by considering the effect of orbital smearing, that is, spectral lines are smeared out due to the change in radial velocity over the course of an exposure. Assuming that v_rad increases linearly from its minimum to its maximum value during half the orbital period, it follows that (376 − 317) km s⁻¹/(P/2) times the exposure times of those observations that have sufficient spectral resolution to resolve line broadening (0.42–0.5 h) has to be lower than or equal to the measured value for the projected rotational velocity of $\mathit{v}sin(i)=18.4_{-2.0}^{+2.2}$ km s⁻¹, which yields P ≥ 0.42 h × 2 × 59/20.6 = 2.4 h.

Figure 3 shows the χ² landscape that results from fitting the observed radial velocities with a Keplerian curve for orbital periods between the aforementioned limits. Instead of simply evaluating the Keplerian curve at the midpoint of the exposure times, we averaged it over the exposure times to account for the fact that they are non-negligible fractions of the orbital period. The parameters of the best fitting configuration are listed in Table 4, and the corresponding phased radial-velocity curve is shown in Fig. 4. The star is part of a very close, short-period single-lined spectroscopic binary system with a large negative systemic radial velocity. Although a non-zero eccentricity, e, is formally not excluded by the current data situation, a circularized orbit is most likely from an evolutionary point of view (see Sect. 3.4). A more precise determination of e requires follow-up observations that cover a few consecutive orbits with a single spectrograph in order to avoid the main sources of uncertainty that this study might be suffering.

Fig. 3. χ2 landscape (“periodogram”), which results from fitting the measured radial velocities with a Keplerian model, as a function of the orbital frequency or period. The step size in the orbital frequency was chosen such that phase shifts are always less than 0.01. A minimum filter has been applied to lower the number of points in the plot to a reasonable amount. The dashed red lines mark the best fit.

Fig. 4. Phased radial-velocity curve. The measurements are represented by black symbols with 1σ error bars, and the best fitting Keplerian model is indicated by the solid red curve. The corresponding orbital parameters are given in Table 4. The ESI spectra are shown as crosses, ISIS spectra as open squares, and X-shooter spectra as open diamonds. The dashed horizontal lines indicate the individual exposure times over which the Keplerian curve is averaged before being compared with the measurements. The residuals, χ, are shown in the lower panel. An orbital phase of 0 corresponds to the pericenter passage.

With the orbital parameters from Table 4, we can now turn the tables and check to what extent our estimate for the projected rotational velocity is actually affected by orbital smearing. Assuming again that v_rad increases linearly from its minimum to its maximum value during half the orbital period, the amount of orbital smearing can be estimated to be 4 × K₁/P ≈ 44 km s⁻¹ h⁻¹ times the relevant exposure times (0.42–0.5 h), which yields values between 18–22 km s⁻¹. Consequently, the measured projected rotational velocity of $\mathit{v}sin(i)=18.4_{-2.0}^{+2.2}$ km s⁻¹ is heavily affected by orbital smearing and its true value is probably much lower, which implies that the object is a slow rotator or that we see the star’s rotational axis under a very small inclination, i.

2.5. Bayesian inference of stellar parameters

Based on Bayesian inference, the radius and mass of the visible component (R₁, M₁), the mass of the unseen companion (M₂), and the orbital inclination (i_o), which need not be the same as the star’s inclination (i), can be constrained from the measured angular diameter, Gaia EDR3 parallax, surface gravity, binary mass function, and the geometric restriction that any star in a binary system has to be smaller than its Roche lobe (RL) radius.

To start with, we note that these four parameters are actually not independent of one another but are instead tightly coupled via the following two observational constraints: first, the value for the surface gravity,

$$\begin{array}{c}\hfill g=G{M}_1/R_1^2,\end{array}$$(2)

as derived from spectroscopy (G is the gravitational constant; see Table 2), and second, the value for the mass function,

$$\begin{array}{c}\hfill f=\frac{M_2{sin}^3(i_{\mathrm{o}})}{{(1+M_1/M_2)}^2}={(1-e^2)}^{3/2}\frac{K_1^{3}P}{2\pi G},\end{array}$$(3)

as derived from the analysis of the radial-velocity curve (see Table 4). Consequently, only two of the four quantities of interest have to be treated as free parameters in the statistical analysis. Because plausible prior information is available for i_o and R₁, we chose those two as independent variables and used Eqs. (2) and (3) to (numerically) solve for the masses M₁ and M₂: M₁ = M₁(g, R₁) and M₂ = M₂(f, M₁(g, R₁),i_o) = M₂(f, g, R₁, i_o) (see Fig. 5). Assuming orbital inclinations to be distributed isotropically, the prior for i_o is simply

$$\begin{array}{c}\hfill P(i_{\mathrm{o}})\propto sin(i_{\mathrm{o}}).\end{array}$$(4)

Fig. 5. Mass and RL radius as a function of orbital inclination. Left ordinate (red): mass of the unseen companion, M2, as a function of the orbital inclination, io. A mass of M1 = 0.21 ± 0.01 M⊙ for the visible component is used to numerically solve the binary mass function given in Eq. (3) for M2. Right ordinate (black-rimmed blue-shaded curves): comparison of the RL radius of the visible component at the pericenter passage based on Eq. (10), RRL, p,1 (dark blue), and the respective stellar radius, R1 (light blue), following from M1 and from the spectroscopically inferred surface gravity, $g=G{M}_1/R_1^2$ (see Table 2; G is the gravitational constant). Physical conditions are met when the light blue band is below the dark blue band, that is, when the star is smaller than its RL radius. The widths of all shaded regions cover all involved 1σ uncertainties.

For the radius, R₁, we were able to construct a tailored prior from the measurements of the star’s angular diameter (Θ = 3.82 × 10⁻¹² and σ_Θ = 0.06 × 10⁻¹²; see Table 3) and its Gaia EDR3 parallax (ϖ = 0.117 − Z mas; Z is the parallax bias or “zero point”; σ_ϖ = 0.073 mas). Lindegren et al. (2021b) advise a global parallax correction of Z = −0.017 mas, while Lindegren et al. (2021a) propose a more complex recipe that yields Z = −0.005 mas for this specific target. Because the difference between both values is significantly smaller than σ_ϖ, our results are insensitive to the choice of the correction scheme and we opted for Z = −0.017 mas. Substituting the inverse distance, 1/d, in the equation for the angular diameter, Θ = 2R₁/d, with the parallax ϖ = 1au/d (au is the astronomical unit) and assuming Θ and ϖ to be distributed like Gaussians yields the following integral for the prior probability of R₁:

$$\begin{array}{cc}& P(R_1)\propto {\displaystyle \int \underset{0}{\overset{\infty }{\int }}exp(-\frac{{(\acute{\mathrm{\Theta }}-\mathrm{\Theta })}^2}{2{\sigma }_{\mathrm{\Theta }}^2}-\frac{{(\acute{\varpi }-\varpi )}^2}{2{\sigma }_{\varpi }^2})\delta (R_1-\frac{\acute{\mathrm{\Theta }}\mathrm{au}}{2\acute{\varpi }})}\hfill \\ \hfill & \times P(\acute{\varpi })\mathrm{d}\acute{\mathrm{\Theta }}\mathrm{d}\acute{\varpi },\hfill \end{array}$$(5)

where δ is the Dirac delta function and P(ϖ) is a prior for the parallax. Using a constant value for P(ϖ) would lead to a long tail of unrealistically large radii because the parallax of J1604+1000 is only poorly constrained (σ_ϖ/ϖ ≈ 0.5) and R₁ ∝ ϖ⁻¹. One popular way of coping with this issue is to suppress small parallaxes by applying an exponentially decreasing distance prior as done, for example, by Bailer-Jones et al. (2018). However, as shown, for instance, by de la Fuente Marcos & de la Fuente Marcos (2019), this approach is “not well suited for analyzing individual halo stars when the parallax error is large.” Therefore, in order to circumvent implausible parallaxes for our target, we had to pursue a different strategy. As outlined in Sect. 1, large distances would render the star gravitationally unbound to the Milky Way and thus imply an extragalactic origin. Already for a single star, the probability of observing such a Galactic intruder would be extremely small, and it is even smaller for a binary system because its high space motion would require such violent acceleration scenarios that the binary system would rather be disrupted than ejected. Consequently, we introduced a prior for the parallax that favors Galactic trajectories that are bound to the Milky Way, that is, for which the difference between the current Galactic rest-frame velocity, v_Grf, and the local Galactic escape velocity, v_esc, is negative:

$$\begin{array}{c}\hfill P(\varpi )\propto \{\begin{array}{cc}1,\hfill & v_{\mathrm{Grf}}-v_{\mathrm{esc}}\le 0\hfill \\ exp(-{(v_{\mathrm{Grf}}-v_{\mathrm{esc}})}^2/2{\sigma }_v^2),\hfill & v_{\mathrm{Grf}}-v_{\mathrm{esc}}>0.\hfill \end{array}\end{array}$$(6)

Choosing this functional form and σ_v = 10 km s⁻¹ ensures a smooth and gradual transition between bound and unbound orbits. The quantity v_Grf − v_esc is a function of celestial coordinates, parallax, and proper motions, all of which are available in Gaia EDR3, as well as systemic radial velocity (see Table 4) and the assumed Galactic mass distribution. For the latter, we employed Milky Way mass Model I from Irrgang et al. (2013), which is consistent with Galactic mass estimates based on Gaia DR2 astrometry (see Irrgang et al. 2018a). Uncertainties of and correlations between the other function arguments are accounted for in a Monte Carlo fashion when Eq. (5) is numerically integrated to obtain the resulting prior distribution for R₁, which is shown as a solid gray line in Fig. A.2. The sharp drop at R₁ ≈ 0.46 R_⊙ is caused by Eq. (6) and corresponds to ϖ = Θau/(2R₁) = 3.82 × 10⁻¹² × 1au/(2 × 0.46 R_⊙) = 0.18 mas, which coincides with the 1σ upper limit of the measured parallax. Consequently, the demand for bound trajectories requires the true parallax to be on the high side of what is currently allowed by Gaia.

Independent of the above finding, there is another aspect that calls for even smaller stellar radii and thus turns out to be the main driver for high parallaxes: the simple geometric constraint that any star in a binary system has to be smaller than its RL radius. Owing to the compactness of the binary system, this requirement plays an important role because it renders low orbital inclinations and large stellar radii less likely. This assertion is illustrated in Fig. 5 for an exemplary mass of M₁ = 0.21 ± 0.01 M_⊙, but it is also valid for other masses. To be as general as possible, we accounted for the slight nominal eccentricity, e, of the binary orbit by focusing in the following on the orbital phase at which the RL is smallest, that is, on the periastron passage. The expression for the projected periastron distance is

$$\begin{array}{cc}& r_{\mathrm{p}}sin(i_{\mathrm{o}})=(1-e)(a_1+a_2)sin(i_{\mathrm{o}})\hfill \\ \hfill & =(1-e)a_{1}sin(i_{\mathrm{o}})(1+a_2/a_1)\hfill \\ \hfill & =(1-e)a_{1}sin(i_{\mathrm{o}})(1+M_1/M_2),\hfill \end{array}$$(7)

where a₁ and a₂ are the components’ semimajor axes and

$$\begin{array}{c}\hfill (1-e)a_{1}sin(i_{\mathrm{o}})=(1-e){(1-e^2)}^{1/2}\frac{K_{1}P}{2\pi }=0.102\pm 0.008R_{\odot }\end{array}$$(8)

is given from the analysis of the radial-velocity curve (see Table 4). Using the shorthand notation q ≔ M₁/M₂ and an approximation formula for the relative RL radius (Eggleton 1983),

$$\begin{array}{c}\hfill r_{\mathrm{RL}}(q)=\frac{0.49q^{2/3}}{0.6q^{2/3}+ln(1+q^{1/3})},\end{array}$$(9)

the RL radius at periastron can be written as

$$\begin{array}{cc}& R_{\mathrm{RL},\mathrm{p},1}(f,g,R_1,i_{\mathrm{o}})=r_{\mathrm{RL}}(q)r_{\mathrm{p}}=r_{\mathrm{RL}}(q)\frac{r_{\mathrm{p}}sin(i_{\mathrm{o}})}{sin(i_{\mathrm{o}})}\hfill \\ \hfill & =r_{\mathrm{RL}}(q)\frac{(1-e)a_{1}sin(i_{\mathrm{o}})(1+q)}{sin(i_{\mathrm{o}})}.\hfill \end{array}$$(10)

The dependences on the mass function, f, the surface gravity, g, and the radius, R₁, enter via q = M₁/M₂, M₁ = M₁(g, R₁), and M₂ = M₂(f, g, R₁, i_o) (see above). Consequently, the geometric constraint that the stellar radius must not exceed the RL radius, R₁ ≤ x R_RL, p,1, yields

$$\begin{array}{c}\hfill h(f,g,R_1,i_{\mathrm{o}})\frac{R_{1}sin(i_{\mathrm{o}})}{r_{\mathrm{RL}}(q)(1+q)x}\le (1-e)a_{1}sin(i_{\mathrm{o}}).\end{array}$$(11)

The factor 0 < x ≤ 1 gives the fraction of the RL that the star is allowed to fill and is equal to 1 if only the current geometry of the binary system is considered. However, as discussed in Sect. 3.2, such an approach would ignore the evolutionary history of the stellar components, which hints at x < 1. To understand in detail how the Bayesian inference is affected by x, we did not marginalize our results over this parameter but explicitly explored different values for it in the following. This has the advantage that we do not bias the outcome by any prior assumption on x, which would be very vague at best.

The goal of the statistical analysis is to find those values of i_o and R₁ – and thus also of M₁ and M₂ – that have the highest probability of obeying the observational constraint given by Eq. (11). Hence, we wished to determine the probability of having a set of parameters {R₁, i_o} given the observation h, that is, we had to find an expression for the conditional probability P({R₁, i_o}|h). Using Bayes’ theorem, we can write

$$\begin{array}{c}\hfill P(\{R_1,i_{\mathrm{o}}\}|h)=CP(h|\{R_1,i_{\mathrm{o}}\})P(R_1)P(i_{\mathrm{o}}),\end{array}$$(12)

where C is a normalization constant and P(i_o) and P(R₁) are the priors over the values of the two independent parameters that are given in Eqs. (4) and (5). The constraint from Eq. (11) is implemented via the following likelihood function:

$$\begin{array}{c}\hfill P(h|\{R_1,i_{\mathrm{o}}\})\propto \{\begin{array}{cc}1,\hfill & h\le \overline{h}\hfill \\ exp(-{(h-\overline{h})}^2/2{\sigma }_h^2),\hfill & h>\overline{h}.\hfill \end{array}\end{array}$$(13)

The quantities $\overline{h}=0.102R_{\odot }$ and σ_h = 0.008 R_⊙ are derived from the observed radial-velocity curve (see Table 4). With all expressions in Eq. (12) defined, it is possible to compute the posterior probability distribution functions (PDFs) for the two independent variables by marginalizing out the other one. To this end, we applied a Monte Carlo procedure for the numerical integration, which allowed us to propagate the uncertainties in the surface gravity, g, and in the mass function, f. Posterior PDFs for various derived quantities, such as M₁ and M₂, were also computed.

The result of this exercise for RL filling factors 0.5 ≤ x ≤ 1 is summarized in Table 5 and illustrated in Fig. A.2. Scenarios with x < 0.5 are not discussed here because they predict stellar radii that are highly unlikely given the prior information on R₁. Interestingly, the mass ratio, the semimajor axis, and the RL radii of the two components are almost insensitive to the assumed value for x. Moreover, independent of x, the unseen companion is most likely not massive enough to ignite hydrogen burning in its core and therefore qualifies as a candidate BD. In contrast, the inferred parameters of the visible component are sensitive to the precise value of x. Here we first consider the case where the star is allowed to fill its entire current RL at periastron, that is, x = 1. As demonstrated in Fig. A.2, the prior and posterior PDFs of i_o and R₁ are then almost identical, which shows that the geometric constraint of Eq. (11) barely affects the outcome in this case. This is expected when looking at Fig. 5: Requiring the RL radius to exceed the stellar radius only renders small orbital inclinations and large stellar radii a little less likely, the latter of which are in any case already unlikely by virtue of Eq. (6). Consequently, R₁ and thus also M₁ are essentially prescribed by the prior on R₁ (see Eq. (5)). We now turn to cases with x < 1. As demonstrated in Fig. A.2, this requirement puts increasingly strong constraints on the four parameters of interest. Generally speaking, the smaller x is, the smaller M₁, M₂, and R₁ are and the larger i_o is, which is not surprising because it simply reflects the fact that a smaller allowed volume necessitates a smaller object. In principle, the light curve could be exploited to constrain x, and thus the precise parameters of the visible component, solely from observations. To demonstrate this, we calculated the expected V-band semi-amplitudes caused by the ellipsoidal deformation of the visible component, A_elli, 1, and by irradiance effects, A_irra, following Eqs. (3) and (6) in Morris & Naftilan (1993). A gravity-darkening coefficient of τ₁ = 0.51 (see Eq. (10) in Morris 1985) and a limb-darkening coefficient of u₁ = 0.34 (Claret & Bloemen 2011) were used to compute A_elli, 1. To calculate a rough upper limit for A_irra, the unknown stellar radius and effective temperature of the unseen companion were approximated by their highest plausible values, that is, by R_RL, p,2 and by assuming a very hot BD of 3000 K (Baraffe et al. 2003), respectively. Furthermore, blackbody behavior was assumed to compute the magnitude difference between the two components, which enters the formula for A_irra. The resulting numbers, which are listed in Table 5, show that ellipsoidal deformation is the dominating effect and that its strength decreases from about 13 mmag to 5 mmag when going from x = 1 to x = 0.5. While the currently available light curve does not allow us to resolve semi-amplitudes below ∼110 mmag (see Sect. 2.3), future follow-up observations might be good enough to measure those subtle variations with the required precision. For the time being, we have to rely on aspects from stellar evolution theory to narrow down x and thus argue that J1604+1000 hosts a proto-He WD (see Sect. 3.2).

2.6. Kinematic analysis

The unusually large negative systemic radial velocity (see Table 4) of this binary system calls for a closer inspection of its kinematic properties. Using the spectroscopic distance and the correction for gravitational redshift from Table 5 (for an exemplary case of x = 0.8), as well as proper motions from Gaia EDR3, which seem to be reliable because the renormalized unit weight error (RUWE; see Lindegren 2018) indicates a well-behaved astrometric solution (RUWE = 0.99), allows the current position and velocity vector of the system to be computed. Based on the object’s location in the Toomre diagram (Fig. 6), we conclude that it is on a highly retrograde, halo-like orbit that is more extreme than that of similar binary systems. The fact that the respective difference between the current Galactic rest-frame velocity and the local Galactic escape velocity, ${\mathit{v}}_{\mathrm{Grf}}-{\mathit{v}}_{\mathrm{esc}}=-{150}_{-80}^{+120}$ km s⁻¹, is well below zero illustrates that the condition x = 0.8 is already more than enough to make the system gravitationally bound to the Milky Way, even without imposing this condition via Eq. (6).

Fig. 6. Position of the star in the Toomre diagram for x = 0.8. The quantity V is the velocity component in the direction of Galactic rotation, U is in the direction towards the Galactic center, and W is perpendicular to the Galactic plane. The star, the Sun, and the local standard of rest (LSR) are marked by a red cross with 1σ error bars, a yellow circled dot (⊙), and a black plus sign (+), respectively. The dashed gray circles centered around the LSR represent boundaries for thin (radius of 85 km s−1) and thick (radius of 180 km s−1) disks following Fuhrmann (2004). A reference sample of 75 binaries that host extremely low-mass He WDs is shown as blue data points, with light blue marking objects that orbit in the Galactic disk (Brown et al. 2020).

3. Discussion

3.1. The nature of the unseen component

The companion of the visible B-type star does not exhibit signatures in the optical spectra. Despite this lack of a direct hint, the available constraints clearly indicate that the object is of very low mass (see Table 5), most probably below the metallicity-dependent hydrogen burning limit of ∼0.07–0.094 M_⊙ (see, e.g., Dieterich et al. 2014 and references therein). Consequently, it is a BD candidate. This is primarily a result of the system’s low mass function (see Table 4) and the prior distribution for i_o, which, according to Eq. (4), favors high inclinations and consequently low companion masses (see Fig. 5). For the M₂ values listed in Table 5, models by Baraffe et al. (2003) predict BD radii of 0.12–0.17 R_⊙ for an age of 0.1 Gyr and 0.09–0.10 R_⊙ for an age of 10 Gyr, all of which are smaller than the corresponding RL radius R_RL, p,2, which is in the range 0.17–0.23 R_⊙ (see Table 5). So even without invoking any restrictions on the dimension of the companion, a consistent picture is obtained when assuming a BD nature for the unseen component.

In principle, the SED could be exploited to further investigate the nature of the low-mass companion, substellar or not, owing to its expected low temperature (∼1500–3000 K) that makes the object shine in the infrared. However, even when assuming the lowest R₁ value from Table 5, $R_1=0.{27}_{-0.06}^{+0.07}{R}_{\odot }$ for x = 0.5, the maximum possible relative radiation area of the unseen component is ${(R_{\mathrm{RL},\mathrm{p},2}/R_1)}^2=0.{38}_{-0.18}^{+0.38}$, which is almost two orders of magnitude lower than what is obtained from modeling the observed infrared excess with a blackbody component (${17}_{-5}^{+16}$; see Table 3). We therefore conclude that the observed infrared excess is much too strong to be explained by the low-mass component, which also renders the prospects of obtaining insights about the companion from SED modeling dim. Instead, optical light curves with higher cadence and precision than what is currently available are eventually necessary in order to search for potential signatures of the companion, such as ellipsoidal modulation or an irradiance effect, which are estimated to be in the millimagnitude regime (see Table 5). The unexpected presence of an infrared excess may be explained as a relic of a previous mass transfer episode (see Sect. 3.4).

3.2. The nature of the visible component

The analysis of the follow-up spectra yielded improved atmospheric parameters with respect to the preliminary analysis outlined in Sect. 1. The effective temperature turned out to be lower and the surface gravity higher. Figure 7 shows the preliminary and revised position of J1604+1000 in the Kiel diagram, demonstrating that the new parameters are no longer consistent with those of a BHB star because they place the star well below the zero-age horizontal branch (ZAHB). Instead, the comparison with evolutionary tracks for stripped helium cores suggests that the star is a proto-He WD of 0.21 ± 0.01 M_⊙, which, according to Table 5, would correspond to 0.5 ≲ x ≲ 0.6. However, as outlined in Sect. 3.4, the application of those tracks has one major caveat for this particular object, which is why we refrain from using them to determine the stellar parameters.

Fig. 7. Preliminary and revised position of J1604+1000 (light red and dark red 1σ error bars) in the Kiel diagram. The black lines are evolutionary tracks for stripped helium cores from Driebe et al. (1998) labeled with their respective masses (stars evolve from top to bottom). For reference, the locus of the zero-age horizontal branch (ZAHB) for [Fe/H]= − 1.48 (Dorman et al. 1993) is shown as a dashed gray line. The meaning of the (light) blue data points is the same as in Fig. 6. The six proto-He WDs that are mentioned in the text – PSR J1816+4510, EVR-CB-001, GALEX 1717+6757, HD 188112, EL CVn, and WASP 0247−25 B – are presented in black.

The evolutionary histories of BHB stars and He WDs are very similar. The horizontal branch (HB) is a sequence of core-helium burning stars that formed after the ignition of helium burning at the tip of the red giant branch (RGB). The morphology of the HB is determined by the mass of the hydrogen envelope: The lower it is, the bluer the HB star appears. Hence, in order to populate the BHB, a considerable fraction of the envelope mass must have been lost during the transition from the RGB to the HB. However, it is conceivable that mass loss occurs even before the progenitor star has reached the tip of the RGB and ignited helium burning in the core. The outcomes of such an early mass loss are stripped stars that cool down and become low-mass He WDs. If they are less massive than 0.3 M_⊙, they are called extremely low-mass (ELM) WDs. The high mass loss that is necessary to strip off the envelope and eventually form He WDs is most likely due to binary mass transfer (Driebe et al. 1998, Althaus et al. 2013, Istrate et al. 2016). In fact, most ELM WDs are found in short-period double-degenerate systems (see Brown et al. 2020 and references therein).

Although the surface gravity of J1604+1000 is currently not high enough to classify it as a WD, its low mass and its status as short-period binary system make it a prime candidate for being the progenitor of a low-mass He WD. This assumption is strengthened by noting that the atmospheric parameters of our star are almost identical to those of the low-mass proto-He WD companion of the millisecond pulsar PSR J1816+4510 (T_eff = 16 000 ± 500 K, log(g) = 4.9 ± 0.3; Kaplan et al. 2013). Another He WD progenitor with similarly low surface gravity but slightly higher temperature is EVR-CB-001 (Ratzloff et al. 2019), which is a binary system composed of a pre-ELM WD (T_eff = 18 500 ± 500 K, log(g) = 4.96 ± 0.04) and an unseen low-mass He WD companion. Finding such unevolved proto-He WDs offers the opportunity to study a key intermediate stage in their evolution (Ratzloff et al. 2019).

Assuming the visible component to be a proto-He WD now provides an argument for why to investigate RL filling factors x < 1 (see Sect. 2.5). Proto-He WDs generally contract with time (see, e.g., evolutionary tracks by Driebe et al. 1998, Althaus et al. 2013, and Istrate et al. 2014a, 2016), and it could well be that the stripping event that formed the proto-He WD did not occur recently but rather several million years ago, giving the stripped star time to shrink by some fraction in the meantime. Because the typical timescales for changes in the binary orbit by the emission of gravitational waves or by orbital interactions with a potential third body, such as a shell or disk, are much longer, we can assume that the orbital configuration has remained more or less unaltered since the stripping event. We may thus not only require that the current stellar radius not exceed the RL radius but also that a somewhat larger object had to have fit into the current orbital geometry in the past. This is equivalent to saying that the current shrunken stellar radius has to be smaller than an unknown fraction, x, of the RL radius.

3.3. Standard formation scenarios for He WDs

Typically, (proto-)He WDs are found in one of the three following systems: double-degenerate systems (see, e.g., Brown et al. 2020 and references therein), where both components are WDs; in millisecond pulsar systems (see, e.g., Istrate et al. 2014b and references therein), where the (proto-)WD orbits a fast spinning neutron star; and in EL CVn binaries (Maxted et al. 2014), where the proto-He WD is the under-luminous B-type star in an eclipsing binary systems with A- or F-type main sequence primaries (see, e.g., Chen et al. 2017, Wang et al. 2020a and references therein). While millisecond pulsar systems and EL CVn binaries result from mass transfer onto the companion through stable RL overflow, the formation of double-degenerate systems depends on the mass ratio of the accreting WD and the donor star, which usually is an RGB star of ∼1–2 M_⊙. If the WD accretor is massive enough to yield a mass ratio that is sufficiently close to unity, mass will be transferred via stable RL overflow. In contrast, if the mass ratio is significantly different from unity, mass transfer will not be stable and a common-envelope (CE) phase will occur (see, e.g., Istrate et al. 2016 and references therein for more details).

Given that CE evolution is one of the most complex and uncertain aspects in stellar evolution theory, it is not surprising that most theoretical studies about the evolution of low-mass He WDs focus on the RL overflow channel, either by explicitly following the evolution of a binary (e.g., Althaus et al. 2013; Istrate et al. 2014a, 2016) or by artificially removing the envelope during the calculation of single-star evolutionary models (see, e.g., Driebe et al. 1998). As long as the initial model for the stripped helium core is correct, the latter approach is justified by noting that the further evolution of the model is independent of the details of the previous mass loss episode (Driebe et al. 1998). Consequently, the two channels, RL overflow and CE ejection, may yield very similar evolutionary tracks in the Kiel diagram, provided that the post-stripping structure is indeed similar. A similar post-stripping structure is, for example, assumed in the binary population synthesis study by Li et al. (2019), who discuss the characteristics of double degenerates resulting from the two formation channels.

However, a key aspect for the further evolution of a recently stripped proto-He WD is M_H, the mass of its hydrogen envelope (see, e.g., Sarna et al. 2000). Calculations with RL overflow – and thus almost all existing evolutionary tracks – robustly predict relatively thick hydrogen envelopes (M_H ∼ 10⁻⁴–10⁻² M_⊙; see, e.g., Driebe et al. 1998, Althaus et al. 2013, Istrate et al. 2016). For CE evolution, the situation is more complex because the details of this dynamical process are far from being understood, which is why M_H remains an unknown parameter in that case. According to Calcaferro et al. (2018), the formation of low-mass He WDs with thin hydrogen envelopes (M_H ∼ 10⁻⁶ M_⊙) can then not be discarded, although strong observational evidence for that is currently missing. The proposed smoking gun for this scenario is the discovery of ELM WDs with effective temperatures below 7000 K (Calcaferro et al. 2018). In the following, we suggest that J1604+1000 has formed with a thin hydrogen envelope, probably providing the as yet missing evidence for this hypothesis.

3.4. The need for a thin hydrogen envelope

Based on the discussions in the previous sections, we propose the following formation scenario for J1604+1000, which has to be fine-tuned by tailored binary evolution calculations that are beyond the scope of this discovery paper. The progenitor system is composed of a low-mass (∼1–2 M_⊙) main sequence star plus a probably substellar companion of much lower mass (≲0.08 M_⊙). Triggered by stellar evolution, the primary expands and, at some point, fills its RL, which leads to the formation of a CE owing to the extreme mass ratio. Friction forces cause the system to spiral in and to finally eject the CE, leaving behind a close binary system that is surrounded by the expanding stripped envelope. The signature of a cool extended source observed as infrared excess in the analysis of the SED (see Sect. 2.2) may then be caused by this circumstellar material. A similar scenario has been suggested for the formation of WDs with BD companions (see Parsons et al. 2017 and references therein).

The star’s position in the Kiel diagram combined with evolutionary tracks with thick hydrogen envelopes (Driebe et al. 1998, Althaus et al. 2013, Istrate et al. 2014b) hints at a mass for the stripped helium core of about 0.21 ± 0.01 M_⊙ (see Fig. 7). However, those tracks also predict log(g)∼3 right after stripping, which would imply that the radius of the exposed core has decreased by a factor of ∼10^{(4.86 − 3)/2} = 8.5 since the stripping event. Because a 0.21 ± 0.01 M_⊙ star with log(g) = 4.86 ± 0.04 is already close to filling its RL (see Fig. 5), it is impossible to make such an object considerably larger in the past, which rules out this scenario. As demonstrated by Calcaferro et al. (2018), a reduction in the thickness of the hydrogen envelope causes the envelope to become denser, which leads to a smaller stellar radius and thus to a larger surface gravity. Consequently, evolutionary tracks with thin hydrogen envelopes may resolve this contradiction and provide a spectroscopic mass estimate in a fully consistent picture.

3.5. Abundance studies of (proto-)He WDs

Regardless of the precise mass of J1604+1000, it is one of the very rare cases that have allowed for a comprehensive abundance analysis of a proto-He WD. Little information on the chemical composition of other (proto-)He WDs is currently available, which remains a crucial barrier to the advancement of our understanding of these objects (see, e.g., Istrate et al. 2016). For the companion of the abovementioned millisecond pulsar PSR J1816+4510, Kaplan et al. (2013) determined super-solar abundances for He, Na, Mg, Si, Ca, and Fe, which, apart from Fe, is completely contrary to what we find here. This is remarkable given that the atmospheric parameters of the two visible components are almost identical. For the ELM WD SDSS J074511.56+194926.5 (T_eff = 8380 ± 120 K, log(g) = 6.21 ± 0.07), Gianninas et al. (2014) measured solar abundances for Mg, Ca, Ti, Cr, and Fe. The most comprehensive abundance studies available are based on ultraviolet spectra of GALEX J1717+6757 (Vennes et al. 2011, Hermes et al. 2014), a pre-ELM WD of similar temperature (14 900  ±  200 K) but higher surface gravity (5.67 ± 0.05 dex) than our star, and the sdB-type pre-ELM WD HD 188112 (T_eff = 21 500  ±  500 K, log(g) = 5.66  ±  0.06; Heber et al. 2003, Latour et al. 2016). Hermes et al. (2014) derived abundances of C, Al, Si, P, S, Ca, Ti, Cr, and Fe as well as upper limits for N, O, Mg, Sc, and Ni, while Latour et al. (2016) determined abundances of Mg, Al, Si, P, S, Ca, Ti, Cr, Mn, Fe, Ni, Zn, Ga, Sn, and Pb as well as upper limits for C, N, and O. The resulting abundance patterns are thought to be produced by an interplay of atomic diffusion and rotational mixing (Hermes et al. 2014, Latour et al. 2016, Istrate et al. 2016). Also, based on ultraviolet spectroscopy, Wang et al. (2020b) recently found evidence for atomic diffusion in the atmosphere of the pre-ELM WD EL CVn (T_eff = 11 890 ± 490 K, log(g) = 4.77 ± 0.02). Istrate et al. (2017) and Heuser (2018) presented an abundance analysis of the EL CVn-type star WASP 0247−25 B (T_eff = 10 870 ± 230 K, log(g) = 4.70 ± 0.12), which turned out to be somewhat rich in He and deficient in O, Mg, Si, Ca, Ti, and Fe, though these deficiencies are less pronounced than those for J1604+1000. Finally, based on a visual inspection of spectra of the ELM survey (Brown et al. 2020), Hermes et al. (2014) concluded that all ELM WDs show Ca in their optical spectra if their surface gravity is lower than log(g) = 5.9.

4. Summary and outlook

Based on a comprehensive investigation that utilizes multi-epoch data from astrometry, photometry, and optical spectroscopy, we infer that J1604+1000 is a short-period (P ≈ 3.4 h) single-lined spectroscopic binary system that contains a visible B-type star (M₁ ≲ 0.52 M_⊙) and an unseen low-mass companion (M₂ ≲ 0.08 M_⊙) that is most likely a BD. The masses of the two components were constrained via Bayesian inference from the measured angular diameter, Gaia EDR3 parallax, surface gravity, binary mass function, and the condition that the stellar radius has to be smaller than its respective RL radius. The combination of effective temperature, T_eff = 15 840 ± 160 K, surface gravity, log(g) = 4.86 ± 0.04, and stellar mass suggests that the visible component is a proto-He WD. It exhibits a peculiar abundance pattern that is indicative of ongoing atomic diffusion processes. In terms of number fractions, He and Mg are below solar by more than 1 dex and Ca by 0.7 dex, while Fe is above solar by about 0.8 dex. Although other chemical species do not exhibit spectral lines in the available optical spectra, upper abundance limits for C, N, O, Ne, Al, Si, S, and Ar can be derived; these elements are also significantly under-abundant (from 0.5 dex for N to 2.5 dex for Si) with the exception of Ar, for which the limit is close to solar. Besides providing an estimate for the angular diameter, Θ = 2R₁/d = (3.82 ± 0.06)×10⁻¹² (R₁ is the radius of the visible component and d the distance to the binary system), the investigation of the SED also reveals an infrared excess that can be empirically modeled by a blackbody component with temperature $T_{\mathrm{bb}}={2300}_{-600}^{+400}$ K and an effective radiation area that is ${17}_{-5}^{+16}$ times the projected surface area of the visible star.

All observational constraints can be consistently explained by proposing that the system is a post-CE binary that probably emerged relatively recently from the mass transfer phase. Its extremely high Galactic rest-frame velocity indicates that it belongs to a very old (halo) population. Therefore, the BD is probably compact, which favored its prospects to survive its engulfment in the envelope of the red giant progenitor. Tighter constraints on the stellar parameters of the two components can be expected from more precise and accurate parallax measurements from future Gaia data releases as well as from optical light curves that are at least ten times more precise than the currently available CSS V-band light curve; this would allow the effect of ellipsoidal deformation to be measured. Moreover, owing to the compactness of the binary and the inferred constraints on the orbital inclination, the system is very likely eclipsing, with primary eclipse depths of roughly 10%.

The reported discovery is particularly interesting for two reasons. First, the derived chemical abundance pattern of the proto-He WD shows heavy signatures of ongoing atomic diffusion processes and may thus help us to better understand the details of those complex processes, especially because quantitative abundance studies of proto-He WDs have been very rare so far. High-resolution ultraviolet spectroscopy would allow for abundance determinations of many more chemical species than what is currently possible based on optical spectra. Second, standard evolutionary tracks for stripped helium cores predict post-stripping stellar radii that are too large to be consistent with the current orbital configuration of J1604+1000, thus making the object a promising test bed for theory. Motivated by the RL overflow scenario, all of those standard evolutionary tracks are computed for thick hydrogen envelopes. However, this approach may not be valid for systems that underwent a CE phase, as we propose for J1604+1000. Following Calcaferro et al. (2018), we argue that models with thin hydrogen envelopes would probably resolve the abovementioned inconsistency because they predict smaller stellar radii and thus larger surface gravities. Consequently, the thickness of the hydrogen envelope should be considered as an additional parameter in the spectroscopic mass and age determinations of (proto-)He WDs in post-CE systems.

Acknowledgments

We thank the anonymous referee for a very careful and constructive report. A. I. and U. H. acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG) through grants IR190/1-1 and HE1356/71-1. S. G. was supported by the Heisenberg program of the DFG through grants GE 2506/8-1 and GE 2506/9-1. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. We thank John E. Davis for the development of the SLXFIG module used to prepare the figures in this paper. Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 0102.D-0092(A). The William Herschel Telescope is operated on the island of La Palma by the Isaac Newton Group of Telescopes in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. The ISIS spectroscopy was obtained as part of W17AN011. Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the US Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University. The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE). This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. The CSS survey is funded by the National Aeronautics and Space Administration under Grant No. NNG05GF22G issued through the Science Mission Directorate Near-Earth Objects Observations Program. The CRTS survey is supported by the US National Science Foundation under grants AST-0909182 and AST-1313422.

References

Appendix A: Additional figures

Fig. A.1. Exemplary comparison of the best fitting model spectrum (red line) with normalized observation (black line; X-shooter spectrum taken on March 23, 2019, i.e., 2 458 565.8693 HJD). Light colors mark regions that have been excluded from fitting, e.g., due to the presence of interstellar or telluric lines. Residuals, χ, are shown as well. The optical spectrum only exhibits lines of H, He, Mg, Ca, and Fe.

Fig. A.1. continued.

Fig. A.2. Corner plot visualizing the correlations between the orbital inclination (io), the radius and mass of the visible component (R1, M1), and the mass of its unseen companion (M2) for an RL filling factor x = 1 (see Eq. (11)). The diagonal panels show the respective posterior PDFs. Maximum values are marked by dashed red lines, and lower and upper bounds of the highest density interval are shown as dotted red lines. Prior PDFs for io and R1 according to Eqs. (4) and (5) are drawn as solid gray lines.

Fig. A.2. continued for x = 0.9.

Fig. A.2. continued for x = 0.8.

Fig. A.2. continued for x = 0.7.

Fig. A.2. continued for x = 0.6.

Fig. A.2. continued for x = 0.5.

All Tables

All Figures

Fig. 1. Comparison of synthetic and observed photometry. The top panel shows the SED. The colored data points are filter-averaged fluxes that were converted from observed magnitudes (the respective full widths at tenth maximum of the filters are indicated by the dashed horizontal lines), while the solid gray line represents the best fitting model, i.e., it is based on the parameters from Table 3, degraded to a spectral resolution of 6 Å. The flux is multiplied by the wavelength to the power of three to reduce the steep slope of the SED on such a wide wavelength range. The individual contributions of the stellar (light blue) and blackbody (light red) components are shown as well. The bottom panel shows the residuals, χ, that is, the difference between synthetic and observed magnitudes divided by the corresponding uncertainties. The photometric systems have the following color code: GALEX (violet; Bianchi et al. 2017), SDSS (golden; Alam et al. 2015), Pan-STARRS1 (red; Chambers et al. 2017), Gaia (blue; Riello et al. 2021), and WISE (magenta; Schlafly et al. 2019).

In the text

	Fig. 2. P-value inferred from χ2 testing the null hypothesis that the observed CSS light curve exhibits a certain amount of variability, A. The lower the p-value, the more unlikely this null hypothesis is. The width of the line represents 1σ uncertainties derived from modeling the unknown shape of the variability via a Monte Carlo simulation.
In the text

Fig. 3. χ2 landscape (“periodogram”), which results from fitting the measured radial velocities with a Keplerian model, as a function of the orbital frequency or period. The step size in the orbital frequency was chosen such that phase shifts are always less than 0.01. A minimum filter has been applied to lower the number of points in the plot to a reasonable amount. The dashed red lines mark the best fit.

In the text

Fig. 4. Phased radial-velocity curve. The measurements are represented by black symbols with 1σ error bars, and the best fitting Keplerian model is indicated by the solid red curve. The corresponding orbital parameters are given in Table 4. The ESI spectra are shown as crosses, ISIS spectra as open squares, and X-shooter spectra as open diamonds. The dashed horizontal lines indicate the individual exposure times over which the Keplerian curve is averaged before being compared with the measurements. The residuals, χ, are shown in the lower panel. An orbital phase of 0 corresponds to the pericenter passage.

In the text

Fig. 5. Mass and RL radius as a function of orbital inclination. Left ordinate (red): mass of the unseen companion, M2, as a function of the orbital inclination, io. A mass of M1 = 0.21 ± 0.01 M⊙ for the visible component is used to numerically solve the binary mass function given in Eq. (3) for M2. Right ordinate (black-rimmed blue-shaded curves): comparison of the RL radius of the visible component at the pericenter passage based on Eq. (10), RRL, p,1 (dark blue), and the respective stellar radius, R1 (light blue), following from M1 and from the spectroscopically inferred surface gravity, $g=G{M}_1/R_1^2$ (see Table 2; G is the gravitational constant). Physical conditions are met when the light blue band is below the dark blue band, that is, when the star is smaller than its RL radius. The widths of all shaded regions cover all involved 1σ uncertainties.

In the text

Fig. 6. Position of the star in the Toomre diagram for x = 0.8. The quantity V is the velocity component in the direction of Galactic rotation, U is in the direction towards the Galactic center, and W is perpendicular to the Galactic plane. The star, the Sun, and the local standard of rest (LSR) are marked by a red cross with 1σ error bars, a yellow circled dot (⊙), and a black plus sign (+), respectively. The dashed gray circles centered around the LSR represent boundaries for thin (radius of 85 km s−1) and thick (radius of 180 km s−1) disks following Fuhrmann (2004). A reference sample of 75 binaries that host extremely low-mass He WDs is shown as blue data points, with light blue marking objects that orbit in the Galactic disk (Brown et al. 2020).

In the text

Fig. 7. Preliminary and revised position of J1604+1000 (light red and dark red 1σ error bars) in the Kiel diagram. The black lines are evolutionary tracks for stripped helium cores from Driebe et al. (1998) labeled with their respective masses (stars evolve from top to bottom). For reference, the locus of the zero-age horizontal branch (ZAHB) for [Fe/H]= − 1.48 (Dorman et al. 1993) is shown as a dashed gray line. The meaning of the (light) blue data points is the same as in Fig. 6. The six proto-He WDs that are mentioned in the text – PSR J1816+4510, EVR-CB-001, GALEX 1717+6757, HD 188112, EL CVn, and WASP 0247−25 B – are presented in black.

In the text

Fig. A.1. Exemplary comparison of the best fitting model spectrum (red line) with normalized observation (black line; X-shooter spectrum taken on March 23, 2019, i.e., 2 458 565.8693 HJD). Light colors mark regions that have been excluded from fitting, e.g., due to the presence of interstellar or telluric lines. Residuals, χ, are shown as well. The optical spectrum only exhibits lines of H, He, Mg, Ca, and Fe.

In the text

	Fig. A.1. continued.
In the text

Fig. A.2. Corner plot visualizing the correlations between the orbital inclination (io), the radius and mass of the visible component (R1, M1), and the mass of its unseen companion (M2) for an RL filling factor x = 1 (see Eq. (11)). The diagonal panels show the respective posterior PDFs. Maximum values are marked by dashed red lines, and lower and upper bounds of the highest density interval are shown as dotted red lines. Prior PDFs for io and R1 according to Eqs. (4) and (5) are drawn as solid gray lines.

In the text

	Fig. A.2. continued for x = 0.9.
In the text

	Fig. A.2. continued for x = 0.8.
In the text

	Fig. A.2. continued for x = 0.7.
In the text

	Fig. A.2. continued for x = 0.6.
In the text

	Fig. A.2. continued for x = 0.5.
In the text