Multiplicity of Galactic Cepheids and RR Lyrae stars from Gaia DR2 - I. Binarity from proper motion anomaly

Classical Cepheids (CCs) and RR Lyrae stars (RRLs) are important classes of variable stars used as standard candles to estimate galactic and extragalactic distances. Their multiplicity is imperfectly known, particularly for RRLs. Astoundingly, to date only one RRL has convincingly been demonstrated to be a binary, TU UMa, out of tens of thousands of known RRLs. Our aim is to detect the binary and multiple stars present in a sample of Milky Way CCs and RRLs. In the present article, we combine the Hipparcos and Gaia DR2 positions to determine the mean proper motion of the targets, and we search for proper motion anomalies (PMa) caused by close-in orbiting companions. We identify 57 CC binaries from PMa out of 254 tested stars and 75 additional candidates, confirming the high binary fraction of these massive stars. For 28 binary CCs, we determine the companion mass by combining their spectroscopic orbital parameters and astrometric PMa. We detect 13 RRLs showing a significant PMa out of 198 tested stars, and 61 additional candidates. We determine that the binary fraction of CCs is likely above 80%, while that of RRLs is at least 7%. The newly detected systems will be useful to improve our understanding of their evolutionary states. The discovery of a significant number of RRLs in binary systems also resolves the long-standing mystery of their extremely low apparent binary fraction.


Introduction
The remarkable correlation of the intrinsic luminosity of classical Cepheids (CCs) (Leavitt 1908;Leavitt & Pickering 1912;Fouqué et al. 2007) and RR Lyrae stars (RRLs) (Catelan et al. 2004;Neeley et al. 2017) respectively with their pulsation period and metallicity makes these two classes of variable stars essential standard candles for Galactic (Drake et al. 2013), globular cluster (Carney et al. 1992), and extragalactic distance measurements (Clementini et al. 2003;Riess et al. 2011Riess et al. , 2016).An analysis of the GDR2 parallaxes of Galactic CCs in the context of the extragalactic distance scale was recently presented by Riess et al. (2018).Classical Cepheids are intermediate-mass stars (typically 5 to 10 M ) and this class therefore comprises a high fraction of binary and multiple stars (Gieren 1982;Szabados 2003a;Evans et al. 2013Evans et al. , 2015;;Sana 2017).RR Lyrae stars Tables A.1 and A.5 are available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr(130.79.128.5) or via http: //cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/ are short-period (P ≈ 0.5 d) low-mass pulsators (m ≈ 0.6 M ) that are abundant in the Galaxy and especially in globular clusters (Bailey & Leland 1899;Pickering et al. 1901).They are old (age ≈ 10 Ga 1 ), horizontal branch stars, with a typical radius of 4 to 8 R (Marconi et al. 2005).Despite their relative faintness compared to CCs, RRLs have been used to measure distances in the Local Group and beyond (Da Costa et al. 2010;de Grijs et al. 2017;Monelli et al. 2018), and their ubiquity makes them important standard candles for Galactic astronomy (Dékány et al. 2018;Contreras Ramos et al. 2018).Approximately two hundred thousand variable stars are classified as RRLs from groundbased surveys or the Gaia DR2 (Soszyński et al. 2014, Clementini et al. 2018, Muraveva et al. 2018, Holl et al. 2018, Rimoldini et al. 2018).It is remarkable, however, that there is to date very little evidence for RRLs in binary systems, as only one case has been convincingly identified: TU UMa (Liška et al. 2016a).
From the analysis of the light curves of a large sample of nearly 2000 RR Lyrae stars, Hajdu et al. (2015) identified 12 binary candidates that display phase shifts of their light curves that can be attributed to light-time effect (LiTE) that point to the presence of an orbiting companion.Liška et al. (2016b) presents a study of 11 systems searching for LiTE in O-C diagrams, and in the RRLyrBinCan 2 database of candidate binary RRLs.A search has also been conducted by Guggenberger & Steixner (2015) in the Kepler mission light curve database, but without detection.Sódor et al. (2017) has interpreted the Kepler light curve phase modulations of KIC 2831097 as being caused by the presence of an orbiting black hole of 8.4 M , but the radial velocity data contradict the binary interpretation.The common presence of the Blazhko effect (Blažko 1907;Jurcsik et al. 2011;Szeidl et al. 2012;Jurcsik & Hajdu 2017) and period drifts (Szeidl et al. 2011) in many RRLs complicates the uniqueness of the interpretation of the observed phase shifts.Soszyński et al. (2011) found a likely candidate for an RRL in a 15.2-day period eclipsing binary system (OGLE-BLG-RRLYR-02792).It was subsequently identified as a peculiar type of "RR Lyr impostor" (Pietrzyński et al. 2012;Smolec et al. 2013), and was included in a new class of binary evolution pulsators (BEP) that are believed to be very rare (Karczmarek et al. 2017).
The companions of CCs and RRLs are important for several reasons (Szabados et al. 2010).They may influence their evolution through mass transfer.They also shift the apparent brightness of their parent stars in a systematically positive way by up to 10% or more in the visible (Gallenne et al. 2013), which affects the zero point of the CC Leavitt law.Anderson et al. (2016a) showed that companions have a limited effect on the observational properties of the brighter, long-period CCs, whose usefulness as distance indicators is thus not affected (see also Anderson & Riess 2018).However, due to their lower intrinsic brightness, short-and intermediate-period CCs are more likely to exhibit a significant relative photometric contribution from the companions, particularly at short wavelengths.Companions of CCs are often hot main sequence (MS) stars, therefore making the CCs appear bluer.This consequently biases the estimate of their color excesses and reddenings.If not taken into account, their orbital displacement affects the trigonometric parallax measurements.
The census of the companions of CCs and RRLs is incomplete, due to the high contrast between the bright pulsators and their companions.Stellar population synthesis models by Neilson et al. (2015) predict that the binary fraction of CCs is likely lower than for their MS progenitor, and that about half of the CCs are products of binary interactions.This would be caused by interactions between the close-in companions and the Cepheid progenitors while they evolve on the red giant branch.Sana et al. (2012) claims that 70% to 100% of O stars have companions, whereas Neilson et al. (2015) predict that only 35% of Cepheids do.The binary fraction of Galactic CCs is thus a key observable to test the intermediate-mass star formation and evolution scenarios.
The goal of the present work is to test for the presence of close-in companions of Galactic CCs and RRLs using the Gaia Second Data Release (hereafter GDR2;Gaia Collaboration et al. 2016, 2018a).Our CC and RRL samples are presented in Sect. 2. In Sect. 3 we use the GDR2 position and proper motion (PM) measurements together with the Hipparcos catalog (Perryman et al. 1997;van Leeuwen 2007) to search for PM anomalies.In the companion Paper II (Kervella et al. 2019) we search the GDR2 for common PM stars located near the CCs and RRLs, and we test the possibility that they are gravitationally bound.We postpone the discussion of individual stars to Paper II.

Selected samples
We present in this section the sample of CCs and RRLs selected for our present PM analysis (Paper I) and for the search for resolved common proper motion companions presented in Paper II.In particular, we detail our choice of parallax values and the systematic corrections that we applied to the different data sets (Sect.2.3).

Cepheids
We chose the sample of 455 Galactic CCs assembled by Berdnikov et al. (2000).We uniformly adopted the CC parallaxes from the GDR2, which we corrected following the procedure detailed in Sect.2.3, except for four stars (U Aql, R Cru, SU Cru, and Y Sgr) for which we adopted the Hipparcos parallax from van Leeuwen (2007).For δ Cep, we adopted the GDR2 parallax of its physical companion δ Cep B, as detailed in Paper II.For RY Vel, whose GDR2 and Hipparcos parallaxes are negative, we adopted the photometric distance of Berdnikov et al. (2000) based on multicolor period-luminosity relations, renormalized to the LMC distance modulus established by Pietrzyński et al. (2013), giving = 0.39 ± 0.06 milliarcseconds (mas).We add to the sample the short-period double-mode pulsator Y Car, which is a known triple system (Evans et al. 2005).We adopt the distance modulus of µ = 10.8 ± 0.3 determined by Evans (1992), corresponding to a parallax of = 0.69 ± 0.10 mas, i.e., with a ±15% uncertainty.Although the membership of Y Car to the open cluster ASCC 60 is listed as inconclusive by Anderson et al. (2013), the distance of the cluster (1.1 kpc) determined by Kharchenko et al. (2016) does not exclude this possibility.The GDR2 parallax value ( = 0.301 ± 0.035 mas) is likely unreliable, possibly due to the astrometric wobble of the center of light of the system.
Out of the 455 CCs present in the Berdnikov et al. (2000) catalog plus Y Car, 254 are present in the Hipparcos catalog and were tested for the presence of a proper motion anomaly (hereafter PMa).The remaining stars are usually fainter than the Hipparcos magnitude limit.

RR Lyrae
We extracted the RR Lyrae type variables from the General Catalogue of Variable Stars (Samus et al. 2017), which comprises 8509 stars.Only 198 of these stars are present in the Hipparcos catalog and therefore suitable for the search for companions from their PMa.We adopt the GDR2 parallaxes of RRLs, uniformly corrected following the procedure detailed in Sect.2.3.For RR Lyr itself, which is absent from the GDR2 catalog, we adopt the TGAS parallax from the GDR1 of [RR Lyr] = 3.64 ± 0.23 mas (Michalik et al. 2014;Gaia Collaboration et al. 2016).
We processed all the stars present in the selected catalogs (456 CCs and 198 RRLs), but the variability class of some targets is incorrect, and we present the results related to these objects separately from CCs and RRLs.

Gaia DR2 basic corrections and quality control
The GDR2 parallaxes are affected by a mean global zero point (ZP) offset (Lindegren et al. 2018).Examples of determinations of the GDR2 ZP include for instance the work by Riess et al. (2018), who derived a value of −46±13 µas specifically for CCs, and Muraveva et al. (2018) who obtained −56 ± 6 µas for RRLs.Arenou et al. (2018) list a statistically identical ZP value to that of Muraveva et al. (2018) for the full sample of GDR2 RRLs (−56 ± 5 µas; their Table 1).However, their ZP for the restricted sample of RRLs present in the General Catalogue of Variable Stars (GCVS; Samus et al. 2009) is −33±9 µas.For CCs, Arenou et al. (2018) obtain a ZP offset of −32 µas.
The choice of ZP does not affect the PM anomaly and consequently the detected binaries in the present paper.It has however an influence on the masses of the companions and their linear orbital radii, which are inversely proportional to the parallax.The determination of the ZP of GDR2 is a complex question that is beyond the scope of the present work.We therefore systematically corrected the GDR2 parallaxes of CCs and RRLs by adding a constant ∆ G2 = +29 µas offset to the catalog values, as recommended by Lindegren et al. (2018) and Luri et al. (2018).This value corresponds to a sky average derived from quasar measurements, and is compatible with the ∆ G2 obtained by Arenou et al. (2018) for CCs and for RRLs.A future revision of the GDR2 ZP to a new value ∆ * G2 can be used to correct the determined companion masses m 2 to new values m * 2 through the simple multiplication where G2 is the uncorrected parallax from the GDR2 catalog.We note that choosing an offset correction of ∆ G2 = +56 µas instead of +29 µas has a negligible impact on all the CC and RRL companion masses (Sect.3.4 and 3.5) within their error bars.
We implemented the correction of the parallax uncertainties described in Eq.A.6 of Lindegren et al. (2018), as recommended by Arenou et al. (2018).We corrected the GDR2 PM vectors for the rotation of the Gaia reference frame (Gaia Collaboration et al. 2018b) reported by Lindegren et al. (2018) using the expressions µ α,corr = µ α + w x sin(δ) cos(α) + w y sin(δ) sin(α) − w z cos(δ), (2) µ δ,corr = µ δ − w x sin(α) + w y cos(α), (3) where w x = −0.086± 0.025 mas a −1 , w y = −0.114± 0.025 mas a −1 , and w z = −0.037±0.025mas a −1 .As discussed by Lindegren et al. (2018) the systematic uncertainty on the GDR2 PM vectors is limited to σ sys (µ) = 66 µas a −1 per component for small separations (see also Arenou et al. (2018) and Luri et al. 2018).This is possibly lower in reality, but we conservatively added quadratically this systematic uncertainty to the stated PM error bars of both RA and Dec axes.Although this is not a requirement for the present Paper I, we corrected the G-band magnitudes using the expression G corr = 0.0505 + 0.9966 G from Casagrande & VandenBerg (2018) in view of the calibration of the resolved companion magnitudes in Paper II.The validity of this correction is demonstrated over the range 6 G 16.5, but we also apply it to fainter stars.The amplitude of the correction is at most 30 mmag for a G = 6 star, and therefore of marginal importance for our purpose.The photometry of the stars with G < 6 is unreliable due to saturation, but we did not find such bright candidate companions.
We tested the GDR2 record of the stars of our samples following the three quality criteria defined by Arenou et al. (2018) (their Sect.4.1): (1) a reduced χ 2 of the Gaia astrometric fit below a limit dependent on the G magnitude (e.g., χ 2 red < 8 for a G = 10 magnitude star), (2) a photometric G BP −G RP flux excess factor within acceptable color-dependent limits, and (3) more than six visibility periods.The stars that do not satisfy these three criteria are flagged with a † symbol in Table A.1.We also computed the reduced unit weight equivalent noise (RUWE4 ; denoted in the following, see also Kervella et al. 2018) of the GDR2 astrometric solution.The RUWE is a combination of the astrometric χ 2 , the number of good observations N, the G magnitude, and the color index C = G BP − G RP .In Table A.1 and the following the stars for which > 1.4 (as recommended by Lindegren) are flagged with a ‡ symbol.
Binary stars present a natural discrepancy in the χ 2 of their astrometric model fit, due to the present assumption in the GDR2 astrometric model that all stars are single.As a consequence, they are more likely to not fulfill quality criterion (1) of Arenou et al. (2018) and exhibit > 1.4.These quality indicators can thus be viewed as de facto indicators, however imperfect, of astrometric binarity.To prevent the rejection of actual binary stars, we therefore kept all the stars in our analysis, including those with quality flags (we provide the flag information).Future Gaia data releases will include the binarity in the astrometric fit, and therefore provide a separate view of the contributions of binarity and instrument noise to the χ 2 .

Proper motion anomaly
The principle of our search for close-in orbiting companions is to look for a difference in PM vector between the mean PM computed from the Hipparcos (1991.25)and GDR2 (2015.5)astrometric (α, δ) positions on the one hand (hereafter µ HG , for Hipparcos-Gaia) and the individual PM vectors µ Hip and µ G2 respectively from the Hipparcos and GDR2 catalogs on the other hand.This approach to compare the long-term to short-term PM vectors has historically been employed by Bessel (1844) to discover the white dwarf companion of Sirius.It was also applied recently to various types of stars by Wielen et al. (1999), Jorissen et al. (2004), Frankowski et al. (2007), Makarov et al. (2008), Brandt (2018), Brandt et al. (2018), Kervella et al. (2018), andSnellen &Brown (2018).Figure 1 shows the definition of the different PM vectors considered in the present work.The combination of Gaia and Hipparcos data has already been used after Gaia DR1 (Lindegren et al. 2016) to produce the TGAS catalog (Michalik et al. 2014(Michalik et al. , 2015)).A description of the sources of uncertainty to take into account in the combination of these two catalogs is presented by Lindegren (2018).
We identify µ HG to the projected velocity vector of the center of mass, while µ Hip and µ G2 represent the projected velocity vector of the photocenter of the system at the Hipparcos and GDR2 epochs, respectively.For single stars that have a linear uniform space motion, these three projected vectors have the same constant direction and norm (neglecting the variable spherical projection effects for such distant stars).The presence of an orbiting companion will displace the photocenter away from the center of mass, due to the difference between the mass ratio and the flux ratio of the two stars.In this case, the photocenter will revolve around their center of mass following a "virtual orbit" with a semimajor axis a where a is the semimajor axis of the physical orbit of the primary star around the center of mass, L 1 its flux, and L 2 the flux of the secondary component.As the Hipparcos and Gaia missions measure the PM of the photocenter, a deviation will appear with the PM of the center of mass.
The photocenter of a binary system comprising a CC is usually very close to the CC due to the high brightness of supergiants compared to their companions, that are usually MS dwarfs (L 2 L 1 ).For RRLs, the flux of the companion stars is also small; although the RRLs are less luminous than CCs, their companions are also significantly fainter (compact objects, very lowmass dwarfs) than their CC counterparts as they are very old stellar systems.In the following, we uniformly assume that the photocenter of the system is coincident with the position of the CC or RRL.This assumption results in a systematic underestimation of the true tangential orbital velocity of the Cepheid by a factor L 1 /(L 1 + L 2 ).
We define the signal-to-noise ratio of the PMa of the Hipparcos/GDR2 measurements with respect to the mean PM µ HG as where C = 2 ρ σ µ Hip/G2 σ µ HG corresponds to the correlation term (with a degree of correlation ρ) between µ HG and the PM vectors from the Hip/GDR2 catalogs.These two quantities are correlated as the astrometric positions (α, δ) in the Hip/GDR2 catalogs, which are used to compute µ HG , and are themselves correlated to the µ Hip/G2 PM vector coordinates.However, since the position uncertainty intervenes in its computation with a divisive factor 24.25 (difference in years between the Hipparcos and GDR2 epochs), C is much smaller than the µ 2 Hip variance in ∆ Hip and than the Hipparcos positional variance in ∆ G2 , and can thus be neglected in both cases.For the identification of candidate binary stars, we consider the maximum of the two values ∆ = max(∆ Hip , ∆ G2 ).In general, for a given star showing a PMa, the signal-to-noise ratio ∆ G2 is significantly higher than ∆ Hip thanks to the higher accuracy of Gaia (typically by one order of magnitude).However, depending on the configuration of the orbit and the orbital phase, a PMa may be detectable at the Hipparcos epoch and not at the GDR2 epoch.

Levels of analysis
Several levels of analysis can be achieved, depending on the available observational constraints: 1. Proper motion anomaly only: Knowing the parallax, ∆µ gives the 2D tangential linear velocity v tan of the target in the center-of-mass referential at the measurement epoch.This is a projection of the true 3D velocity vector of the star on the plane of the sky (i.e., the plane perpendicular to the line of sight containing the star), and therefore its norm is a lower limit of its orbital speed.With an a priori estimate of the mass of the target (Sect.3.2.2) and an additional hypothesis on the mass ratio q of the binary (Sect.3.2.3),we derive a range of maximum semimajor axes and orbital periods using the expressions Here we implicitly assume that the orbit is circular.For this simplified analysis, we considered only the PMa vectors determined from the GDR2.The PMa vectors from the Hipparcos catalog generally have one order of magnitude lower accuracy than those computed from the GDR2, and therefore provide limited constraints on the orbital radii and orbital periods of the companions.2. Proper motion anomaly and radial velocity: When the parameters of the spectroscopic orbit (P, e, ω, K) are known, the knowledge of the orbital radial velocity v r of the target at the same epoch as the PMa vector, together with the parallax, gives its complete 3D velocity vector The availability of two orbital velocity vectors, v Hip and v G2 , gives access to the inclination i of the orbital plane and the longitude of the ascending node Ω through a cross product This resolves the sin(i) degeneracy, thus allowing us to derive the full set of orbital parameters.With a prediction of the P. Kervella et al.: Multiplicity of Galactic Cepheids and RR Lyrae stars from Gaia DR2 -I primary mass (Sect.3.2.2), the mass m 2 of the secondary star can then be determined.
We refer in the following to the level of analysis that we can achieve on a given target using, e.g., "level 2" to designate the systems with two PMa vectors and the corresponding radial velocities.
The Hipparcos observations were conducted between 7 November 1989 and 18 March 1993 (Perryman et al. 1997), i.e., covering 1227 d.The GDR2 catalog values are based on data collected between 25 July 2014 and 23 May 2016 (Gaia Collaboration et al. 2018a), i.e., covering 668 d.This means that the PM vectors from these two catalogs are not instantaneous, but represent a weighted average over these observing windows that depends on the distribution of the individual observed transits.For binary systems with orbital periods shorter than these observing windows, the measured PMa is still a valuable tracer of binarity, but due to the integration of more than one orbital cycle, the PMa vector coordinates are smoothed by the observing window.As a consequence, the determination of the orbital parameters together with the spectroscopic orbit may be biased (see, e.g., the case of S Mus discussed in Sect.3.3).For short-period companions, the error bars of the Hipparcos and GDR2 PM vectors incorporate the residual wobble due to the several orbital cycles covered during the observing window (Kervella et al. 2018).The smoothing of the PMa signal also results in a significant decrease in the sensitivity of this indicator to orbiting companions for orbital periods 1000 days.The inclusion of binary fitting in the astrometric solution of future Gaia data releases will allow this limitation to be waived (see, e.g., the recent work by Snellen & Brown 2018 using Hipparcos epoch astrometry).
Cepheid binary systems with fully determined orbital parameters from the combination of visual (from classical imaging or optical interferometry) and spectroscopic orbit are still rare.To date only V1334 Cyg has a fully determined, high-precision set of orbital parameters (Gallenne et al. 2018b) including the masses of both components to 3% accuracy and their distance to 1% accuracy.This favorable configuration provides a stringent test of the reliability of the PMa analysis.In addition, Gallenne et al. (2018a) recently obtained high-accuracy interferometric astrometric orbits from interferometry for U Aql and S Mus.We briefly discuss them in Sect.3.3 together with V1334 Cyg.

A priori mass estimates
The masses of the CCs were approximated using a combination of the theoretical period-luminosity-radius relation for fundamental mode pulsators by Caputo et al. (2005) and the periodradius relation calibrated by Gallenne et al. (2017).We did not "fundamentalize" the periods of the first overtone pulsators.This is a very simple approach, but it provides sufficiently accurate estimates of the Cepheid masses for our purpose (companion mass and escape velocity estimates, see Paper II).For the short-period, first overtone pulsator V1334 Cyg (P = 3.33 d), the agreement between the prediction (4.6 M ) and the determined mass by G18 (4.29 M ) is satisfactory (8%).For the fundamental mode long-period pulsator Car, the agreement is good between the predicted value (8.4M ), the range of 8 − 10 M defined by Neilson et al. (2016) and the 9 M estimate given by Anderson et al. (2016b).For Polaris Aa, which is a first overtone pulsator, the predicted mass is m[Polaris Aa] = 4.8 M , significantly lower than the 7 M estimate by Anderson (2018).The estimate, however, was derived assuming the HST/FGS parallax from Bond et al. (2018), which is underestimated (Engle et al. 2018), and therefore the mass is likely overestimated.From the astrometric monitoring of the orbit of Polaris B, Evans et al. (2018) derive a value of m[Polaris Aa] = 3.45 ± 0.75 M , which is lower than our estimate but compatible within the uncertainties.In the following analysis, we adopt a conservative ±15% uncertainty on the predicted masses of the CCs of our sample.
For RRLs, we adopt a uniform mass of 0.6 ± 0.1 M (±15%) independent of the period.This conservatively covers the full range of possible masses predicted by the mass-metallicity relation of Jurcsik (1998): (8)

Mass ratio
Mass ratios of known multiple CCs are reviewed by Evans et al. (2015).They are distributed mostly uniformly between 0 and 1.The CCs with determined companion masses usually have mass ratios q < 1 for MS companions, as the CC has to be more evolved than the companion.However, mass ratios larger than one are exceptionally possible when the companions are themselves binary systems, for example AW Per (Evans et al. 2000;Griffin 2016).
For RRLs, the only binary known with confidence is TU UMa, for which Liška et al. (2016a) estimate m 1 = 0.55 M and obtain a minimum mass for the companion of m 2 = 0.34 M , which corresponds to a mass ratio q = m 2 /m 1 = 0.6.We note however that the true mass of TU UMa B is significantly higher than this minimum value (Sect.3.5).
In absence of spectroscopic orbital parameters, we assume q = 0.5 ± 0.3 for CC and RRL companions in the following discussion, with the pulsating star being the more massive.

Validation on V1334 Cyg
V1334 Cyg is a short-period (P = 3.33 d) first overtone CC (Evans 2000;Gallenne et al. 2013) that is a known spectroscopic and interferometric binary system (Evans 1995(Evans , 2000;;Gallenne et al. 2013).The full set of orbital parameters, masses and distance of V1334 Cyg have been determined with very high accuracy by G18.This therefore provides us with an excellent test system (see Sect. 3.2) to validate our approach based on PM anomalies.
Figure 2 shows the orbits of the two components around the barycenter from G18, as well as the Hip and GDR2 tangential velocity vectors (PM anomalies).The agreement in position angle of the PM vectors with respect to the expected directions is satisfactory.The GDR2 vector is consistent in terms of norm, but the Hip vector's norm is slower than expected.To conduct a blind analysis we considered as input parameters only the spectroscopic orbital parameters determined by Evans (2000).The parameters that we derive are listed in Table 1, together with the corresponding values found by G18 for comparison (in parentheses).The agreement is good on the inclination of the orbital plane i (1.2σ) and the longitude of the ascending node Ω (0.8σ).The determined companion mass m 2 is also in good agreement (0.3σ).
The i and Ω parameters are directly determined from the radial and PMa vectors at the Hip and GDR2 epochs.They are usually impossible to estimate without spatially resolving the system.This good consistency of the results between two fully independent approaches demonstrates the high potential of the PMa signal to determine orbital parameters from Gaia measurements.(2018b).The virtual orbit of the photocenter of the system is shown as a gray ellipse.The measured tangential velocity vector (proper motion anomaly) is represented at the Hipparcos and Gaia epochs.
It is interesting to note that we considered in this analysis that the photocenter is perfectly coincident with the Cepheid.For V1334 Cyg, we know from Gallenne et al. (2018b) that its virtual orbit (Fig. 2, gray ellipse) is ≈ 20% smaller than that of the CC (orange ellipse).Correcting a posteriori for this offset results in an increase of the companion mass to m 2 = 3.9 ± 0.6 M , within 0.2σ of the true mass of V1334 Cyg B (4.0 M ).The inclination i is increased to 119 ± 6 • , which is also within 1σ of the value determined by Gallenne et al. (2018b).This confirms that our hypothesis that the orbit of the photocenter is identical to that of the CC results in a systematic underestimation of the mass of the companions (Sect.3.1).However, V1334 Cyg is an extreme case as the relative brightness of V1334 Cyg B in the visible is not negligible (≈ 10%).For most of the CC binaries considered here, the companions are much fainter, and the bias on the determined companion masses is negligible.Gallenne et al. (2018a) recently reported an orbital solution for U Aql and S Mus based on astrometric measurements obtained by interferometry.They derived companion masses of m 2 = 2.2±0.2M and 4.0±0.2M , respectively for the two CCs.For U Aql the agreement with our estimate of m 2 = 1.9 ± 0.3 M (0.8σ; Table 2) is good.For S Mus we obtain a companion mass of m 2 = 2.2 ± 0.3 M , significantly lower (5σ) than the Gallenne et al. (2018a) value.This discrepancy arises from the orbital period of this system, which is significantly shorter (P orb = 506 d) than the GDR2 and Hipparcos measurement windows.As discussed in Sect.3.2.1, this biases the corresponding PMa estimate, as testified by the classification of S Mus as "preliminary" in Table 2.

PM anomalies of Cepheids
Tables A.1 to A.3 list the result of the search for PM anomalies on our selection of CCs.We identify 31 stars with a high ∆ > 5, 26 stars with 3 < ∆ < 5, and 75 additional CCs with indications Table 1.Parameters of the V1334 Cyg system from the combined analysis of the spectroscopic orbit of Evans (2000) and the proper motion anomaly vectors.The high accuracy values derived by Gallenne et al. (2018b) are given for each parameter in parentheses.cally related to the CC (Paper II), and we therefore removed this star from our binary star count.We note that SU Cas and δ Cep are likely triple and quadruple systems, respectively, and Polaris, which is the nearest CC (but is not part of our sample), is a triple system (Paper II).

Adopted parameters
Considering the 100 nearest CCs in our sample (with > 0.56 mas), 32 stars show ∆ > 3, and 31 others are known binaries from the Szabados (2003b) database.Four CCs in this sample (TV CMa, ER Car, V0532 Cyg, and V0950 Sco) are not classified as binaries by Szabados (2003b) and have ∆ < 3, but we report resolved gravitationally bound companions in Paper II.Altogether, we therefore obtain a minimum binary fraction of P = 67% for this sample of 100 nearby CCs.Another way to estimate the binary fraction is to rely on an estimate of the completeness level r of the PMa binary detections within our 100 CC sample.An approximation of r is provided by the fraction of CCs with ∆ > 3 among the known binary CCs.We obtain a value of r = 43% for our sample that characterizes the mean efficiency of the PMa analysis to detect known CC binaries.Applying this ratio to the detected PMa with ∆ > 3 gives an extrapolated binary fraction of P = 31/0.43= 72%, which is consistent with the overall minimum binarity previously determined.The smoothly decreasing shape of the binary fraction curve shown in Fig. 4 (right panel) is due to the the sensitivity of the PMa technique in terms of companion mass being a linear function of the parallax.This can be observed, for instance, by restricting our sample to the 50 closest CCs.This sample contains 24 stars with ∆ > 3, and we derive r = 49%, hence an extrapolated binary fraction of P = 98% (in agreement with Fig. 4 for nearby stars).
Considering the minimum P values above, and taking into account the decreasing sensitivity of the PMa analysis with distance, we conclude that the actual binary fraction of CCs is probably above 80%.This fraction is consistent with the estimate by Szabados (2003b) (P 80%), but higher than observed by, among others, Anderson et al. (2016a)  The orbital parameters and companion masses that we derive from the level 2 analysis of the systems that have spectroscopic orbits are presented in Table 2.The derived companion mass m 2 is inversely proportional to the adopted parallax of Fig. 5. Left: Histogram of the mass ratios q = m 2 /m 1 of the Cepheid systems with spectroscopic orbits.Right: Distribution of the companion masses m 2 as a function of the orbital semimajor axis.The "preliminary" systems in both panels correspond to the systems with uncertain parameters (lower part of Table 2).
the system.As a consequence, a revision of the Gaia parallaxes in the future data releases will result in a revision of the companion masses.We note a good agreement of our values of Ω for FF Aql and W Sgr with the HST-FGS determinations by Benedict et al. ( 2007), but a difference in the inclinations i, which may come from a different definition of this parameter.Anderson et al. (2015) recently announced the discovery of a close-in companion of the prototype CC δ Cep.Adopting their spectroscopic orbital parameters, the mass m 2 = 0.72 ± 0.11 M that we obtain for this companion is in good agreement with the range of 0.2 − 1.2 M estimated by these authors.We note, however, that due to its high brightness, the reliability of the GDR2 astrometric solution of δ Cep is uncertain, and its PM vector may therefore be biased.This companion mass should therefore be confirmed using a more accurate PM vector from a future Gaia data release.We list in Table 2 the results for all targets, but the shorter periods ( 1000 d) should be considered preliminary due to the PM vector smearing over the Hipparcos and GDR2 observing windows (Sect.3.2.1).A histogram of the mass ratios q = m 2 /m 1 is presented in Fig. 5 (left panel).We observe a high frequency of relatively low-mass companions, with a median mass ratio q = 0.4.This is consistent with the statistical estimates by Moe & Di Stefano (2017) for CCs in binaries, which are based on the observational results by Evans et al. (2013) and Evans et al. (2015).The distribution of the masses as a function of the orbital semimajor axis is shown in Fig. 5 (right panel).The limited number of systems with a semimajor axis larger than 10 au shows that the radial velocity technique has a higher sensitivity to short orbital periods.For the same companion mass, the astrometric detection technique is more sensitive to long periods and is therefore very complementary.
Selected properties of binary systems with CCs were compared with the synthetic population of 30 000 solar metallicity CCs in binaries, generated with the population synthesis code StarTrack (Belczynski et al. 2002(Belczynski et al. , 2008)).The simulation confirms that the mass ratios (from less to more massive stars) form a uniform distribution over the entire Cepheid mass range (5 − 10 M ), with an average mass ratio of 0.5.The distribution of the companion masses as a function of the orbital semimajor axis (Fig. 5) is in good agreement with the population synthesis model; the model not only confirms that at separations up to 7 au (1500 R ) the companions have masses smaller than 6 M , it also Table 2. Orbital parameters for Cepheid systems that have spectroscopic orbital parameters.The lower part of the table lists the Cepheids with orbital periods shorter than 1000 days for which the orbital parameters and companion mass are poorly constrained (see Sect. 3.2.1).A null value for the eccentricity e and the argument of periastron ω indicates that the spectroscopic orbit was assumed to be circular.predicts more diverse companion masses (with the upper limit of 10 M ) for larger separations.Within the synthetic population, 99% of all companions to CCs are MS stars, and only 1% of the companions are giant stars.Among MS companions 75% have spectral types B and A, which introduce a non-negligible photometric contribution, particularly at short wavelengths.On average, the difference in magnitude between the binary system (CC with companion) and the CC alone is as large as 0.375 mag in the U-band, and decreases with increasing wavelength: 0.127 mag (B), 0.053 mag (V), 0.037 mag (R), 0.028 mag (I), 0.022 mag (J), 0.020 mag (K).The difference is larger for more massive MS companions and, naturally, for companions with larger physical radii (giant stars), which in extreme cases can contribute as much as 50% to the total luminosity of the system (Karczmarek et al., in prep.).
When radial velocities are not available, we estimate upper limits to the semimajor axis and orbital periods of the Cepheid systems (see Table A .4).These parameters can be used to determine the feasibility of the search for companions using classical imaging, adaptive optics, or interferometry.

PM anomalies of RR Lyrae
Table A.5 presents the results of our search for PMa in our sample of RRLs.Out of our list of 8509 RRLs, only 189 are present in the Hipparcos catalog and therefore suitable for a PMa analysis.These 189 stars gave five detections with a high ∆ > 5, and 8 stars with 3 < ∆ < 5 giving a minimum binary fraction for RRLs of P = 7% (the known non-RRL stars were removed from this count).In addition, 61 stars show suspected level PM anomalies at 2 < ∆ < 3. Figure 6 shows the statistics of the sample of RRLs with detected anomalies as a function of parallax.We do not include in this plot the RRLs with resolved candidate companions (Paper II).As was true for CCs, the fraction of detected PM anomalies decreases with the parallax, due to the decreasing sensitivity of the search technique.For the nearest targets, the binary fraction is approximately 0.4, which is probably a reasonable approximation of the true mean binary fraction of our RRL sample.
TU UMa is the only RRL that has been convincingly shown to be a member of a binary system (Szeidl et al. 1986;Kiss et al. 1995;Wade et al. 1999;Liška et al. 2016a).We observe a strong GDR2 PMa at ∆ G2 = 6.1 and also a significant Hipparcos PMa ∆ Hip = 2.8 that confirms the presence of an orbiting companion.Adopting the spectroscopic orbital parameters determined by Liška et al. (2016a) allows us to conduct a level 2 analysis and determine its complete orbital parameters.The result is presented in Table 3.The mass that we obtain for the companion (m 2 = 1.98 ± 0.33 M ) is high, and is due to the high inclination of the retrograde orbit of the system.This may imply that the companion of TU UMa is a massive white dwarf (a hypothesis already proposed by Kiss et al. 1995 andLiška et al. 2016a) or possibly a neutron star.Considering the old age of the RRL, a white dwarf companion will be cool and difficult to detect by imaging or interferometry, particularly as the angular separation with the primary is only on the order of 10 mas.It is important to note, however, that the PM vector from Hipparcos is imprecise for this star, with uncertainties larger than 1 mas a −1 on both axes.The orbital parameters will improve with the future Gaia data releases.
Estimates of the upper limits of the semimajor axes and orbital periods of a selected sample of RRLs are presented in Table A.6.Most of the detected RRL binaries are likely on very long-period orbits, with the exceptions of AT And, CZ Lac, and AR Ser.
Table A.7 presents the results of the PMa analysis for the stars that were incorrectly classified as CCs or RRLs.

Conclusion
We detected a significant number of new candidate companions of CCs and RRLs from the signature of their orbital motion on the proper motion vector of the targets.CCs have long been Table 3. Parameters of the TU UMa system from the combined analysis of the spectroscopic orbit of Liška et al. (2016a) (their Model 2) and the proper motion anomaly vectors.1.98 ± 0.33 M known to have a high binary fraction, and our survey of nearby CCs indicates that their binary fraction is likely above 80%; in addition, there is a significant fraction of triple or quadruple systems (e.g., Polaris, U Aql, W Sgr, AW Per, δ Cep).

Adopted parameters
The very small number of known binaries in the RRL class has long been a puzzle, but we detect significant PM anomalies ∆ > 3 for 7% of the 189 nearby RRLs that we surveyed, indicating that they are likely binaries.The massive companion of TU UMa, likely a white dwarf, points at the possibility that a significant fraction of RRL companions may be compact objects, which complicates their detection.
The presence of PM anomalies is an efficient way to determine the binary status of a large number of stars, and provides a valuable constraint on population synthesis models.The most interesting candidates can easily be identified for further characterization.In future Gaia data releases, the availability of timeresolved dynamical PM and radial velocity measurements opens the possibility, together with the parallax, of determining 3D linear velocity vectors of the photocenters of a massive number of binary and multiple systems.This will allow us to improve our understanding of their physical properties and of the role of binarity in stellar evolution.The combination of the future Gaia data releases with targeted spectroscopic observing campaigns will enable the thorough characterization of a large number of binary and multiple stars of all types with astrophysically interesting properties, at a reasonable cost in observing time.This will deeply improve our understanding of their physical properties and of the role of binarity in stellar evolution.is the GDR2 parallax (Sect.2.3), except when marked with * (Hipparcos; van Leeuwen 2007) and + (Berdnikov et al. 2000, for RY Vel only).The PM vectors at the Hipparcos (µ Hip ) and GDR2 (µ G2 ) epochs are compared to the mean PM computed using the Hipparcos and GDR2 positions (µ HG ).The observed differences are listed in terms of signal-to-noise ratio ∆ Hip and ∆ G2 .When the Arenou et al. (2018) quality parameters are not all satisfied, the star is marked with † after ∆ G2 and when the RUWE > 1.4 they are marked with ‡.The minimum ∆ anomaly is set to S/N=5 for a strong detection ( ), S/N=3 for a detection (•) and S/N=2 for a suspected binary (•).The binary type as listed in the database by Szabados (2003b)  Notes.The binary type is indicated using this code: B -spectroscopic binary (: when confirmation is needed); b -photometric companion, physical relation should be investigated; O -spectroscopic binary with known orbit; V -visual binary.

Fig. 1 .
Fig. 1.Principle of the search for a proper motion anomaly.µ Hip designates the Hipparcos proper motion vector (epoch 1991.25),µ HG the mean proper motion vector between the Hipparcos and Gaia DR2 positions, and µ G2 is the Gaia DR2 proper motion vector (epoch 2015.5).

Fig. 2 .
Fig.2.Orbits of V1334 Cyg A (orange ellipse) and its companion B (light blue) around their common center of mass fromGallenne et al. (2018b).The virtual orbit of the photocenter of the system is shown as a gray ellipse.The measured tangential velocity vector (proper motion anomaly) is represented at the Hipparcos and Gaia epochs.

Fig. 3 .Fig. 4 .
Fig. 3. Histogram of the detected PMa signal-to-noise ratios ∆ of CCs (right panel) and fraction of known binaries with respect to the total number of Cepheids per bin (left panel).The dashed red lines mark our binary detection threshold of ∆ = 3.

PFig. 6 .
Fig. 6.Left: Histogram of the RR Lyrae with detected proper motion anomalies (∆ > 3) as a function of parallax.Right: Binary fraction as a function of parallax.The error bars represent the binomial proportion 68% confidence interval.

Table A .
1. Proper motion anomalies of Galactic Cepheids.
is provided in the "Bin.type" column.The table is available at CDS.

Table A .
2. Continued from Table A.1.