© The Authors 2024

1. Introduction

Classical Cepheids (henceforth: Cepheids) are evolved, luminous, intermediate-mass (∼3–11 M_⊙), radially pulsating stars. For more than a century, Cepheids have been crucial calibrators of the extragalactic distance scale (Hertzsprung 1913) thanks to the famous period-luminosity relation (PLR or Leavitt law, Leavitt & Pickering 1912). The unceasing interest in Cepheids is not only due to their cosmological relevance, but also because these are crucial populations to understand stellar physics and they are important tracers in galactic studies (Luck & Lambert 2011; Genovali et al. 2014; Poggio et al. 2021; Kovtyukh et al. 2022, and references therein).

Given the importance of Cepheids as distance tracers, many astrophysical effects have been explored as potential biases of the Leavitt law, including multiplicity (e.g., Mochejska et al. 2000; Kiss & Bedding 2005). For certain individual Cepheids, light contributed by companions can indeed affect distance estimates; the most extreme cases of this include V1334 Cyg in the Milky Way (Gallenne et al. 2018) and a sample of Cepheids with double-lined spectroscopic companions in the Large Magellanic Cloud (which are indeed selected using their Leavitt law outlier nature Pilecki et al. 2021). However, Anderson & Riess (2018) showed that statistically speaking, Cepheid companions cannot bias the extragalactic distance scale or the Hubble constant measured therefrom. A similar result was recently obtained using population synthesis (Karczmarek et al. 2022). However, orbital motion on timescales of up to several years may affect the accuracy of parallax measurements of Cepheids (Anderson et al. 2015, 2016a; Benedict et al. 2022).

Parallaxes of Cepheids provided by the third data release (DR3) of the ESA mission Gaia do not yet account for orbital motion due to companions. However, this is expected to be available in future data releases, and DR3 already provided a large number of non-variable stars with non-single star solutions (Gaia Collaboration 2023). Combining astrometric and radial velocity data to determine orbits is particularly useful because it allows the total mass of a system to be measured (e.g., Gallenne et al. 2019) and thus provides key information for constraining stellar models in the context of the mass discrepancy problem (e.g., Caputo et al. 2005; Anderson et al. 2016b).

In the past few decades, numerous complementary techniques have been used for the detection of Cepheid binarity and multiplicity. Direct spectroscopic identification of low-mass companions at any separation is possible using X-ray and ultraviolet (UV) observations (Evans 1992; Evans et al. 2013, 2022). The detection of orbital motion through radial velocity (RV) variations is effective for a wide range of companion masses and spectral types; however, it necessitates long temporal baselines (Evans et al. 2015; Anderson et al. 2016a). Conversely, photometric binary detection methods, such as those using photometric amplitudes (Coulson & Caldwell 1985; Klagyivik & Szabados 2009) or photometric colors and ratios of amplitudes of different colors (Gieren 1985), require confirmation of the binary detection through subsequent spectroscopic investigations. Cepheid companions can also be detected using optical long-baseline interferometry, if the contrast and angular separation between the Cepheid and its companion are sufficiently small (Gallenne et al. 2014a, 2015, 2016, 2019).

Gaia proper motions have been introduced as indicators of multiplicity (Kervella et al. 2019a,b, 2022). While proper motion analysis is a powerful tool to detect resolved companions, it cannot be applied to the brightest stars that saturate the Gaia detectors. Furthermore, Gaia astrometric quality flags such as the reduced χ² statistic are used as indicators of binarity assuming that the orbital motion leads to a photocenter wobble (Belokurov et al. 2020; Stassun & Torres 2021). In fact, several cuts on Gaia quality flags such as the RUWE and astrometric excess noise have been defined to filter the potentially multiple systems (Gaia Collaboration 2023). However, these cuts are established using stars where astrometric processing is simpler than in high-amplitude chromatically variable stars, such as classical Cepheids. It remains unclear if the same cuts can be applied to Cepheids, which are brighter than most stars used for calibrating Gaia systematics (cf. Khan et al. 2023 and references therein), and whose variability can cause excess noise in astrometric measurements.

The traditional approach of RV variations provides an unambiguous detection of a binary companion and helps the orbital parameters to be constrained. Dedicated RV monitoring campaigns have unraveled the orbital elements’ distribution of Cepheids (e.g., Evans et al. 2015; Szabados 2003). Low-mass companions such as the one of δ Cep (Anderson et al. 2015) or companions of fainter Cepheids (m_V > 8 mag) require significant observational effort and precision to be detected, and this can only be achieved by long-term RV monitoring.

The VELOcities of CEpheids (VELOCE) project provides unprecedented RV time series of classical Cepheids (Anderson et al. 2024; Paper I henceforth). Taking advantage of the 11-year baseline of VELOCE, here we report the results of a systematic study of single-lined spectroscopic binaries (SB1s). Where the data allow it, we provide estimates of orbital periods, or describe the orbital motion using different methods. The paper is structured as follows. The observations and data analysis are presented in Section 2, in Section 3 we describe the different methods we used to analyze the SB1 Cepheids, followed by results in Section 4, which is divided into subsections for the results from each method used to describe the binary. In Section 5, we present the correlation between derived orbital parameters and compare our results with the Gaia binarity indicators. We briefly discuss and summarize this work in Section 6. The RV data as well as any fitted models, including orbital solutions and polynomial trends, used here are made publicly available as part of Paper I. In the appendix here (Appendix A) and in the on-line supplementary m aterial database¹, we have presented supporting figures and tables. Lastly, in Appendix B, we provide detailed additional notes on individual systems.

2. Sample selection and available data

VELOCE comprises of precise Cepheid RV time-series data collected from both hemispheres using the high-resolution spectrographs Coralie (Queloz et al. 2001) and Hermes (Raskin et al. 2011). We used a minimum of 30 RV measurements per star, extending to up to 300 RV measurements for some cases. These observations cover a temporal baseline of nearly eleven years. Most of the VELOCE observations targeted signal-to-noise ratios (S/Ns) of at least 20, more commonly between 25 to 30, per pixel near 5500 Å. With such S/N we have typical RV uncertainties in the order of 50 ms⁻¹. A detailed description of the observations and the RV data can be found in Paper I.

We identify SB1 candidates from 258 classical Cepheids observed by the VELOCE project (cf. Paper I for details) via time-variable pulsation-averaged velocities, v_γ. Orbital motion in Cepheids introduces Doppler shifts on timescales of typically one year or longer that are mostly well separated from the pulsational timescales (of order days to months). Using VELOCE data, orbital motion can be detected in several different ways, depending on the availability and quality of the measurements (Section 3). The cleanest examples in our sample are Cepheids where clear, and possibly repeating, variations of v_γ are discernible by eye in the pulsation (Fourier series) fit residuals. Such cases generally have orbital semi-amplitudes of several km s⁻¹ and orbital periods of ≲3 yr.

While the selection function of Cepheids in VELOCE is rather complex, cf. Paper I, we consider VELOCE neither particularly biased toward the detection of SB1 systems, nor against it. For example, VELOCE initially specifically targeted Cepheids with no prior RV information and Cepheids residing near star clusters (Anderson et al. 2013). Later, Cepheids bright and nearby enough to measure interferometric radii were added (e.g. Anderson et al. 2016c; Breitfelder et al. 2016), and long-period Cepheids were targeted specifically for their use in distance ladder calibration (Anderson et al. 2016a). VELOCE also followed Cepheids with the goal of detecting modulated variability (e.g. Anderson 2014, 2016, 2019, 2020). Of course, stars found to exhibit orbital motion were followed up with the goal of characterizing orbits. However, only few stars were added to specifically search for suspected orbital signals.

Table 1 contains the list of stars considered in the present paper, ordered by whether their nature as spectroscopic binary was first reported here, was previously discussed in the literature. In Figure 1, we present the on-sky distribution of the sample SB1 Cepheids and the range in Gaia G mag that they span.

Fig. 1. Location of the sample SB1 Cepheids on the sky in the equatorial (top panel) and galactic (bottom panel) coordinate system. The symbols are color coded in the Gaia G magnitude. The symbols surrounded with blue circles are systems for which we provide an orbital solution in the current work.

3. Data treatment and methods

3.1. RV curve models adopted

In most Cepheids, the RV variability due to pulsation dominates the RV time series with pulsational peak-to-peak amplitudes typically ranging from 20 − 70 km s⁻¹. It is customary to model the pulsational variability using a Fourier series with a number of harmonics adapted to match the complexity of the RV variability. Paper I describes this process for all Cepheid RV curves in VELOCE. In short, the RV signal due to pulsation is modeled as:

$$\begin{array}{c}\hfill v_r(t)=v_{\gamma }+{\displaystyle \sum _n}{a}_{n}cos(2n\pi \varphi )+b_{n}sin(2n\pi \varphi ),\end{array}$$(1)

where n represents the order of the harmonic, t is time (in heliocentric Julian Date) and ϕ = (t − E)/P_puls is the pulsation phase with epoch²E and the pulsation period, P_puls. As described in Paper I, the optimal model is determined by fitting each RV time series using up to 20 Fourier harmonics and then performing model comparisons using the F-test and the Bayesian Information Criterion (BIC).

The key evidence of the spectroscopic binary nature presented here is a long-term modulation of v_γ = v_γ(t), typically on timescales > 1 yr and much longer than the pulsation period. We thus identify binaries by the correlated fit residuals they produce (Eq. 1). Orbital motion is then either modeled using a Keplerian orbit or approximated via linear, quadratic, or even higher order polynomials if the orbit is not sufficiently sampled to determine the full orbital solution.

As noted by Anderson (2014), Cepheids can exhibit modulated RV curves that can be misinterpreted as showing evidence of orbital motion. We note that all targets discussed here are bona fide spectroscopic binaries as identified by the absence of modulated line shape variability (Anderson 2016), although a few of the binaries presented here also show additionally modulated spectral variability. A detailed description of modulated RV variability of VELOCE Cepheids will be presented in future work.

3.2. Keplerian fitting

In the case of spectroscopic binaries, v_γ = v_γ(t) in Eq. (1) varies with time. The degree to which such variation can be used to infer the orbital elements depends largely on the sampling, baseline, and precision. Depending on these factors, we employed two methods for fitting Keplerian orbits to the data.

When fitting pulsational variability and orbital motion simultaneously, we fitted the Fourier series (Eq. (1)) and a Keplerian model simultaneously using the nonlinear least squares routine leastsq available in scipy.optimize (Virtanen et al. 2020). Following, for example, Hilditch (2001), the orbital motion is represented as:

$$\begin{array}{c}\hfill v_{\gamma }(t)=v_{\gamma ,0}+K[cos(\omega +\theta )+ecos\omega ],\end{array}$$(2)

where ω denotes argument of periastron, K semi-amplitude (in km s⁻¹), e eccentricity, θ true anomaly, and v_γ, 0 the systemic RV relative to the solar system barycenter. We note that v_γ (without subscript 0) refers to the pulsation averaged velocity (cf. Paper I), which coincides with v_γ, 0 if an orbital solution is determined or if the star does not have a companion. The uncertainties on the parameters are determined using the covariance matrix returned by the fitting routine, which incorporates statistical information on the reliability of the parameters estimated. Two examples of combined Fourier series plus Keplerian orbits for classical Cepheids are illustrated in Figures 2 and 3.

Fig. 2. Pulsation and orbital fit of the naked-eye Cepheid prototype δ Cephei. Left top: pulsation RV variability where the systemic velocity, vγ, 0 is indicated by a dashed line. Right top: residuals after fitting the Fourier series and Keplerian orbit against observation date. Left bottom: phase-folded orbital RV variation against orbital phase. Right bottom: orbital RV variation versus the observation date.

Fig. 3. Pulsation and orbital fit of one of the faintest binaries in our sample, NT Puppis. The figure description is the same as Figure 2.

To understand and mitigate any dependence of our results on starting values, as well as to avoid convergence to a local minimum, we also implemented a Markov chain Monte Carlo (MCMC) method for determining Keplerian orbits. The MCMC fits further allow one to explore the posterior distribution of the orbital parameters to identify possible correlations and to rigorously quantify uncertainties (e.g. Ford 2005, 2006). However, we have not yet implemented an MCMC fit to simultaneously model pulsational and orbital variations; this is foreseen for future work. Thus, all MCMC orbital solutions presented here use template fitting residuals of VELOCE (cf. Section 5.1 of Paper I for more details) as well as literature radial velocity data, where pulsational variability has been removed (cf. Section 3.4). The MCMC orbits presented in this work utilize an extended baseline, achieved through the combination of VELOCE with literature datasets (V+L).

After testing different orbital parameterizations, we found that the following implementation converged the fastest (cf. also Fulton et al. 2018):

$$\begin{array}{c}\hfill \mathit{x}=(v_{\gamma ,0},ln{P}_{\mathrm{orb}},lnK,e,\omega ,{\varphi }_0),\end{array}$$(3)

where ϕ₀ is the orbital phase to the next pericenter passage at t = E. We opted for a uniform prior for all parameters, and set sharp boundaries far away from possible solutions so that they would not alter the posterior distribution. Initial guesses were obtained by computing the Maximum A Posteriori (MAP)³ probability estimates from a grid of starting values of the parameters. The resulting orbits allowed us to visualize the possible different sets of parameters that could explain our data and have a sense of the stability of each solution by counting how many grid points converged to the same orbital parameters. We then chose the best parameters as the initial guess for our MCMC (an example of MCMC orbital fitting is provided in Figure 4).

Fig. 4. Orbit of U Vul determined by fitting vγ epoch residuals of literature and VELOCE data.

We used the MCMC implementation provided by the python package emcee (Foreman-Mackey et al. 2013) and, as suggested by the authors, used the autocorrelation time, τ_ACT, to determine whether the chain was sufficiently long. In particular, the convergence criterion consisted of the following inequalities:

$$\begin{array}{cc}& ⟨{\tau }_{\mathrm{ACT}}⟩<N_{\mathrm{iteration}}/60,\hfill \end{array}$$(4)

$$\begin{array}{cc}& \frac{\mathrm{\Delta }{\tau }_{\mathrm{ACT}}}{{\tau }_{\mathrm{ACT}}}<0.01.\hfill \end{array}$$(5)

After running the MCMC, we visually inspected the resulting posterior distributions to check the absence of significantly non-Gaussian distributions (Figure 5). We also verified that the mean and standard deviation were compatible with the 16th, 50th, and 84th percentiles. We use the mean and standard deviation as the final results for the orbital estimates and its uncertainties. Finally, we also investigated circular orbits in a few cases where e was either very small or not significantly detected. In this case:

$$\begin{array}{c}\hfill \mathit{x}=(v_{\gamma ,0},ln{P}_{\mathrm{orb}},lnK,{\varphi }_0),\end{array}$$(6)

Fig. 5. Joint posterior distributions derived from the MCMC sampling for U Vul. The histograms show marginalized posterior distributions of the parameters.

whereby ϕ₀ represents the orbital phase to the next outgoing node at t = E.

Using orbital elements, we calculate the projected semimajor axis as:

$$\begin{array}{c}\hfill asini[\mathrm{au}]=\frac{1}{2\pi G}K{P}_{\mathrm{orb}}\sqrt{1-e^2},\end{array}$$(7)

and the mass function as

$$\begin{array}{c}\hfill f_{\mathrm{mass}}[{\mathrm{M}}_{\odot }]=\frac{1}{2\pi G}{K}^3{P}_{\mathrm{orb}}{(1-e^2)}^{\frac{3}{2}},\end{array}$$(8)

where G denotes the gravitational constant and the mass of the Sun was assumed to be M_⊙ = 1.98840987 × 10³⁰ kg (as used in astropy.constants).

3.3. Orbital trends represented using polynomials

Cepheids whose pulsation residuals computed using Eq. (1) indicated the presence of incompletely sampled orbital motion were modeled using the sum of Eq. (1) and a low-order polynomial of the form:

$$\begin{array}{c}\hfill v_{\mathrm{trend}}(t)={\displaystyle \sum _{k=1,\dots ,i}}{p}_k·{(t-E)}^k.\end{array}$$(9)

The polynomial degree i adopted depends on the structure observed in the pulsation residuals (e.g., linear or nonlinear) and the overall v_γ(t) variation over the VELOCE temporal baseline. We prioritized low values of i, typically < 3, and sought to obtain a reasonably close fit to the data. Additional data are required to determine full orbital solutions in these cases. It is important to note that in the present study, Fourier+polynomial fitting was exclusively employed for incompletely sampled orbits within VELOCE. This differs from the approach in Paper 1, where polynomials were utilized to represent modulated variability of any origin.

We note that v_trend(t) represents the change of the pulsation averaged velocity v_γ(t) over time, so that v_γ in Eq. (1) becomes the pulsation-averaged velocity near the midpoint of the RV time-series. Introducing such polynomials allows one to faithfully recover the pulsational variability and to partially quantify the orbital motion already sampled by the available data. However, it should be noted that in this case v_γ is specific to the time-series and not generally valid for the star. Specifically, v_γ then does not refer to the center of mass line-of-sight velocity of the system with respect to the Sun.

3.4. RV template fitting (RVTF)

We investigated variations of v_γ over timescales exceeding the VELOCE baseline using RV datasets from the literature. To this end, we adopted the literature RV zero-point differences from Paper I and the RV template fitting (RVTF) approach developed in Anderson et al. (2016a), Anderson (2019). In this approach, the pulsational RV curve is assumed to be fully characterized by the high-quality pulsation fits obtained in Paper I using VELOCE data of each star. The template fitting technique then fits the VELOCE pulsation models to the available literature data while solving for two parameters, an offset in v_γ(t) with respect to VELOCE, and a phase shift Δϕ. While the former is required to investigate binarity, the latter allows one to efficiently deal with pulsation periods changing over the course of several decades without having to assume a functional form for such period changes (cf. Section 5.2 of Paper 1).

In Paper I, we developed a semi-automated approach for clustering the available data sets using a Kernel Density Estimation algorithm whose bandwidth was varied according to data availability and quality. For binaries, we double checked the most appropriate bandwidth to use to avoid fitting data exhibiting significant variations of v_γ as one epoch. Furthermore, we visually inspected all template fits to make sure they covered a sufficiently broad range in phase to be informative. Template fits restricted to an insufficiently narrow range in phase were discarded. For each clustered template fit and per literature reference, we obtained the mean date, range of dates, fit parameters Δv_γ, Δϕ, their uncertainties, as well as the time-series epoch residuals (denoted as ϵv_γ henceforth).

To investigate the significance of orbital signatures, we calculated the S/N between velocity offset and square-summed uncertainties as follows:

$$\begin{array}{c}\hfill \mathrm{S}/{\mathrm{N}}_{\mathrm{SB}1}=\frac{|\mathrm{\Delta }{v}_{\gamma }|}{\sqrt{{({\sigma }_v)}^2+{({\sigma }_{v_{\gamma }})}^2}},\end{array}$$(10)

where Δv_γ is the difference between the zero-point-corrected⁴v_γ obtained by template fitting and the pulsation averaged velocity from Paper I. σ_v denotes the uncertainty on v_γ, and σ_vγ is the uncertainty on the instrumental zero-point correction v_γ. Using Eq. (10), we searched for stars exhibiting deviations of more than S/N_SB1 > 3 and rejected any spurious v_γ variations among closely neighboring cluster sets following visual inspection, cf. Section 4.3.

Where feasible, we determined full orbital solutions using our MCMC algorithm (cf. Section 3.2) using the ϵv_γ from template fitting as input. In practice, this often worked better than trying to fit a combined Fourier series and Keplerian model, since period changes could be effectively removed using the template fitting algorithm. The v_γ and ϵv_γ of literature datasets used in the orbital fits with MCMC algorithm and in the RVTF analysis are included in the VELOCE data files published as part of Paper I.

Please note that a thorough analysis comparing Gaia DR3 RVs with VELOCE RVs was performed in Paper I. The main arguments against including Gaia data in the current work were as follows: Firstly, their temporal baselines mostly overlapped with existing data, and no significant improvements in orbital coverage were identified. Secondly, we observed that certain stars displayed temporal trends in Gaia radial velocities (RVs) over timescales where VELOCE RVs showed no variations in pulsation-averaged velocity. However, a clear pattern regarding when and for which stars this occurred could not be discerned. In light of these disparities, Gaia RVs were omitted as a literature source for the spectroscopic binary search employing the method described above. With extended baselines, increased observation frequency, and improved data processing, future data releases with Gaia DR4 and DR5 are poised to excel in detecting SB1 systems and refining orbital solutions.

4. Results

This section describes the results of our systematic search for spectroscopic binary Cepheids using VELOCE data. An overview of the results is presented in Table 2. Tables 3 and 4 list orbital parameters for 18 Cepheids whose orbital parameters are being reported for the first time and 15 Cepheids whose orbits had previously been reported in the literature, respectively. The full versions of Tables 3 and 4 and figures illustrating the orbital solutions can be found with the supplementary m aterial available at the VELOCE zenodo repository⁵. Overall, we identified 76 bona fide SB1 Cepheids and 14 additional SB1 candidates in the VELOCE sample. 32 of these are reported as SB1 systems for the first time. Unless otherwise specified, we here discuss only bona fide SB1 Cepheids.

We first discuss fully determined orbits in Section 4.1. Sections 4.2 and 4.3 present Cepheids with incompletely sampled orbits represented by polynomials and investigated using template fitting, respectively. We discuss the spectroscopic binary fraction in Section 4.4 and provide brief notes on individual Cepheids in Appendix B.

4.1. Full orbital solutions

This section presents the results obtained by fitting full orbital solutions to the RV time-series. We determined complete orbital solutions for 33 stars (including three tentative orbital fits). Section 4.1.1 begins with results based on VELOCE data alone, followed by results based on combined VELOCE and (zero-point offset-corrected) literature RVs in Section 4.1.2. Tables 3 and 4 distinguish these solution types using the tags “V” and “V+L”, respectively.

4.1.1. Orbits from VELOCE data alone (V)

Thanks to the high RV precision and data homogeneity, orbits derived using VELOCE data alone are typically very accurate. However, the available 11 yr baseline is not always sufficient to sample orbits fully, and occasional gaps in the time series may unfortunately coincide with particularly important orbital phases, preventing a definitive determination of all orbital parameters. All systems presented here were fitted using the combined Fourier series plus Keplerian model, which simultaneously solves for pulsational and orbital variability. The resulting orbital periods range from 237 to approximately 4000 days, that is, the baseline available from VELOCE DR1.

Figures 2 and 3 illustrate two examples of combined Fourier series plus Keplerian orbits for the naked-eye prototype of classical Cepheids, δ Cephei, and the faint southern target NT Puppis, cf. Appendix B for additional details. All the other combined fits can be found online⁶.

The choice of starting values for model fit parameters, notably P_orb and e, can impact the determination of orbital solution. We therefore adopted reasonable starting values estimated by inspecting the Fourier series fit residuals. Where starting values could influence the result, we cross-checked the results using the MCMC template fitting algorithm, cf. Section 3.2. For most stars, however, the data set sampled the orbital cycles sufficiently well to allow an initial guess close to the correct orbital period.

The best fit orbital parameters obtained here generally agree with previous orbital solutions published in the literature despite improved uncertainties. However, we obtain significantly different orbital periods for six Cepheids δ Cep (cf. Figure 2), AX Cir, MU Cep, V1334 Cyg, W Sgr, and XX Cen. These stars, among others, are discussed in detail in Appendix B.

4.1.2. Orbits from VELOCE and literature data (V+L)

We determined orbital solutions using combined VELOCE and literature data using our MCMC algorithm applied to the template fitting epoch residuals; such fits are labeled as V+L in Tables 3 and 4. Orbits of five newly discovered binaries (FO Car, MY Pup, R Cru, VY Per, and VZ Pup) and nine previously known binaries (DL Cas, FF Aql, S Mus, S Sge, SY Nor, TX Mon, U Vul, V1334 Cyg, and YZ Car) were obtained using VELOCE data alone as well as in combination with literature RVs. The V and V+L orbits for 12 of these 14 agree within the uncertainties, although V+L orbits tend to have lower uncertainties, mostly due to better orbital sampling. The exceptions are MY Pup and VZ Pup. The VELOCE baseline only marginally covers the orbital period of MY Pup, hence, incorporating literature data results in a modification of the period. The P_orb of VZ Pup differs by 2.4σ between V and V+L orbits due to the P_orb being on the order of the available VELOCE baseline. The V+L result is thus clearly preferred, hence, whenever available we use V+L orbital estimates in any sample descriptions presented here. Lastly, for two stars, AX Cir and SU Cyg, we could only obtain a combined fit using VELOCE and literature data. In the case of AX Cir, due to its very long orbital period and for SU Cyg due to the incomplete sampling of orbital phase with VELOCE data alone.

4.1.3. “Tentative” orbits using VELOCE and literature data

We determined first tentative orbital solutions for three stars: BP Cir, R Mus, and V0659 Cen. We consider these orbits tentative because the literature data were not sufficiently constraining, while the VELOCE data alone were not sufficient to sample a full orbit. Additional observations are required to ascertain the true orbital solutions for these stars. Please note that the uncertainties associated with the orbital parameters of the tentative orbits listed in Table 3 may not accurately reflect the limitations arising from the insufficient sampling of these long orbits. The stated uncertainties are formal uncertainties on the model given the data. However, in cases of poor or incomplete orbital sampling, it is possible that the true orbital period is outside the range of these formal errors. This is why we present these cases separately as tentative.

4.2. Incompletely sampled orbits modeled as polynomial trends

We represented 16 newly detected SB1 systems and 14 previously reported SB1 systems using a combined Fourier series plus polynomial trend. Table 5 quantifies the measured range of the orbital RV variation, A_p2p, over the temporal baseline (denoted by ΔT_p2p), which is the baseline of the observations over which the A_p2p was calculated. Among the 30 systems examined with Fourier plus polynomial fit, 15 Cepheids exhibit nonlinear trends. Among the 15 stars with linear trends, 5 demonstrate a decreasing trend in v_γ, whereas the remaining 10 feature increasing v_γ. Additional observations are required to determine these orbits.

Figure 6 presents the distribution of A_p2p versus ΔT_p2p and the pulsation period. A clear trend of increasing A_p2p as we extend to longer temporal baselines is visible. For a fixed value of K, this trend arises because longer baselines increasingly sample the full (albeit incomplete) peak-to-peak variation of the orbit. However, there is significant dispersion in the value of K at fixed P_orb due to different inclinations and mass ratio. The smallest trends of between 0.2–0.35 km s⁻¹ are found for V1162 Aql, DR Vel, FR Car and SX Vel (over ∼1.8, 4.5, 5.8 and 5.9 yr respectively). By comparison, the smallest K among completely solved orbits is 0.67 km s⁻¹ (MY Pup, 5.7 yr). At the long timescale limit, we find A_p2p = 1.2 km s⁻¹ variation for T Mon over the course of nearly 10.5 yr. By comparison, the lowest K = 2.1 km s⁻¹ for a completed orbit and a similar P_orb of 11 yr is found for VZ Pup.

Fig. 6. Distribution of Ap2p, maximum amplitude variation in vγ, for the stars in our sample where we could fit a polynomial to the vγ trend (see Section 3.3 for details). In the left panel, we have the Ap2p vs the baseline in time over which it was calculated and in the right panel, the same against the pulsation period. The newly discovered SB1s from our sample are plotted in filled squares, while the literature-known SB1s are plotted in filled circles. The color bar represents the degree of polynomial fitted in the Fourier+Polynomial model.

4.3. Template fitting for stars without orbital signatures in VELOCE data

Paper I presented a list of SB1 Cepheids exhibiting time-variable average velocities detected using VELOCE data alone. SB1 Cepheids with very long P_orb could, however, remain undetected based on VELOCE data alone. Hence, we investigated the temporal stability of v_γ for all Cepheids in VELOCE using the well-defined RV templates from Section 3.4. To this end, we computed S/N_SB1 (Eq. 10) and considered a star a likely binary if any template fit resulted in a 3-σ deviation (S/N_SB1 > 3) from a constant value of v_γ.

Applying this method to 174 bona fide VELOCE Cepheids with available literature data resulted in the discovery of four additional SB1 candidates: β Dor, SZ Aql, V0402 Cyg, and X Cyg, cf. Figure 7. All four exhibit visually convincing trends in v_γ on the order of 1 − 3 km s⁻¹ on timescales exceeding 15 000 − 20 000 d.

Fig. 7. Zero-point corrected Δvγ for new SB1s detected through RVTF. The dotted lines mark the standard deviation of the vγVELOCE. This figure only includes stars where we have a 3-σ detection through our analysis. The color bar represents the S/NSB1 of the various vγ measurements where the S/NSB1 was calculated as described in the text.

Additionally, we found 8 bona fide SB1 Cepheids (CD Cyg, η Aql, RS Ori, SS CMa, SZ Cyg, V0340 Ara, V0916 Aql and VY Sgr) among a sample of 32 SB1 candidates previously reported in the literature (Szabados 2003) where the available baseline of VELOCE data alone was insufficient to detect orbital motion. Figure 8 illustrates these cases. CD Cyg and SS CMa were already discussed as potential SB1s from the RVTF analysis of Anderson et al. (2016a) using VELOCE data. Here, we confirm that the deviations in v_γ indicate of binarity.

Fig. 8. Literature-known SB1 candidates detected through RVTF. The figure description is the same as Figure 7.

For the remaining 24 Cepheids among the 32 candidates, we could not confirm unambiguous orbital motion despite large overlaps in the data used. Ten of these 24 Cepheids do not exhibit a S/N_SB1 > 3, over a temporal baseline of 10 000–20 000 days. For the remaining 14 cases, individual v_γ variations are found, although they should be considered spurious variations because closely neighboring v_γ estimates clearly show the absence of a trend at these epochs. In the supplementary m aterial, we present all these non-SB1 cases and we have tagged these stars with RVTF-No in the Binarity column of Table 1.

Given the high accuracy of the VELOCE template fitting procedure, which effectively removes any pulsational variability and period variations, and since instrumental zero-point differences have been corrected using nonbinary Cepheids (cf. Paper I), we conclude that these 24 Cepheids should not be considered spectroscopic binaries.

4.4. SB1 fraction among VELOCE Cepheids

The 76 bona-fide SB1 Cepheids detected as part of VELOCE yield a lower limit to the spectroscopic binary fraction of 29.6 ± 3.4% that is fully consistent with the result by Evans et al. (2015) of 29 ± 8% based on the 40 brightest Cepheids. If all 14 SB1 candidates from Paper I should be confirmed to be bona-fide SB1 systems, then the VELOCE binary fraction would increase to ∼35%.

Using the subset of Cepheids with fully determined orbital solutions from Tables 3 and 4 and assuming that all undetermined orbits have very long orbital periods, we find a lower limit to the fraction of SB1 Cepheids with orbital periods of less than 10 yr of 15.2 ± 2.4%, slightly lower albeit consistent with the reported 20 ± 6% from Evans et al. (2015). Although our stringent template-based investigation (Section 4.3) questioned the veracity of the SB1 nature of 24 Cepheids for which such evidence had been previously discussed (cf. references in Szabados 2003), the fraction of SB1 systems among Cepheids we determined remains similar because we also identified new SB1 systems. However, the more accurate determination of time-variable v_γ alongside detailed zero-point corrections affords much greater confidence in the determination of the SB1 nature of individual Cepheids. Of course, the SB1 fraction of Cepheids is merely a lower limit of the total fraction of Cepheids with companions due to the possibility of very large separations, unfavorable inclinations, and extreme mass ratios.

Figure 9 illustrates the magnitude distribution of newly discovered and detected SB1 Cepheids alongside the distribution of binary systems reported in the online compilation by Szabados (2003), as well as the distribution of magnitudes for which orbits were determined here. To ensure consistency in the comparison, we have exclusively included the bona-fide SB1 binaries from Szabados (2003) in Figure 9. Therefore, other types of binaries (such as photometric, visual, etc.) from their catalog are not accounted for in this analysis. Please note that Figure 9 does not include 24 systems that are not detected as SB1 in our study.

Fig. 9. Stacked histogram presenting the new VELOCE SB1 detections (left panel) and the new VELOCE orbits (right panel) along with the same thing from Szabados (2003). SY Nor’s orbit is provided by us as well as Bersier (2002); however, it is missing in the Szabados (2003) database.

As the right panel of Figure 9 shows, we here double the number of orbital solutions for Cepheids fainter than ∼8-th magnitude, and there are only 14 Cepheids with known orbital solutions that are not included in this work (one Cepheid among these 14, namely T Mon has incompletely sampled orbit in this study). Conversely, we determine the first orbital solutions for 18 Cepheids, most of which fainter than ≳8 th magnitude.

5. Sample analysis

5.1. Correlations among orbital elements

Figure 10 shows the correlation between different orbital elements (ω, K, e, and a sin i) with P_orb for all orbital solutions reported here. This includes first-time orbital estimates of 18 stars and 15 VELOCE orbital estimates for Cepheids with literature-known orbits. Newly discovered SB1 systems populate the full range of orbital parameters, indicating that our sample selection was not biased towards the literature-known orbital properties of SB1 Cepheids. Semi-amplitudes tend to decrease as P_orb increases. Thanks to the precision of VELOCE RV measurements and long-term monitoring, we were able to discover many new SB1 systems with very small semi-amplitudes and long orbital periods. A noticeable correlation exists between a sin i and P_orb, and the dispersion along this correlation might stem, at least partially, from fluctuations in the combined masses of the Cepheid and its companion.

Fig. 10. Distribution of different orbital elements (Upper left: argument of periastron ω; upper right: eccentricity e; lower left: semi-amplitude K of the orbital RV variations; and lower right: projected semi-major axis a sin i) versus orbital period Porb for Cepheids with the orbits determined here. SB1s with newly reported orbits in this study are plotted using open squares, while VELOCE orbital estimates of SB1s with literature-known orbits are plotted in open circles. Symbols with double boundaries highlight V+L orbits determined using VELOCE and zero-point corrected literature RVs. R Cru, SY Nor and YZ Car are absent in the top left panel due to insignificant eccentricity.

The eccentricity versus orbital period distribution (upper right panel in Figure 10) is particularly interesting. We find that the maximum eccentricity increases with orbital period. Despite a large range of eccentricities at most P_orb values, an absence of very high e < 0.6 is clearly noticeable at P_orb ≲ 3 yr. The two Cepheids with the shortest P_orb have very low or no noticeable eccentricity. We also notice an absence of circular orbits at P_orb ≳ 3 yr, which appears to become an exclusion region at longer P_orb. Since such distant companions have most likely evolved in isolation, this could point to a primordial feature of intermediate-mass (B-star) binaries. However, NBODY simulations show that cluster dynamics increase e for long orbital periods over time (Dinnbier et al., in prep.). Additional long-period orbits are required to investigate this effect.

Most of these trends qualitatively agree with the empirical results by Evans et al. (2015), as well as with predictions from the Cepheid binary population synthesis study by Neilson et al. (2015) and Karczmarek et al. (2022), as well as with the dynamical NBODY6 simulations by Dinnbier et al. (2024). However, we do find several Cepheids at high 5 < P_orb < 10 yr and relatively low eccentricities (e < 0.2) where Cepheid binaries have been predicted to have a low probability of occurring (Neilson et al. 2015). Additionally, R Cru is very close to the minimum P_orb consistent with the synthetic populations by Neilson et al. (2015) and Karczmarek et al. (2022).

Incompletely sampled orbits described in Section 4.2 and illustrated in Figure 6 provide partial information on K and P_orb. Comparing A_p2p with K (in Figure 10) as a function of P_orb, we notice the inverse trends. Shorter orbital periods tend to have larger K, whereas we find larger A_p2p over longer baselines. This suggests that most short-P_orb orbits with K ≳ 0.5 km s⁻¹ of Cepheids observed by VELOCE have been determined. Very low A_p2p values are likely to continue growing over longer baselines.

5.2. Cepheids exhibiting proper motion anomaly

Kervella et al. (2022) used Gaia to detect astrometric binaries using the proper motion anomaly (PMa). This method compares the proper motion vector determined by the ESA HIPPARCOS mission (Perryman et al. 1997) and various Gaia data releases (Gaia Collaboration 2016a,b, 2018, 2021; Michalik et al. 2015; Lindegren et al. 2018, 2021). PMa quantifies the significance of a long-term rotation of the proper motion vector due to orbital motion. Stars with PMa signals exceeding 3σ were thus labeled as binaries by Kervella et al. (2022). Proper motion probes the on-sky projection of orbital motion, orthogonal to the line-of-sight variations measured using RVs. Therefore, the PMa and RV methods function as complementary tools for binary detection (cf. Kervella et al. 2019b for details on the sensitivity function of PMa). Consequently, it’s typical for some stars identified as SB1 to lack discernible deviations in their proper motions which would categorize them as PMa-True candidates. This could either be if the Cepheids are located at large distances or possess orbital periods too short for detection. On one hand, companions with orbital periods shorter than ∼3 years are unlikely to be identified via PMa due to the time window smearing inherent in Gaia DR3. On the other hand, extremely long-period companions, spanning hundreds of years, are similarly challenging to detect, representing a limitation shared by VELOCE as well. Importantly, a PMa-False flag does not raise concerns about the SB1 detection. Among the 76 binary Cepheids in our sample, 19 have been reported to exhibit PMa signals (with binary flag “BinH2EG3b” in Kervella et al. 2022 as 1)⁷. Evidence of astrometric orbital signals in form of the PMa has been reported for 10 of 33 SB1 (30%) Cepheids with orbital solutions presented here (including tentative ones) as well as for 8 of 30 (25%) SB1 systems with long-term v_γ trends exhibit PMa (Kervella et al. 2022). This also includes 2 of the newly discovered SB1 Cepheids, GX Car and V0659 Cen.

Two systems with reported PMa, RZ Vel and SV Per, deserve a specific mention. For these two stars, we found no evidence of time-variable v_γ. However, SV Per is known to have a very nearby (0.2″) companion that has been resolved using HST/WFC3 spatial scans (Riess et al. 2018) and which is most likely unresolved by Gaia. Time-variable contrast differences due to the Cepheid’s pulsation would affect the measurement of the photocenter for SV Per, which may lead to a spurious PMa signal.

Finally, we note that none of the Cepheids considered to be non-SB1s in VELOCE–stars without evidence of time-variable v_γ-have been reported to exhibit a significant PMa. Specifically, this is true for all stars without evidence of time-variable v_γ based on VELOCE data alone as well as our RVTF analysis (Sect. 3.4), including stars for which our results contradict previous claims of time-variable v_γ.

As this comparison shows, RVs from VELOCE and literature data standardized to the VELOCE zero-points provide the most complete evidence of orbital motion for Cepheids on timescales of up to a few decades. Nonetheless, we caution that signals on the order of 1 km s⁻¹ and lower can be introduced also by modulated variability (Anderson 2014). CCF shape parameters, such as the bisector inverse span and full width at half maximum can be informative to this end (Anderson 2016) and will be considered in future work.

5.3. Gaia DR3 astrometric quality flags of SB1 Cepheids

Gaia DR3 reports several astrometric quality flags and parameters to aid the interpretation of the reported astrometric measurements. The Renormalized Unit Weight Error (RUWE) has received particular attention and is often considered an indicator of stellar multiplicity. RUWE is calculated as the square root of the normalized χ-square of the astrometric fit to the along-scan observations (Lindegren 2018), and the threshold of RUWE < 1.4 is frequently used to indicate well-behaved astrometric solutions (Lindegren et al. 2018, 2021). However, several factors can result in elevated RUWE, which does not identify the origin of the excess residuals and includes all possible sources of error in the fit to the astrometric model. The latter include any signals affecting photocenter stability, notably involving marginally resolved visual binaries and unmodeled astrometric orbital signals.

In Cepheids, the high-amplitude chromatic variability of Cepheids could also contribute to RUWE. For example, Gaia may automatically select different gating schemes and window classes depending on the momentaneous brightness of Cepheids that vary by ∼1 mag in optical bands. Additionally, chromatic variability may lead to noise due to time-variable chromatic diffraction properties. Given these added difficulties, it is likely for RUWE to be elevated for Cepheids compared to non-variable stars. We therefore investigated to what degree excess RUWE and other astrometric quality indicators reported by Gaia may serve as an indicator of multiplicity for Cepheids.

Figure 11 shows the comparison of Gaia DR3 RUWE against the projected semimajor axis and P_orb. The figure shows a significant excess RUWE near a sin i values close to 1 au as 5 out of 8 stars with the highest RUWE (RUWE > 2) fall into this category (stars located within the dashed red lines in the left panel of Figure 11). Interestingly, a corresponding feature is not so clearly seen near P_orb = 1 yr, suggesting that the size of the orbital ellipse on the sky more directly affects the astrometric solution than the periodicity of the signal. The 6300 d binary with the highest RUWE of 7.8 is AX Cir, whose hot main sequence (likely B6V) companion has been resolved interferometrically (Gallenne et al. 2014b). A photometric contamination is the potential origin of the elevated Gaia RUWE.

Fig. 11. Distribution of Gaia DR3 RUWE with the semi-major axis (left) and orbital period (right) of our sample systems. The stars tagged as binaries from PMa (Kervella et al. 2022) are marked with crosses. The red solid line indicates a sin i of 1 au, and the dashed red lines indicates the same at 0.5 and 2 au.

Figure 12 shows RUWE in Gaia DR3 vs the amplitude of trends, A_p2p, determined using VELOCE data only in Sect. 4.2. Given that the majority of these orbital signals occur on timescales longer than the VELOCE baseline, it is not surprising that RUWE < 1.4 for the vast majority of these stars. The most significant outlier at RUWE ∼ 8 is RW Cam, which was shown to feature a significant UV excess indicative of a hot main sequence companion (Stepien 1968). The star with the second highest RUWE of 1.7 is T Mon, host to a B9.8V companion detected directly from UV spectroscopy (Evans & Lyons 1994) but remained undetected interferometrically in H-band (Gallenne et al. 2019). These stars collectively indicate that elevated Gaia RUWE values in Cepheids may be attributed to unresolved companions rather than orbital motion.

Fig. 12. Distribution of Gaia DR3 RUWE with the maximum amplitude of the (linear or nonlinear) vγ variation for SB1 systems with Trend. The stars tagged as binaries from PMa (Kervella et al. 2022) are marked with crosses.

Figure 13 shows the probability density function of RUWE for SB1 and non-SB1 Cepheids in VELOCE. Approximately 25% of VELOCE SB1s have a RUWE > 1.4. While there is a slightly higher probability for an SB1 Cepheid to have excess RUWE compared to non-SB1 Cepheids, the difference is not very pronounced. Qualitatively, among the newly identified SB1 systems (RVTF), Figure 13 suggests similar probabilities near RUWE ≈ 1.7−1.8. However, among these, β Dor’s extreme brightness (m_V = 3.8 mag) likely contributes substantially to the high RUWE of 4.5.

Fig. 13. Probability density function of Gaia DR3 RUWE for VELOCE SB1 and non-SB1 systems.

Finally, we also investigated the Gaia astrometric quality indicators astrometric_gof_al and astrometric_excess_noise as indicators of multiplicity of Cepheids. However, no significant trends appeared in this comparison aside from a mild increase in astrometric_excess_noise near a sin i ∼ 1 au.

In summary, we find excess RUWE mainly for Cepheid binaries with projected semimajor axis close to 1 au. The parameter astrometric_excess_noise is also increased in this range, though the correlation is less clear. In addition to specific orbital configurations, elevated RUWE can arise from photocenter variations in unresolved binaries and in very bright stars, for example. Additionally, the majority of binary Cepheids are not identified by RUWE > 1.4.

5.4. Amplitude ratios as indicators of companion stars

Photometric amplitude ratios can be used to identify Cepheid companion stars in case the photometric contrast is sufficiently low (e.g. Madore 1977; Madore & Fernie 1980; Evans & Udalski 1994). Since most Cepheid companions are hot early-type stars, shorter-wavelength data is generally more sensitive to companions than longer-wavelength photometry. A noticeable exception are double-lined spectroscopic binaries (SB2) in the Large Magellanic Cloud that consist of a Cepheid and a giant companion (Pilecki et al. 2021).

In the Milky Way, Klagyivik & Szabados (2009) considered low ratios between photometric and RV amplitudes as an indicator for the presence of companion stars. Thanks to Gaia’s multiband photometry, such a comparison is possible here using a large sample of stars. Figure 14 thus shows the ratio between Gaia photometric and VELOCE RV amplitudes in three bands (G_Bp, G, and G_Rp) as a function of the logarithmic pulsation period. The higher ratios at shorter wavelength are readily explained by the larger photometric amplitudes. SB1 Cepheids identified in VELOCE are shown as filled circles, and stars which are not SB1 systems are shown as open circles. As a consequence of the Cepheid period-color relation, the contrast between Cepheid and companion should increase toward longer periods, rendering photometric detection more difficult. Additionally, any signal in the amplitude ratio should be enhanced at shorter wavelengths. However, the comparison between SB1 Cepheids and all others observed in all three bands reveals no indication that bona fide SB1 Cepheids have reduced amplitudes compared to other Cepheids (Figure 14). The only potential signature is a slightly reduced range of amplitude ratios for SB1 Cepheids, which, however, clusters in the center of the distribution. The comparison between photometric amplitudes in G_Bp vs G_Rp (not shown) yields the same impression. Hence, we find that Gaia’s photometric amplitudes are not good indicators of Cepheid companions.

Fig. 14. Photometric amplitude ratios as indicators of companion stars. First three panels from left to right: photometric-to-RV amplitude ratios from Gaia DR3 and VELOCE against the logarithmic pulsation period, with wavelength of the photometric band increasing from left to right. Right panel: ratio of Gaia photometric amplitudes in GBp to GRp. Filled circles show SB1 Cepheids reported in VELOCE. Other single-mode Cepheids observed by VELOCE are shown as open circles. There is no clear difference between SB1 Cepheids and all others in any of the amplitudes considered, nor a dependence of such a difference on log Ppuls.

6. Discussion and summary

We have presented the largest homogeneous investigation of 76 SB1 Cepheids to date based on the unprecedented RV data set provided by the VELOCE project (Paper I). An additional 14 SB1 candidates have been identified in Paper I, but are not discussed here due to their nature as candidate binaries. VELOCE data allow us to extend the search for binary Cepheids to fainter magnitudes, adding 32 new SB1 systems, and 18 first orbital determinations of Cepheids, 11 of which are fainter than m_V > 8 mag (cf. Figure 9). Thanks to these improvements, we sharpened the uncertainty on the SB1 fraction of Cepheids by a factor of ∼2.3.

We present 30 definitive and three tentative orbital solutions, which range from the discovery of the shortest P_orb in the Milky Way for R Cru (238 d) to ∼17 yr (AX Cir), and also include a correction to the orbital solution of δ Cep. For 18 Cepheids, we estimate orbital solutions for the first time. VELOCE observations are particularly sensitive to P_orb ≲ 5 yr, and the available temporal baselines differ from star to star. The ensemble of orbital solutions suggests two peculiar features in the e − log(P) diagram, which are a possible exclusion zone of low-eccentricity systems (0 ≤ e < 0.2) at long orbital periods (P_orb ≳ 10 yr) as well as a dearth of stars with orbital period in the range of 2.5 − 3.5 yr. Further observations and a more detailed analysis of detection limits will be needed to confirm these features.

Based on our search for SB1 Cepheids and a detailed investigation of evidence of orbital signals using literature, we find a SB1 fraction of 29.6 ± 3.4%. This is a lower limit since 14 SB1 candidates reported in Paper I were not considered in this fraction. Furthermore, we estimate the fraction of SB1 Cepheids with P_orb < 10 yr of 15.2 ± 2.4%, half of which have P_orb of less than 2.5 years. Our binary fraction estimates agree with previous reports, albeit with 2 − 3× lower uncertainty. Of course, SB1 Cepheids are detectable mainly for orbital periods of < 100 yr, and other methods must be considered to estimate the total binary fraction of Cepheids. Considerable effort is currently underway to comprehensively determine the binary fraction of Cepheids, employing a diverse array of techniques such as astrometry, spectroscopy, interferometry, and observations across various wavelengths. For example, Evans et al. (2022) estimated a total binary fraction of 57 ± 12% based on SB1 Cepheids, and Cepheids with evidence of companions from X-ray and UV observations. It is worth noting that future Gaia data releases, particularly DR4 and DR5, will contribute significantly with an abundance of astrometric orbits, enriching our understanding of these systems even further.

Several of the SB1 Cepheids studied here are part of triple systems, including AW Per (Evans et al. 2000), FF Aql (Udalski & Evans 1993), S Sge (Evans et al. 1993), V1334 Cyg (Abt & Levy 1970; Evans 1994), V0659 Cen (Evans et al. 2013, 2022), and W Sgr (Evans et al. 2009). Here, we further identify FO Car, RV Sco, RY Sco, and UX Per as triple systems with an inner SB1 and outer visual companion. Comparing the binary fraction and incidence of triples among Cepheids with dynamical NBODY simulations that assume an initial binary fraction of 100% (Dinnbier et al. 2024) suggests that a significant fraction of Cepheid progenitors (B-type stars) must be formed as triples or even higher-order systems.

Last, but not least, we find that the Gaia DR3 parallax quality indicator RUWE is elevated for Cepheid binaries whose projected semimajor axis, a sin i, is close to 1 au. However, upon assessing RUWE for both SB1 and non-SB1 from VELOCE, we ascertain that Gaia RUWE is not a reliable indicator for identifying Cepheid spectroscopic binaries and can be noticeably affected by photometric systematics, such as marginally resolved binaries or high apparent brightness. Additionally, RUWE tends to be slightly higher for Cepheids than for non-variable stars.

An examination of PMa in relation to various orbital properties of our stars does not reveal any correlation between the orbital orientation and PMa detection. Furthermore, all VELOCE Cepheids with previously reported PMa detections are also identified to be SB1 systems.

Precision radial velocities provide the most sensitive diagnostic of multiplicity of Cepheids to date. Future Gaia data releases will allow us to determine systemic masses of Cepheids using the combined astrometric and spectroscopic orbital signals. In Gaia DR3, only one Cepheid was reported with a non-single-star astrometric solution. Suspiciously, the orbital period of RX Cam was given as 33 months, matching exactly the observational baseline of Gaia DR3. However, RX Cam was unfortunately not observed as part of the VELOCE project. The upcoming fourth Gaia data release in combination with VELOCE data will be a treasure trove for determining accurate Cepheid masses and to establish a well-sampled mass-luminosity relation required to elucidate the physics of Cepheids, notably with respect to mixing processes, such as convection and rotation, required to explain the mass discrepancy problem (e.g. Prada Moroni et al. 2012; Anderson 2014).

Data availability

Supplementary material is available at https://zenodo.org/records/12818503

Acknowledgments

RIA, SS, and GV acknowledge support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 947660). RIA, and SS further acknowledge support through a Swiss National Science Foundation Eccellenza Professorial Fellowship (award PCEFP2_194638). SS would further like to extend their acknowledgement to the Research Foundation-Flanders (grant number: 1239522N). 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. The lead authors would like to acknowledge the very useful compilation of information on binary Cepheids provided by Laszlo Szabados at Konkoly Observatory (Szabados 2003, https://cep.konkoly.hu/intro.html).

References

Appendix A: Guiding camera Image

When available, the Coralie and Hermes guiding camera images were utilized to examine for triple systems. In Figure A.1, we present a guiding camera image of FO Car obtained from Coralie.

Fig. A.1. Coralie guiding camera image of FO Car. The center of the Cepheid is dark because of the fiber placement during the integration. The object at the top left is an outer (visual) companion, which renders the SB1 Cepheid FO Car a triple system.

Appendix B: Notes on individual binaries

We have presented evidence of the SB1 nature of 76 Milky Way Cepheids on timescales spanning up to 40 years, including 32 stars whose SB1 nature was not previously known. The following provides additional background for specific stars of interest, notably if our results disagree with previous studies. We note that several Cepheids were first reported to be SB1 systems based on early VELOCE results, such as FN Vel and MU Cep (Anderson 2013), BG Vel and RZ Vel (Szabados et al. 2013b), δ Cep (Anderson et al. 2015), and AQ Car, CD Cyg, KN Cen, SS CMa, SZ Cyg, XZ Car, VY Car and X Pup (Anderson et al. 2016a). Anderson et al. (2016a) noticed tentative evidence of v_γ variations in VY Car, AQ Car, SZ Cyg and X Pup. However, they ultimately deemed these findings inconclusive due to imprecise data in the literature and the occasional substantial fluctuations in pulsational variability. In the current work we successfully detect these v_γ variations for SZ Cyg (when VELOCE data is combined with zero-point corrected literature datasets, see Figure 8) and AQ Car (cf. Table 5, Figure B.1). Nevertheless, in the case of VY Car and X Pup, the v_γ values, when adjusted for zero-point, do not reveal any consistent variation indicative of binarity (see Section 4 in the supplementary m aterial).

α UMi (Polaris) was one of the first Cepheids whose nature as an SB1 system was realized (Moore 1929), and much literature exists on this star. Polaris was extensively studied using pre-2019 VELOCE data by Anderson (2019), who determined a full orbital solution and studied the stability of RV amplitudes as well as the presence of additional periodic signals. Given the very long ∼30-year orbital period, we did not update the orbital solution here and instead refer to the aforementioned publication as well as Torres (2023) for orbital parameters.

δ Cephei is the archetype of classical Cepheids and was discovered to be an SB1 system using an earlier subset of VELOCE data. Anderson et al. (2015) estimated an orbit with $P_{\mathrm{orb}}={2201}_{-6.31}^{+5.73}$ d, $K=1.5_{-0.080}^{+0.239}$ km s⁻¹, and $e=0.{647}_{-0.021}^{+0.038}$ by combining new RV data with RV measurements from the literature. A main difficulty of this approach was the correction of zero-point differences for datasets with short (∼1 yr) baselines. Here, we report a correction to the orbital solution, with P_orb = 3450.9 ± 19.6 d, K = 3.01 ± 0.02 km s⁻¹, and e = 0.745 ± 0.006, based purely on VELOCE data and hence not subject to the uncertainties associated with literature zero-point offsets. This updated orbit of δ Cep was already presented at the 2022 RR Lyr-Cep conference in La Palma (Shetye 2022). It agrees reasonably well with the results recently presented by Nardetto et al. (2024) based on a combination of High Accuracy Radial velocity Planet Searcher for the Northern hemisphere (HARPS-N) data and existing literature. Their estimated orbital period of 3401.8 d and eccentricity of 0.71 differ by 2−σ and 1.7−σ from ours. These small differences may result from their use of a heterogeneous dataset, or the better coverage of the periastron passage by the HARPS-N data. As Figure 2 shows, the orbit is now completely sampled and contained within the VELOCE baseline (3620 d). Importantly, the extreme ends of the orbital velocity have been observed and tightly constrain e and K. Unfortunately, including literature data for a V+L orbit does not improve the orbital solution due to the same issues that affected the prior analysis. Extending the baseline with future Hermes observations will further improve the orbital period. We note that the mass function and the high brightness contrast imply the companion to be a low-mass star, possibly a white dwarf (Gallenne et al. 2016). A similar deduction was made in the recent X-ray study of Cepheids by Evans et al. (2022). A companion of G spectral type or later would be expected to be detected in X-rays, which puts a lower limit on the mass of the companion. Finally, we note that the single-star astrometric solution provided in Gaia DR3 clearly indicates poor fit quality⁸. It thus appears likely that Gaia’s next data release could deliver an astrometric orbital solution, as predicted by Anderson et al. (2015).

η Aql was recently studied by Benedict et al. (2022) using spectroscopic, photometric and astrometric (from the Hubble Space Telescope (HST) Fine Guidance Sensor) observations, including unpublished VELOCE data. No significant variations in the RV residuals were detected after removing the Cepheid pulsational RV signature, leading Benedict et al. (2022) to conclude that either η Aql has a very long orbital period (longer than their baseline of ∼38 years) or that it has a face-on orbit. We considered an extended temporal baseline by including RV data from Lloyd Evans (1980), Barnes et al. (1987), and Wilson et al. (1989) and detect a 3-σ variation in the Δv_γ (Figure 8). If the older datasets are to be believed despite large RV uncertainties of ∼1 km s⁻¹, then η Aql is indeed an SB1 Cepheid with a very long orbit of more than 50 years. Verifying whether this spectroscopic companion corresponds to the one identified in the UV imaging and spectra of η Aql by Evans et al. (2013), at a separation of 180 au, necessitates a follow-up of the long orbital period.

AQ Car was reported to show tentative signs of time-variable v_γ using an earlier set of VELOCE and literature data. In Paper I, we confirmed the SB1 nature of this star. Although weak, ∼0.7 km s⁻¹, the signal is clearly present in the VELOCE data thanks to the baseline of ∼4 yr. Although literature datasets appear imprecise to usefully inform the v_γ variations (cf. Figure B.1), in actuality, the fluctuations observed in the literature data surpass a 3-σ threshold. As σ is calculated using VELOCE data under the assumption of no orbital motion, consequently, the deviations from a stable v_γ are consistently underestimated in this scenario.

Fig. B.1. Temporal variation of vγ for AQ Car based on RVTF. The data points after BJD > 56000 are from VELOCE.

AQ Pup was reported to have a photometric companion (Evans & Udalski 1994) that was not previously confirmed spectroscopically (Vinko 1991). Here, we report low-amplitude v_γ variations (∼ 3.3 km s⁻¹) on a timescale of 2500 days. Whether this signature corresponds to the photometric companion remains unclear.

ASAS J064540+0330.4 – We present the first orbital fit for ASAS J064540+0330.4 with a P_orb of 901.18 ± 2.81 days and an eccentricity of 0.44 ± 0.01. The pulsational residuals of ASAS J064540+0330.4 indicate additional signals beyond the orbital ones. It is important to note that further observational follow-up is necessary to sample the orbital phase better and establish this orbit more firmly. The combined orbital and pulsational fit can be found in the supplementary m aterial.

ASAS J103158-5814.7 (GDR3 5351436724362450304) is a particularly interesting star, cf. Figure B.2 and Paper I. It exhibits three signals at once: very short-period overtone pulsation (classified as first overtone by Gaia DR3 SOS, P_puls = 1.11936 d), high-amplitude orbital motion exceeding 13 km s⁻¹, cf. Table 5, and strongly modulated pulsation amplitudes reminiscent of the “Blazhko" Cepheid V0473 Lyrae (Burki et al. 1982; Molnár et al. 2017). At present, the orbital variations do not allow the estimation of an orbital period and we will thus continue monitoring this intriguing star. Unfortunately, this star was not included in the list of Cepheids with time-series RV data in Gaia DR3.

Fig. B.2. Phase-folded Coralie RV of ASASJ103158-5814.7, the same star that exhibits orbital motion, the variable amplitude, and a variable pulsation period.

AX Cir is a known SB1 and its orbit was first derived by Petterson et al. (2004). Our P_orb estimate of ∼6285 days differs by 5% (10-σ) from that of Petterson et al. (2004). While our combined dataset, comprising VELOCE and literature data, incompletely samples orbital phases, we obtain an eccentricity that is in total agreement with the previously reported orbital eccentricity for AX Cir.

FN Vel was discovered to be an SB1 by Anderson (2013) who provided a first orbital solution. Kovtyukh et al. (2015) showed evidence of the companion using the Ca H & K lines. Here, we provide an updated orbital solution.

BP Cir is known to host a hot main-sequence companion (Arellano-Ferro & Madore 1985; Evans et al. 1992) based on IUE spectroscopy. Petterson et al. (2004) reported v_γ variations exceeding 5 km s⁻¹ and that the star could possibly have a long orbital period. Using VELOCE alone we have v_γ variations of ∼2 km s⁻¹, however, when we combined VELOCE with literature datasets from Petterson et al. (2004) and Borgniet et al. (2019) we found v_γ variations of ∼ 4 km s⁻¹ at a temporal baseline of 9000 days. We estimated BP Cir’s orbital period using VELOCE and literature datasets, and found a tentative solution with P_orb ∼ 4300 days and e ∼ 0.8. Nevertheless, due to limitations in the available data, we are unable to definitively determine whether a period of 4300 days or twice that duration provides a more accurate representation.

FO Car’s observations began in March 2016, and we find a Keplerian orbit of P_orb ∼ 1660 d. To the best of our knowledge, no previous publication⁹ has presented the case of FO Car. However, the online database of binary Cepheids (Szabados 2003) lists it as a spectroscopic binary, crediting Szabados & Berdnikov (in preparation). Russo et al. (1981) found photometric evidence of a companion using a phase shift technique involving optical and UV data. Based on the visual companion identified during the spectroscopic observations (cf. Figure A.1) we consider this as the most likely source of the phase shift. Given the angular separation, the visual companion is clearly not the cause of the 4.5 yr orbit, and hence, FO Car is a triple system composed of the Cepheid with an inner spectroscopic and outer visual companion.

MU Cep – The spectroscopic binarity of MU Cep (not to be confused with the much brighter red supergiant μ Cep) was discovered by Anderson (2013). Here we present an updated orbital solution. The previously reported P_orb ∼ 1500 ± 350 d exceeded the observational baseline. Thanks to continued follow-up of this faint star (⟨m_V⟩ = 12.3 mag), we have achieved nearly complete coverage of MU Cep’s orbital phase and obtain tight orbital constraints, with P_orb = 2028 ± 2.7 d.

R Cru is listed as an SB1 in the database by Szabados (2003), which credits unpublished work by L. Szabados and D. Bersier with this discovery. More recently, Evans et al. (2020) reported a visual companion of R Cru at a separation of 1.9 arcseconds based on HST imaging. We have acquired nearly 7 yr worth of very high-precision RV data for this bright Cepheid. The data clearly reveal low-amplitude orbital motion with P_orb ∼ 238 d, rendering R Cru the Cepheid with the shortest orbital period in the Milky Way. The template fitting routine further finds evidence of nonlinear period changes (see Figure 15 of Paper I), cf. also Csörnyei et al. (2022). The relative semimajor axis of a_rel > 274 R_⊙ implies that the Cepheid and its companion are not in contact, assuming a Cepheid radius of 47 R_⊙ computed using a P-R relation (Anderson et al. 2016b) and a nearly circular orbit. For a Cepheid mass of around 5 M_⊙ and a 2 M_⊙ A2 companion (upper limit from the IUE spectrum, Evans 1992), the mass function for a circular orbit yields a minimum inclination of > 3.3 deg, improving slightly over the limit from the mass function alone (1.6 deg). For the range of inclinations 3.3 − 90 deg, and assuming a 5 M_⊙ Cepheid, we thus find the allowed mass range for the companion to be between 0.09 − 2.0 M_⊙.

RV Sco – Proust et al. (1981) list RV Sco as a visual multiple system, with a non-variable companion of m_V ∼ 13 mag at a separation of 6.0 arcseconds. Evans (1992) did not detect any photometric companion for RV Sco within their IUE spectral analysis. Szabados (1989) suspected its binarity through its variable γ-velocity. Here, we confirm the finding as we detected significant (A_p2p ∼ 5.3 km s⁻¹) amplitude variation in RV resulting into a 4-σ detection of binarity. RV Sco is therefore a triple (SB1 + visual companion) system.

RY Sco is listed as a visual binary in Szabados (2003), with references to Proust et al. (1981), Evans & Udalski (1994). We here detect the stars’s SB1 nature for the first time, with very low-amplitude orbital RV variations of ∼0.6 km s⁻¹. Combining VELOCE with literature datasets and extending to longer baseline, we notice a v_γ variation of ∼2.5 km s⁻¹ on a timescale of 41 years. RY Sco is thus most likely a triple system with both a visual and a spectroscopic companions.

RZ Vel was reported to be an SB1 Cepheid by Szabados et al. (2013b). Szabados et al. (2013b) used a subset of VELOCE data (between BJD 2455650 – 2456400). We extended this VELOCE baseline to BJD 2459600 and combined it with literature datasets from Lloyd Evans (1980), Coulson & Caldwell (1985), Bersier (2002), Borgniet et al. (2019). Although, data from Bersier (2002) was also included in Szabados et al. (2013b), we do not detect 3-σv_γ variations for this or any other recent dataset after the year 2000.

S Mus was first identified as a binary using photometric observations (Lloyd Evans 1968). A first orbital solution was presented by Lloyd Evans (1982), with several others since. In particular, Petterson et al. (2004) provided a precise update to this solution, and Binnenfeld et al. (2022) determined the periodicity using partial periodograms directly computed from the spectra rather than RV measurements. We gathered 78 high-precision Coralie RVs over a total baseline of 6.3 yr. These data alone yield an orbital solution in agreement (to within 1 day) with Petterson et al. (2004) and Gallenne et al. (2019). As shown by Evans et al. (2015, 2020), the companion is an X-ray active B3V star.

S Nor was detected as an SB1 by Szabados (1989), where they estimated tentative orbital periods of 3300 or 6500 days. Groenewegen (2008) modeled S Nor’s RV data using a circular orbit with a fixed P_orb of 3600 days and concluded that more RV data are needed to better determine this orbit. Here, we do not find this orbital estimate to be consistent with VELOCE data. We do not detect any v_γ variations with VELOCE data alone, however, when VELOCE is combined with literature datasets we do obtain a 3-σ detection for the oldest RV datasets. However, the only cluster that exceeds 3-σ is based on only three observations and the uncertainty of this particular cluster is likely underestimated. Consequently, we doubt the validity of the only > 3-σ outlier and cannot consider the star a bonafide binary. Nevertheless, if the RVTF results are taken at face value, then the orbit would have to be significantly eccentric. Lastly, the v_γ residuals indicate that a potential orbit clearly has high eccentricity.

SU Cyg is a solitary case wherein the pulsation model was derived using a combination of VELOCE data and zero-point corrected RV data from Borgniet et al. (2019). While VELOCE data provided strong constraints on the orbit, they were insufficient on their own to create a satisfactory Fourier series model. To address this, we included the Borgniet et al. (2019) measurements to formulate the RV template. This template was then used to compute v_γ¹⁰. Finally, employing the MCMC method, we obtained the V+L orbit of SU Cyg with P_orb of ∼548 ± 0.1 d and eccentricity of ∼0.34 ± 0.005. These orbital estimates are in agreement with ΔP_orb < 1 d (all orbital elements within 3-σ difference) the ones in the literature from Evans (1988), Imbert (1984), Groenewegen (2008).

TX Mon was first reported as an SB1 by Szabados & Pont (1998), who estimated a short P_orb ≈ 470 ± 30 d. We here confirm the SB1 nature of the star and determine the orbit, albeit at much longer P_orb = 1607 ± 1 d, with K = 11.9 ± 0.04 km s⁻¹.

V0659 Cen is host to a hot B6.0V companion detected directly using IUE UV spectroscopy (Evans 1994). V0659 Cen also has a high X-ray flux potentially originating from its companion (Evans et al. 2022). Here we present the first RV evidence of the Cepheid’s orbital motion with a peak-to-peak RV variation of ≳5 km s⁻¹. Combining VELOCE and literature data, we estimate a preliminary orbital period P_orb ∼ 9300 ± 130 d with mild eccentricity (e = 0.32). We designate this orbit as tentative due to the ambiguity regarding whether the orbital period is approximately 6000 or 9000 days. It should be emphasized that the current uncertainties, which are around 100 days, do not account for this ambiguity. Additional HST UV spectroscopy will be particularly valuable to measure the mass of the Cepheid (Evans et al. in prep.).

V1334 Cyg is among the few sample cases where there is a disagreement between V and V+L estimates, specifically of 0.8% (3-σ) in P_orb for V1334 Cyg. Even after combining the VELOCE and literature datasets, constraints on the ascending part and extremities of the orbital phase remain limited. As the literature dataset from Gorynya et al. (1992) and Borgniet et al. (2019) add extra constraints on the descending branch of the orbital phase, we prefer the V+L orbital estimate for V1334 Cyg. Furthermore, both V and V+L orbital estimates are in 0.5% (3-σ) disagreement with the ones from Gallenne et al. (2018) and 3-σ and 2-σ, respectively, with the ones from Evans (2000). This maybe attributed to the aforementioned poor sampling at some orbital phases. Future observations, with overall improved sampling will help to elucidate the discrepancy. Lastly, Gallenne et al. (2018) demonstrated that Gaia DR2 parallax is off by 3.6-σ compared to the orbital parallax determined from the visual orbit + RVs of both components. This could be another reason for the discrepancy between our orbital parameters and the ones from Gallenne et al. (2018).

VY Car was considered a tentative SB1 candidate by Anderson et al. (2016a). Despite improved treatment of instrumental zero-point offsets, the orbital signal has neither disappeared nor clarified. We note the RVTF analysis suggests a trend with time. However, none of the epochs deviate from the average by more than 3-σ. Hence, VY Car remains an unconfirmed SB1 candidate with a potentially very long orbital period.

VY Per was reported to exhibit v_γ variations exceeding 5 km s⁻¹ by Szabados (1992) over a baseline of two years. We confirm the SB1 nature of this star and obtain a tentative orbital solution, with a P_orb ∼ 806 d and e ∼ 0.4. Additionally, we also derived its MCMC (V+L) orbital estimate The MCMC analysis indicates that the true orbital period could be a multiple of 805 d. Unfortunately, the result remains ambiguous due to the sampling of the data, which also does not include the orbit’s descending branch.

W Sgr’s P_orb using VELOCE data differs by 6-σ, 20-σ and 18-σ from that of Babel et al. (1989), Petterson et al. (2004), and Gallenne et al. (2019) respectively. We would like to note that in our work the eccentricity is not sufficiently constrained, hence, the orbit is probably more eccentric and shorter than our model as well as the P_orb is potentially more uncertain than the fit suggests. Unfortunately, the current dataset does not provide sufficient information for a more accurate fit, and additional observations will be necessary. In light of these considerations, we emphasize that Petterson et al. (2004) currently offer a more reliable orbital solution.

X Cyg was reported to exhibit an orbital signal via the light time travel effect by Szatmary (1990), who estimated P_orb ≈ 55 yr. Evans et al. (2015) and Hintz et al. (2021) did not detect long-term v_γ variations based on ∼ 11-year baselines. We here identify X Cyg’s v_γ variations using much longer temporal baselines of ∼ 40 yrs by combining VELOCE data with the literature (Barnes et al. 1987; Wilson et al. 1989; Kiss 1998; Barnes et al. 2005; Gorynya et al. 1992; Storm et al. 2004; Bersier et al. 1994; Borgniet et al. 2019). We thus find a total v_γ difference of ∼ 3 km s⁻¹ and a lower limit of 40 yr on P_orb, consistent with Szatmary (1990).

X Pup was mentioned by Anderson et al. (2016a) to exhibit tentative signs of time-variable v_γ based on a difference between VELOCE data and the literature. However, the improved corrections of RV zero-points from VELOCE-I together with the clustered RVTF analysis do not reveal a significant signature of orbital motion.

XX Cen – Groenewegen (2008) determined an orbital solution with P_orb = 924 ± 1.1 d and K = 4.47 ± 0.28 km s⁻¹, and e fixed to 0. We find a significantly different orbit for XX Cen based on VELOCE data alone and determine an orbital solution with P_orb = 722.44 ± 0.42 d, K = 5.4 ± 0.03 km s⁻¹, that is significantly eccentric with e = 0.295 ± 0.004.

All Tables

All Figures

	Fig. 1. Location of the sample SB1 Cepheids on the sky in the equatorial (top panel) and galactic (bottom panel) coordinate system. The symbols are color coded in the Gaia G magnitude. The symbols surrounded with blue circles are systems for which we provide an orbital solution in the current work.
In the text

Fig. 2. Pulsation and orbital fit of the naked-eye Cepheid prototype δ Cephei. Left top: pulsation RV variability where the systemic velocity, vγ, 0 is indicated by a dashed line. Right top: residuals after fitting the Fourier series and Keplerian orbit against observation date. Left bottom: phase-folded orbital RV variation against orbital phase. Right bottom: orbital RV variation versus the observation date.

In the text

	Fig. 3. Pulsation and orbital fit of one of the faintest binaries in our sample, NT Puppis. The figure description is the same as Figure 2.
In the text

	Fig. 4. Orbit of U Vul determined by fitting vγ epoch residuals of literature and VELOCE data.
In the text

	Fig. 5. Joint posterior distributions derived from the MCMC sampling for U Vul. The histograms show marginalized posterior distributions of the parameters.
In the text

Fig. 6. Distribution of Ap2p, maximum amplitude variation in vγ, for the stars in our sample where we could fit a polynomial to the vγ trend (see Section 3.3 for details). In the left panel, we have the Ap2p vs the baseline in time over which it was calculated and in the right panel, the same against the pulsation period. The newly discovered SB1s from our sample are plotted in filled squares, while the literature-known SB1s are plotted in filled circles. The color bar represents the degree of polynomial fitted in the Fourier+Polynomial model.

In the text

	Fig. 7. Zero-point corrected Δvγ for new SB1s detected through RVTF. The dotted lines mark the standard deviation of the vγVELOCE. This figure only includes stars where we have a 3-σ detection through our analysis. The color bar represents the S/NSB1 of the various vγ measurements where the S/NSB1 was calculated as described in the text.
In the text

	Fig. 8. Literature-known SB1 candidates detected through RVTF. The figure description is the same as Figure 7.
In the text

	Fig. 9. Stacked histogram presenting the new VELOCE SB1 detections (left panel) and the new VELOCE orbits (right panel) along with the same thing from Szabados (2003). SY Nor’s orbit is provided by us as well as Bersier (2002); however, it is missing in the Szabados (2003) database.
In the text

Fig. 10. Distribution of different orbital elements (Upper left: argument of periastron ω; upper right: eccentricity e; lower left: semi-amplitude K of the orbital RV variations; and lower right: projected semi-major axis a sin i) versus orbital period Porb for Cepheids with the orbits determined here. SB1s with newly reported orbits in this study are plotted using open squares, while VELOCE orbital estimates of SB1s with literature-known orbits are plotted in open circles. Symbols with double boundaries highlight V+L orbits determined using VELOCE and zero-point corrected literature RVs. R Cru, SY Nor and YZ Car are absent in the top left panel due to insignificant eccentricity.

In the text

	Fig. 11. Distribution of Gaia DR3 RUWE with the semi-major axis (left) and orbital period (right) of our sample systems. The stars tagged as binaries from PMa (Kervella et al. 2022) are marked with crosses. The red solid line indicates a sin i of 1 au, and the dashed red lines indicates the same at 0.5 and 2 au.
In the text

	Fig. 12. Distribution of Gaia DR3 RUWE with the maximum amplitude of the (linear or nonlinear) vγ variation for SB1 systems with Trend. The stars tagged as binaries from PMa (Kervella et al. 2022) are marked with crosses.
In the text

	Fig. 13. Probability density function of Gaia DR3 RUWE for VELOCE SB1 and non-SB1 systems.
In the text

Fig. 14. Photometric amplitude ratios as indicators of companion stars. First three panels from left to right: photometric-to-RV amplitude ratios from Gaia DR3 and VELOCE against the logarithmic pulsation period, with wavelength of the photometric band increasing from left to right. Right panel: ratio of Gaia photometric amplitudes in GBp to GRp. Filled circles show SB1 Cepheids reported in VELOCE. Other single-mode Cepheids observed by VELOCE are shown as open circles. There is no clear difference between SB1 Cepheids and all others in any of the amplitudes considered, nor a dependence of such a difference on log Ppuls.

In the text

	Fig. A.1. Coralie guiding camera image of FO Car. The center of the Cepheid is dark because of the fiber placement during the integration. The object at the top left is an outer (visual) companion, which renders the SB1 Cepheid FO Car a triple system.
In the text

	Fig. B.1. Temporal variation of vγ for AQ Car based on RVTF. The data points after BJD > 56000 are from VELOCE.
In the text

	Fig. B.2. Phase-folded Coralie RV of ASASJ103158-5814.7, the same star that exhibits orbital motion, the variable amplitude, and a variable pulsation period.
In the text