Synthetic population of binary Cepheids. I. The effect of metallicity and initial parameter distribution on characteristics of Cepheids' companions

The majority of classical Cepheids are binary stars, yet the contribution of companions' light to the total brightness of the system has been assumed negligible and lacked a thorough, quantitative evaluation. We present an extensive study of synthetic populations of binary Cepheids, that aims to characterize Cepheids' companions (e.g. masses, evolutionary and spectral types), quantify their contribution to the brightness and color of Cepheid binaries, and assess the relevance of input parameters on the results. We introduce a collection of synthetic populations, which vary in metal content, initial parameter distribution, location of the instability strip edges, and star formation history. Our synthetic populations are free from the selection bias, while the percentage of Cepheid binaries is controlled by the binarity parameter. We successfully reproduce recent theoretical and empirical results: the percentage of binary Cepheids with main sequence (MS) companions, the contrast-mass ratio relation for binary Cepheids with MS companions, the appearance of binary Cepheids with giant evolved companions as outlier data points above the period-luminosity relation. Moreover, we present the first estimation of the percentage of binary Cepheids in the Large Magellanic Cloud and announce the quantification of the effect of binarity on the slope and zero-point of multiband period-luminosity relations, which will be reported in the next paper of this series.


INTRODUCTION
Classical Cepheids (hereafter referred to as Cepheids) are among the most famous and widely used cosmic distance calibrators; their high luminosity and characteristic light curves make them easily recognizable, while the period-luminosity relation (PLR) that they follow is considered universal and of superior accuracy to that offered by other types of radial pulsators. Still, Cepheids' PLR suffers from several systematic uncertainties, related to, e.g. binarity, metallicity, number of crossing of the instability strip, reddening, that hinder achieving a subpercent precision in distance determination.
Binary Cepheids, in particular, can be a potent source of systematic errors since they constitute 60 − 80% (and likely even more) of all Galactic Cepheids. This high Cepheid-binary fraction is supported by both theoretical (Neilson et al. 2015;Mor et al. 2017) and empirical studies (e.g. Szabados 2003;Kervella et al. 2019a). Both approaches have their limitations; theoretical results are strongly dependent on the input parameters, while empirical ones suffer from the selection bias. Indeed, an ultraviolet (UV) survey of binary Cepheids with hot mainsequence (MS) companions in the Milky Way (MW) was limited to stars with V < 8 mag (Evans 1992), leaving fainter binaries and binaries with cooler companions undetected. In the Large Magellanic Cloud (LMC), 25 spectroscopic binaries with Cepheids have been reported so far (Pilecki et al. 2021;Szabados & Nehéz 2012, and references therein), five of them being eclipsing binaries (Pilecki et al. 2018). In the Small Magellanic Cloud (SMC), only nine Cepheid binaries have been reported so far, among which two are spectroscopic binaries, another two show eclipsing variations, and the remaining five are firm candidates for Cepheid-Cepheid binaries (Szabados & Nehéz 2012, and references therein). Such scarcity of binary Cepheids in the LMC and SMC relative to the MW indicates a strong selection bias in the Magellanic Clouds, which favors Cepheids with giant (and possibly pulsating) companions and highly inclined orbits.
While resolved binary Cepheids are remarkable tools to determine geometrical distances and companions' dynamical masses with an astonishing 1% accuracy (Pietrzyński et al. 2010;Gallenne et al. 2018), unresolved binary Cepheids can bias the measurements in a number of undesired ways. For example, the presence of a companion can affect the radial-velocity curve of a Cepheid, impeding an accurate radius determination (Gieren et al. 1998). Spectroscopic analysis can yield inaccurate stellar parameters and abundances if single-star models are fitted to the combined spectrum of unresolved binaries (El-Badry et al. 2018). Astrometric solutions and parallaxes provided by the Gaia space mission are yet to be corrected for the variability of binary and pulsating stars; until this happens, Gaia parallaxes for binary Cepheids should be inferred from resolved companion stars (Kervella et al. 2019a).
Furthermore, contribution of a companion's light to the total brightness of the system causes Cepheids with unresolved companions to seem brighter than their single counterparts; this effect is largest in the near-infrared domain if the companion is a red giant (RG Pilecki et al. 2018), and in the UV domain if the companion is a hot MS star (Evans 1992). As a result, Cepheids with unresolved companions can alter the slope and the zero-point of the PLR, which has been predicted and described in a qualitative way (e.g. Szabados & Klagyivik 2012) but still lacks quantification. Binary Cepheids are expected to have different color indices, especially if the binary components have very dissimilar effective temperatures, which leads to an incorrect estimation of the reddening values toward the system. Luminous companions diminish the observed pulsation amplitudes of Cepheids (Pilecki et al. 2021), and may be in part responsible for the scatter in the pulsation period-amplitude relations (Klagyivik & Szabados 2009). Observed luminosities and amplitudes of unresolved binary Cepheids are unreliable points of reference for theoretical models of stellar pulsations and evolution, which operate within a framework of single stars. Moreover, Evans et al. (2005) reported that at least 44% of all binary Cepheids are in fact triple systems and Dinnbier et al. (2022) suggested that around half of all Cepheids form in triple and quadruple systems, instead of binaries. This information is crucial for mass determination of (assumed) binary components, as the presence of a third component leads to inaccurate estimations.
In order to address some of the aforementioned issues, we employ a theoretical approach called binary population synthesis. This method relies on approximate formulas that govern the evolution of single and binary stars, which makes it fast and efficient. Results of population synthesis are independent of reddening, blending, and selection bias, and therefore provide invaluable insight in the characteristics of binary Cepheids. For example, Mor et al. (2017) estimated that 68% of all Cepheids should have companions, and Neilson et al. (2015) reported that 35% of MW Cepheids should be detectable as spectroscopic binaries. Anderson & Riess (2018) assessed the impact of the photometric bias from MW Cepheids in wide binaries (a > 400 au) and open clusters on the value of the Hubble constant. While Anderson et al. (2016b) did not resort to the population synthesis to assess the excess light of a Cepheid binary with a MS star as a function of a Cepheid's pulsation period (log P − ∆M ), they recognized that population synthesis is required in order to characterize this relation more thoroughly.
In this paper, we present the most extensive study of synthetic populations of binary Cepheids up to date. We take full advantage of the population synthesis method that treats metallicity and binarity percentage as free parameters, which can be set to any arbitrary value or a grid of values. We create synthetic populations of three metallicities Z = 0.004, 0.008, 0.02, reflecting the metal content of classical Cepheids in the SMC, LMC, and MW, respectively, and with binarity percentages of 0%, 25%, 50%, 75%, and 100%. Such an approach gives us full control over the binarity and metallicity, and allows us to study the impact of these parameters on the observed characteristics of resolved and unresolved binary Cepheids.
Moreover, for every combination of metallicity and binarity, we test four sets of initial conditions and their effect on the outcome. We examine the entire collection of synthetic populations of binary Cepheids for similarities and differences between the variants, and compare our theoretical predictions with features observed in binary Cepheids in the MW and LMC. Following papers in this series will focus on the detailed analysis of PLRs and the quantification of the expected shift of their zeropoints due to binarity (Paper II), and the quantification of the shift in Cepheids' color indices, caused by companions' dissimilar effective temperatures, which affects the reddening toward binary Cepheids, and presents an opportunity to detect companions on the color-color diagram (Paper III).

SYNTHETIC POPULATIONS
The StarTrack population synthesis code (Belczynski et al. 2002(Belczynski et al. , 2008 is based on the revised formulae from Hurley et al. (2000Hurley et al. ( , 2002, fitted to detailed singlestar models with convective core overshooting across the entire Hertzsprung-Russell diagram (HRD), created by Pols et al. (1998). A number of enhancements implemented to the StarTrack code account for wind accretion through the Bondi-Hoyle mechanism, atmospheric Roche lobe overflow (Ritter 1988) and wind Roche lobe overflow (Mohamed & Podsiadlowski 2012;Abate et al. 2013).
We used StarTrack to generate populations of 200 000 binaries on the zero-age main sequence (ZAMS) of three metallicities Z = 0.004, 0.008, 0.02, and helium abundances Y = 0.248, 0.256, 0.280, reflecting the environmental properties of young stellar populations in the SMC, LMC, and MW, respectively. This metallicity comes from adopting the equation  . Throughout our study, we refer to the primary component (A) as the more massive star on the ZAMS while the secondary (B) is the less massive one. The maximum evolutionary age considered for each binary is 14 Gyr. In order to create a synthetic population that consists purely of binaries with Cepheids, we created a filtering algorithm to test if either component in a (synthetic) binary has stellar parameters characteristic of Cepheids, described below. Only binaries with at least one component that passed the test were included in the final sample.
First, we selected binaries with a component that crosses the instability strip (IS). Figure 1 shows two examples of stars that evolve through a Cepheid stage while inside the IS. The first crossing (IS1) happens when a star traverses IS in the Hertzsprung Gap, the second crossing (IS2) happens at the stage of core helium burning (blue loop) when a star traverses the IS toward the blue edge, and the third crossing (IS3) happens when it makes a blue loop toward the red edge. If a star makes its blue loop and turns around while still . Examples of evolutionary tracks of primary components in the log T −log L plane. Companions were omitted for clarity. Solid, dashed, and dotted lines indicate evolutionary stages before, during, and after core helium burning, respectively. Asterisk indicates the moment of helium ignition in the core. The instability strip loci and shapes are adopted from Anderson et al. (2016a).
inside the IS, the turning point divides its evolutionary stage into IS2 and IS3. We excluded all systems that experienced substantial mass transfer (MT), i.e. the mass lost or gained due to the MT constituted more than 10% of the star's initial mass 1 . More than 10% of mass lost/gained before the IS crossing would disrupt the physical structure of binary stars, raising a question of whether a Cepheid variable would retain its pulsation properties after the MT, and if so, whether it be justified to label it a Cepheid once its internal structure and evolutionary status has changed (Karczmarek et al. 2017). However, the mass loss due to stellar winds is common among Cepheids and ranges from 10 −10 to 10 −6 M /yr (Deasy 1988;Matthews et al. 2016); we included every star that experienced mass loss due to stellar winds below this upper limit. Another possible event is a supernova explosion, which however may disrupt the orbit of a binary. Since we focus on Cepheids in stable and uneventful binaries, we did not follow the evolution of disrupted binary components, meaning that all such stars were excluded from our sample.
Every binary selected so far hosts a star that crosses the IS at some point during its evolution, but it can only be regarded as a Cepheid if its birth time stamp equals its age in a Cepheid stage. Therefore, the last step was to assign birth time stamps to our synthetic binaries in order to determine their age and therefore the evolutionary status of the components. This procedure is described in detail in Section 2.3.
Constructing a synthetic population of binary Cepheids is a complex task with an ambiguous out-1 We chose a value of 10% to match criteria set by Neilson et al. (2015), so that their and our results can be compared.
come, because the output depends strongly on the input parameters. The results can be considered reliable if they agree with the observations or if they do not change much upon uncertainties introduced by the input physics relevant for the formation of binary Cepheids. We selected three areas that we consider are of most relevance to the reliability of the results: (i) distribution of initial parameters; (ii) shape and location of the IS; (iii) star formation history (SFH) that impacts Cepheids' birth time stamps. In the rest of this section, we detail the alterations introduced in these three areas, and we describe the process of augmentation of our synthetic data with pulsation periods, amplitudes, and multiband magnitudes, which allow for a comparison with observed binary Cepheids, and provide a frame of reference for future discoveries.

Initial distributions
A synthetic population is generated by StarTrack based on four initial parameter distributions: mass of the primary, mass ratio (secondary to primary), orbital separation (semi-major axis), and eccentricity. By default these parameters are independent and for every binary are drawn from the following distributions: • broken power-law initial mass function (Kroupa & Weidner 2003) with the value of a slope −2.35 for the mass of the primary M A from 2.5 to 12.0 M 2 , • flat distribution of mass ratio of secondary to primary q = M B /M A (Kobulnicky & Fryer 2007) in a range from q min to 1, where q min is the mass ratio that results in a secondary with M B = 0.08 M , • flat distribution of the logarithm of semi-major axis of binary orbit (Abt 1983) in a range from a min to −10 5 R , where a min is twice the sum of component's radii at periastron, • thermal eccentricity distribution f (e) = 2e (Heggie 1975) in range from 0 to 0.99.
The above set of initial parameter distributions is called "set A". Alternative distributions were published by Duquennoy & Mayor (1991), and recently by Moe & Di Stefano (2017). Duquennoy & Mayor (1991) reported normal distribution of mass ratios N (4.8, 2.3 2 ) , log-normal distribution of orbital periods log N (0.23, 0.42 2 ) , and a mixture of normal 2 The upper limit was chosen based on evolutionary models of massive Cepheids, which above 11 − 12 M fail to present a blue loop or their blue loop is erratic. The lower limit was chosen based on evolutionary models of low-mass stars, which below 2.5 − 3.0 M fail to cross the IS as post main-sequence objects.
N (0.27, 0.13 2 ) (for 10 d < P orb < 1000 d), thermal (for P orb 1000 d), and uniform (for P orb 10 d) distributions of eccentricities, with M A , P orb , and q independent, and e dependent on P orb . They did not state the distribution of M A , which we decided to keep as in set A. Their initial distributions are used in our work as "set B". Moe & Di Stefano (2017) reported e and q distributions that follow power laws with different exponent values, depending on P orb , and M A , while keeping M A and P orb independent. They did not describe distributions for P orb and M A . Thus, for the distribution of M A , we used the one from set A, while for P orb we used two already presented variants: log-uniform (as in set A) and log-normal (as in set B), creating two more sets, C and D, respectively. Table 1 summarizes the four sets of initial distributions, and Figure A in Appendix A supplements Table 1 with triangle plots of initial parameters for all four sets.
Binary populations created from the four sets differ especially in the distributions of orbital periods and masses of the companions. In Section 4, we show how the choice of a set impacts the distributions of orbital periods and effective temperatures of Cepheids' companions.

Shape and location of the instability strip
All Cepheid binaries were extracted from synthetic populations if they met a basic criterion: an evolved component (either primary or secondary, or both) had an effective temperature and luminosity that placed it inside the IS. We adopted two variants of the IS, shown in Figure 2: (i) simplistic, metallicity-independent parallel IS from Jeffery & Saio (2016, their Figure 1); (ii) metallicity-dependent and wedge-shaped IS from Anderson et al. (2016a). Although we did not set explicit upper and lower luminosity limits, the upper and lower mass limits imposed limits on luminosity via the mass-luminosity relation (see Section 3.1). As a result, the luminosities of the Cepheid components are 2 < log(L/L ) < 4.5.
The shape and location of the parallel IS is historically driven. The pioneering numerical investigations of the Cepheid IS explored only the blue edge, because of the insufficient computational power and knowledge about the convective processes inside Cepheid variables that impact the location of the red edge (e.g. Baker & Kippenhahn 1962). As a result, the red edge was estimated by assuming an ad hoc efficiency of convective transport or by shifting the blue edge redward by a fixed temperature, and thus was parallel to the blue edge. The wedgelike shape of the IS was first reported by Fernie (1990) based on empirical data of about 100 classical Cepheids, Note-Orbital periods and semi-major axes are interchangeable. This means that when one of them was drawn from a distribution, the other one was calculated. and later supported by a number of numerical studies, which included the effect of convection (e.g. Bono et al. 2000b;Anderson et al. 2016a). In this work, the choice of the IS plays an important role in determining whether a star of a given effective temperature is classified as a Cepheid. For instance, a MW star having T eff = 3.8 K and log(L/L ) = 3.5 will be found inside the parallel IS but outside the Anderson IS. The location of Cepheids on the HRD also impacts their pulsation periods and brightnesses, and therefore affects the PLRs; this effect is detailed in Section 3.

Star Formation History
The SFH is crucial to determine how many stars are observable as Cepheids now, given that they were born at a specific time in the past and have evolved to the point of crossing the IS. In this sense, the distribution of Cepheids' birth time stamps is a function of lookback time, with the current time (now) being 0. For instance, if a star that becomes a Cepheid at the age of 200 Myr is assigned a birth time stamp of 100 Myr, that means that at the time of observations (now) the star is only 100 Myr old and has not become a Cepheid yet. Consequently, such a star is observed as a nonpulsator. On the other hand, if the birth time stamp assigned to the same star was 200 Myr, this star would be observed as a pulsator, because it would cross the IS at exactly 200 Myr old.
Within the frame of population synthesis, we assigned to every binary a birth time stamp drawn from one of two variants of SFH: (i) uniform; (ii) based on Cepheids' ages. The uniform SFH means that all birth time stamps are equally probable for all metallicity environments 3 . The SFH based on Cepheids' ages was calculated from a period-age relation of Bono et al. (2005) using the pulsation periods of fundamental-mode Cepheids observed in the MW (Skowron et al. 2019) and the Magellanic Clouds (Soszyński et al. 2015). As a result of this calculation, three age distributions for the SMC, LMC, and MW were created, with age bins of width of 10 Myr. From these distributions, we drew time stamps of birth with the probabilities related to the heights of the bins. If the time stamp of birth (and therefore age) of a system was between the time of entering and exiting the IS for either binary component, then such a system was marked as a Cepheid binary and added to the final sample. Cepheid's location in the IS (closer to the edge or closer to the middle) and all related stellar parameters (i.a. T eff , L, R, P pul ) were linearly interpolated on the basis of its age relative to its time of IS entrance and exit. For instance, if a star's age is close to the time of IS entrance/exit, its location on the HRD is closer to the IS edge, and if a star's age is closer to the mean of IS entrance and exit times, it resides in the middle of the IS. This selection of stars based on their birth time stamps was repeated until the final sample reached 10,000 systems.
Our samples are much larger than the real samples of Cepheids; the Magellanic Clouds have approximately 5000 Cepheids each with completeness of virtually 100% (Soszyński et al. 2015), while 3352 MW Cepheids reported so far constitute a somewhat incomplete sample (completeness of about 88% down to a magnitude G = 18, Pietrukowicz et al. 2021). By keeping our synthetic samples this large, we allow binary Cepheids with more exotic and less probable companions to occur, and by keeping them equal in size, we assure that statistical errors affecting the sample size remain the same.
Our attempt to create a realistic mixture of Cepheids of various ages from the period-age relation was only partially successful, meaning that our synthetic Cepheids in the SMC, LMC, and MW are no older than 200, 170, and 100 Myr, respectively, while the Cepheids' ages calculated from the period-age relation (Bono et al. 2005) are as old as 350, 260, and 200 Myr for the SMC, LMC and MW, respectively. One of the reasons might be that we used the period-age relations of Bono et al. (2005) who assumed nonrotating progenitors on the MS.
3 Although popular, the uniform SFH has been recently challenged by e.g. Olejak et al. (2020) who presented that synthetic SFH in the MW depends not only on metallicity but also on Galactic component (disk, bulge, and halo). Their SFH at lookback time 0−5 Gyr, i.e. the age span of classical Cepheids, remains uniform.
On the contrary, rotating progenitors can experience enhanced internal mixing and therefore spend, on the MS, even twice as long as their nonrotating counterparts before they evolve into Cepheids (Anderson et al. 2016a). Consequently, Cepheids evolved from rotating progenitors are older. For the consistency sake, we calculated ages for our populations, which consist of nonrotating stars, using formulas of Bono et al. (2005). The other reason might be the fact that the theoretical models of stellar evolution fail to render extensive blue loops for low-mass (and therefore older) stars, and as a result, the sample consists of young Cepheids (either massive IS2+3 crossers or IS1 crossers). We provide more comment on this caveat in Section 3.2. Nevertheless, we proceed with our study bearing in mind that our results are relevant only for young and massive Cepheids and should be interpreted with caution in the context of older and low-mass ones.
The end result of all the above computations is 16 variants of synthetic populations of binary Cepheids for each of the three metallicity environments (SMC, LMC, MW), which were created as permutations of four variants of initial parameters distributions, two variants of the IS, and two variants of the SFH (4 × 2 × 2 = 16).

Magnitudes, amplitudes, and pulsation periods
In order to calculate multiband photometry, we used the online YBC database 4 of stellar bolometric corrections (Chen et al. 2019). We chose ATLAS9 model atmospheres of Castelli & Kurucz (2003) and derived UB-VRIJHK magnitudes in the Bessell & Brett photometric system (Bessell 1990;Bessell et al. 1998). We also created a reddening-free quasi-magnitude Wesenheit index, following the formula from Udalski et al. (1999): (1) V -band maximum peak-to-peak amplitudes were estimated based on the log P -A plot of Klagyivik & Szabados (2009, their Figure 1). We traced the envelope of highest V -band amplitudes as a function of log P , including the dip at log(P/d) ≈ 1. Next, we used theoretical predictions of Bhardwaj et al. (2017) for amplitudes in U , B, V , I, J and K bands in order to calculate ratios of amplitudes in these bands relative to the V -band averaged amplitude. The ratio A H /A V was interpolated linearly between A J /A V and A K /A V . The results of these calculations are as follows: As the last step, we multiplied the envelope of highest Vband amplitudes as a function of log P by the above constants. The amplitudes in the Wesenheit index were created from V , I magnitudes and amplitudes, following Eq. 1. Our estimated maximum amplitudes do not depend on metallicity.
Because StarTrack was not tailored to check whether a star is dynamically unstable and prone to pulsations, we assumed that all stars found in the IS pulsate as fundamental-mode Cepheids and calculated the pulsation periods using two external and independent sets of formulas. The first one was taken from Bono et al. (2000b): The other set of formulas for fundamental periods was created for the purpose of this study using Warsaw Pulsational Code (Smolec & Moskalik 2008) with the homogeneous envelope and convective parameters α = 1.5, From the grids, we extracted models on the blue and red edge of the IS, and using the information about their fundamental-mode pulsation periods, we created the period-luminosity-temperature relations, separately for the IS1 Cepheids:  (16) Bono's and our prescriptions for fundamental periods yield similar results, which agree best for short-period Cepheids and diverge for long-period ones (P pul 40 d) with Bono's periods being systematically longer by 4 − 6 d. Such long-period Cepheids are however rare and their contribution to the sample is minuscule.

SANITY CHECK
We performed three sanity checks in order to evaluate the agreement of selected parameters of our synthetic populations with the observed and literature data. These tests were executed on single Cepheids only, so that the reliability of our sample is endorsed before we introduce the next level of complexity to the analysis, i.e. Cepheids' companions.

Mass-luminosity relation
Mass-luminosity (ML) relation for Cepheids, presented in a form log L = α log M + β, depends on i.a. the metallicity, helium content, rotation, and/or overshooting on the MS, as well as the mass-loss rate (Chiosi et al. 1993;Alibert et al. 1999;Bono et al. 1999Bono et al. , 2000aSzabó et al. 2007;Anderson et al. 2014). In general, high rotation rate, large overshooting, and low metallicity cause the parameter β to increase, making such Cepheids more luminous by as much as log L = 0.25. In Figure 3 the ML relation for our exemplary synthetic population (set D, Anderson prescription for the IS, the SFH based on Cepheids' ages) is compared with ML relations from the literature, showing fair agreement. Noticeably, synthetic Cepheids in IS1 and IS2+3 obey different ML relations, in a sense that IS2+3 Cepheids are systematically brighter. These different ML relations have been taken into account while calculating pulsation periods using Warsaw Pulsational Code, as described in Section 2.4, and have resulted in slightly different PLRs for IS1 and IS2+3 Cepheids.

Proportions of Cepheids in the first, second, and third crossing
Since different SFHs favor different time stamps of birth, the choice of the SFH impacts the number of Cepheids on different IS crossings. The shape of the IS also impacts the number of Cepheids on different IS crossings, because it affects the time that stars spend inside the IS, but this effect is much smaller. Finally, the metallicity correlates with the sizes of blue loops and, consequently, the proportions of Cepheids on different IS crossings. The combination of all three factors results in different percentages of IS1, IS2 and IS3 Cepheids, shown in Figure 4. Not only primaries (blue areas in the plot) but also secondaries (red areas) can become Cepheids, although this scenario is much rarer. In such cases, the companions to secondary Cepheids are more-evolved asymptotic giant branch (AGB) stars or compact objects: white dwarfs (WD) or neutron stars (NS), and they are described in more detail in Section 4.2.
Cepheids on their first crossing, i.e. traversing the Hertzsprung Gap, are short-lived and therefore expected to be extremely scarce in the observational data. In our synthetic populations, the percentage of the IS1 Cepheids can vary from negligible to prevalent, depending on the choice of SFH, IS, and metallicity. In some variants (e.g. a MW population generated from parallel IS and uniform SFH), IS1 Cepheids-despite being short-lived-are much more abundant than IS2+3 Cepheids, and therefore could be observed more frequently.
However, this phenomenon arises from the fact that the theoretical models of stellar evolution for metal-rich stars of masses below ∼ 3.5 M fail to render extensive blue loops. Consequently, such stars cross the IS only once in the Hertzsprung Gap, which actually creates a deficiency of IS2+3 crossers. This issue has been approached by many authors (e.g. Xu & Li 2004, and references therein), but remains unresolved despite the evidence for the existence of short-period IS2+3 Cepheids (Turner et al. 2006;Rodríguez-Segovia et al. 2022).
In the majority of variants, the percentage of IS2 Cepheids is 40 − 60% (except for the MW Cepheids generated from the uniform SFH, where IS1 Cepheids dominate the samples, and thus the IS2 Cepheids are only 15 − 30%). In general, variants with the paral-lel IS produce fewer IS2 Cepheids than variants with Anderson IS, meaning that slightly more IS2 Cepheids constitute the samples if the red edge of the IS is shifted toward lower temperatures; this observation can again be explained by the aforementioned blue loop issue of theoretical models.
First-, second-, and third-crossers can be distinguished based on their rates of period changes (negative for IS2, positive for IS1 and IS3, and larger for IS1 than for IS3), because their pulsation periods decrease as they cross the IS toward the blue edge, and increase as they cross the IS toward the red edge. Turner et al. (2006) measured rates of period change for 200 MW Cepheids and reported that only 33% were negative (belonged to IS2 Cepheids). Poleski (2008) estimated that only 15% from 655 analyzed LMC Cepheids show consistent period change, and from those ∼ 57% have negative period changes. Theoretical estimations of the percentage of IS2 Cepheids with the metallicity of Z = 0.02 yielded only 10 − 15% (Neilson et al. 2012;Miller et al. 2020), and increased to 40 − 45% only after a significant initial rotation was introduced (Miller et al. 2020). Recently the most comprehensive study of pulsation period change of LMC Cepheids (Rodríguez-Segovia et al. 2022) shows that among 1303 objects 43% are IS2 crossers, 53% are IS3 crossers, and the remaining 4% are objects with inconclusive period changes, and two candidates for IS1 crossers. Our results agree with Poleski (2008) and Rodríguez-Segovia et al. (2022) on the percentage of IS2 crossers, but we find IS3 crossers underrepresented, and IS1 crossers overrepresented with respect to the results of Rodríguez-Segovia et al. (2022).

Multiband Period-Luminosity relations
Having calculated pulsation periods and magnitudes, we created multiband PLRs for all variants of our synthetic populations. Figure 5 illustrates the results for the variant with initial parameters from set D, Anderson IS, SFH based on Cepheids' ages, and our prescription for the pulsation periods. For a clearer comparison between the three metallicities, the absolute magnitudes are provided instead of observed ones.
The scatter of PLRs reflects the fact that Cepheids populate the entire width of the instability strip. Madore & Freedman (2012)  Chiosi (1993) Alibert (1999) Bono (1999) Bono (2000) Szabó (2007) Anderson (2014) 1st crossing 2nd crossing 3rd crossing We also notice that, in variants with parallel IS, the scatter remains constant for all values of log(P/d), but it grows with larger log(P/d) in variants with Anderson IS; this effect is especially visible in the B band, where the scatter is the largest in general. The varying scatter as a function of log(P/d) reflects the wedge-like shape of the Anderson IS. Empirical PLRs for short wavelengths (e.g. Musella et al. 1997;Bhardwaj et al. 2016) do not show larger scatter for larger log(P/d), which either supports the parallel variant of the IS over Anderson's wedge-like variant or suggests that not enough Cepheids have been observed to populate the PLR on the long-period end. Synthetic Cepheids cluster at log(P/d) ≈ 0.5, 0.6, 1.0 for the SMC, LMC, and MW, respectively, which is the minimum pulsation period for IS2+3 Cepheids, while all stars with shorter periods are IS1 Cepheids. Observational data support the existence of IS2+3 Cepheids with shorter periods (Turner et al. 2006;Poleski 2008), but theoretical models fail to render extensive blue loops for such stars (see Xu & Li 2004, and discussion in Section 3.2 of this paper).
IS1 Cepheids show a slightly steeper slope for their PLR than IS2+3 Cepheids, which was expected since we used different period formulas for IS1 (Eqs. 11-13) and IS2+3 Cepheids (Eqs. 14-16). This difference in slopes is best visible in the NIR passbands and W V I , and could be potentially used to distinguish between IS1 and IS2+3 Cepheids and reveal the number of IS1 Cepheids in the population, provided a large sample with accurate magnitudes. Figure 6 presents the comparison between a selection of Cepheid PLR slope values taken from the literature (listed in Table 2) and our synthetic Cepheid populations, for the B, V , I, J, H, and K bands, and the W V I Wesenheit index. The slopes of the synthetic Cepheid populations were calculated by the linear least-squares fit to joint samples of IS1+2+3 Cepheids, as seen in Figure 5. We used two different prescriptions for calculating pulsation periods: Bono's (Eqs. 2-4) and ours (Eqs. 11-16), which result in two different sets of slopes, marked in Figure 6 as circles and triangles, respectively. Four variants of synthetic populations are marked with different colors. These are all combinations of two IS prescriptions (Anderson's and parallel) and two SFH formulas (uniform and based on the ages of Cepheids). We found that the initial parameters (sets A-D) have no impact on the slope for single Cepheids, and are therefore omitted.
Our slopes (triangles) tend to be slightly steeper than Bono's (circles), and appear to fit the values of empirical slopes better, especially for longer wavelengths. However, the large scatter of the slopes for both the observed and synthetic populations prevents us from excluding or approving any variant of the synthetic population.
Synthetic slopes for a given passband have different values in the three metallicity environments, but the differences get smaller with the longer wavelengths, similarly to the empirical slopes of PLRs reported by Breuval et al. (2021). This result partially validates the assumption that the slopes are metallicity-independent and, in practice, can be therefore fixed for the distance determinations to Cepheids in farther galaxies (e.g. Wielgórski et al. 2017;Gieren et al. 2018), as long as near-infrared PLRs are used.

Characteristics of binaries
The distribution of three binary characteristics, orbital period log(P/d), eccentricity e, and mass ratio q = M B /M A , are presented in Figures 7, 8, and 9, respectively. These figures show the results for the set A  Table 2.
of the initial parameters, while the results for the sets B-D are shown in Figures B1, B2, and B3 of Appendix B.
Distributions of log P are very similar to their original initial distributions, which means that the orbital periods did not change significantly during the binary evolution. An exception is a peak at log(P/d) ≈ 3 for IS2+3 Cepheids (navy color), which hints that some of binary Cepheids shortened their orbital periods prior to their blue loop. In contrast, IS1 Cepheids (abundant in variants with uniform SFH and color coded as gray), tend to spread uniformly across the entire available range of orbital periods. Such clustering of IS2+3 Cepheids at log(P/d) ≈ 3 was caused by tidal interactions on the red giant branch (RGB), which led to the shrinkage and circularization of their orbits 5 . Indeed, distributions of eccentricities of binary Cepheids peak at e = 0 even though the initial distributions (sets A, B, D) did not favor this value. Apart from the peak at e = 0, the eccentricity distributions remain similar to their initial values.
Comparison with a population synthesis study of binary Cepheids in the MW, carried out by Neilson et al. (2015), shows a satisfactory agreement for distributions of orbital periods but considerable discordance for distributions of eccentricities. Indeed, Neilson et al. (2015) did not amplify the tidal forces, causing the values of 5 For particularly large eccentricities (e > 0.8), binary components tend to rendezvous at a very close proximity at the periastron, which amplifies the tidal forces. For binary stars with large convective envelopes (e.g. RGB stars), tidal forces in the Star-Track code are amplified by a factor of 50, which was calibrated based on the orbital separations and eccentricities of binaries in the Hyades open cluster ((Belczynski 2022, private communicatio). As a result, a considerable fraction of IS2+3 Cepheids, i.e. after the RGB evolutionary phase, have their orbits shrunk and circularized.
eccentricities and orbital periods to remain virtually unchanged throughout binary evolution. Eccentricities of observed MW binary Cepheids (Evans et al. 2005) favor the synthetic population of Neilson et al. (2015), suggesting that the amplification factor of the tidal forces in our population synthesis is too large for the population of Galactic binary Cepheids. Analogously to the orbital periods and eccentricities, the mass ratios of binary IS2+3 Cepheids preserve similar distributions as initially injected to the StarTrack code. For sets A and B, an excess of systems with a mass ratio of 0.2 − 0.3 can be observed, being larger and more clustered for MW binary Cepheids, and less visible for the SMC. Within one metallicity environment, the four combinations of IS and SFH prescriptions do not show noticeable differences for IS2+3 Cepheids. On the other hand, IS1 Cepheids show distributions of mass ratios that tend to cluster around 0.5 − 0.6 for the SMC, 0.4 − 0.5 for the LMC, and 0.3 − 0.4 for the MW, regardless of the set of initial parameters. We caution the reader that the characteristics of binary IS1 Cepheids have not been studied in detail yet, and therefore our results need further investigation with more precise tools, like the Modules for Experiments in Stellar Astrophysics (Paxton et al. 2019).
The key message from the above analysis of physical and orbital parameters on binaries is that the output of the binary population synthesis method is strongly affected by the input values and distributions. Therefore, one needs to be cautious not to trust the result based on only one set of input parameters, but instead investigate different sets and their outputs, in order to reliably assess the impact of initial parameters on the results.  Figure 7. Distributions of orbital periods in 12 variants of synthetic populations, for set A of the initial parameters. IS2+3 Cepheids (navy blue) and IS1 Cepheids (gray) are presented separately. The values on all y-axes were scaled linearly from 0 to 1, and then omitted for a clearer comparison of the shapes of the distributions. Images for sets B, C, and D can be found in Appendix B. spectral type. They are presented in Figures 10, 11, and 12, respectively. These figures present the results for a set A of the initial parameters, while the figures for sets B, C, and D are specified in Appendix B.

Characteristics of companions
In most cases (70 − 90%), companions to Cepheids are MS stars in every set of initial parameters (Figure 10).  The highest percentage of MS companions is in variants with the SFH determined by Cepheids' ages, for all metallicity environments. Evolved companions are as follows: red giants or horizontal branch stars (cumulatively denoted as RG+HB), AGB stars, WDs, and NS are also possible. We report no Cepheid binaries with  . Mass ratios in 12 variants of synthetic populations, for set A of the initial parameters. IS2+3 Cepheids (navy blue) and IS1 Cepheids (gray) are presented separately. The values on all y-axes were scaled linearly from 0 to 1, and then omitted for a clearer comparison of the shapes of the distributions. Images for sets B, C, and D can be found in Appendix B. black holes. Such systems have disrupted their orbits in a supernova event, preceding the creation of a black hole, and as such were excluded from our sample at the stage of data filtering (see Section 2). Among evolved companions, WDs dominate in scenarios with the uniform SFH for all sets of initial parameters and all metallicity environments, while in scenarios with the SFH based on the observed-age distribution of Cepheids, all types of evolved companions are similarly numerous. Notably, RG+HB companions are more common in the sets A, C, and D (3 − 5%) but extremely scarce in the set B (0.8%).
The presented distributions support the idea that binary Cepheids are indeed common, but hard to observe using photometric methods, since the majority of companions are MS stars and their contribution to the overall luminosity of the system is minuscule. In the case where IS1 Cepheids are the primary components, their companions are MS stars by default, because they are less massive and therefore less evolved. Cepheids as secondary components can have companions at any evolutionary stage, but they constitute a marginal fraction of the sample (see Figure 4).
Cepheids' companions have diverse spectral types, and no spectral type is strongly favored over others. Stronger preference for early-type companions (O, B) is visible in sets A and B, but it is only moderate in sets C and D. On the other hand, sets C and D tend to favor K-type companions more than sets A and B. Minuscule differences in spectral types for different metallicities and combinations of IS and SFH variants suggest that the companions' spectral types are solely correlated with initial parameters, among which the most important is the initial mass ratio distribution. Indeed, in sets C and D, the initial mass ratio distribution shows a peak at q = 1, which favors components of a similar mass and therefore a similar evolutionary stage, increasing the probability of late-type companions. Figure 11 shows that 20−40% of MW Cepheids have a companion of spectral type B, which agrees with Evans (1992) who reported that at least 20% of all binary Cepheids in the MW have a companion of a spectral type earlier than A0V.
Surface temperatures of companion stars, although coarsely encoded by their spectral types, are worthy to be investigated on their own. We found that log T eff ranges from 3.5 to 4.5, but in the case of the SMC and LMC (top and middle panel), it tends to cluster around 3.6, 3.8, and 4.2 for B, C, and D sets, respectively (the equivalent spectral types are K5, F6, and B4, respectively Figure 10. Evolutionary stages of companions in 12 variants of synthetic populations, for set A of the initial parameters. Images for sets B, C, and D can be found in Appendix B.
IS1 Cepheids tend to present more erratic distributions, rarely coinciding with that of the IS2+3 group.

Binaries with two Cepheids
In a handful of cases, both components cross the IS simultaneously, and create a Cepheid-Cepheid (C-C) binary. Only two 6 such systems have been reported so far:   Figure 11. Spectral types of companions in 12 variants of synthetic populations, for set A of the initial parameters. Images for sets B, C, and D can be found in Appendix B. Indefinite spectral types belong to Cepheids' companions that evolved into neutron stars. AB Berdnikov 1990;Kervella et al. 2019b) and OGLE-LMC-CEP-1718 (Soszynski et al. 2008;Gieren et al. 2014;Pilecki et al. 2018), suggesting that C-C binaries are extremely rare. Our synthetic populations support this conclusion, as presented in Table 3. In almost all cases, C-C binaries constitute less than 1% of all Cepheid binaries but usually are much more scarce. Values in parentheses show the median duration of the simultaneous Cepheid phase. This phase can occur for pairs of Cepheids on the same IS crossing (e.g. IS2+2, IS3+3) or two different ones (e.g. IS2+3, IS1+3  and short-lasting; thus the majority of C-C binaries are second-or third-crossers.

CE Cassiopeiae (CE Cas
Masses of IS2+2, IS2+3, or IS3+3 pairs are virtually the same; the relative differences in mass are smaller than about 2%. Relative differences in radii are smaller than 25%, and relative differences in log(L/L ) are up to 5%. These values mark upper limits for expected relative differences of masses, radii, and luminosities in C-C binaries in the LMC. Indeed, OGLE-LMC-CEP-1718 has a relative mass difference of 1.2%, relative radius difference of 19%, and relative log(L/L ) difference of 4.6% (Pilecki et al. 2018). CE Cas, which belongs to the MW open cluster NGC 7790 (Berdnikov 1990), is a visual C-C binary on an extremely wide orbit; its projected separation of log(a/R ) ≈ 6.176 corresponds to the orbital period of about 5000 yr (Kervella et al. 2019b). A large 2.3 arcsecond angular separation of the components enables it to collect photometric and spectroscopic data for each star and investigate their properties individually (Berdnikov 1990;Majaess et al. 2013). Unfortunately, this large separation precludes from collecting astrometric and/or velocimetric data, leaving this C-C binary impossible to compare with the synthetic results on the level that was possible for OGLE-LMC-CEP-1718.
Careful inspection of percentages and median lifespans of C-C binaries in Table 3 shows that in the majority of cases C-C binaries at metallicities characteristic to the SMC live longer and are slightly more abundant than C-C binaries at metallicities characteristic to the LMC and MW. One could therefore expect to observe more C-C binaries in the SMC than in the LMC and MW; however, with just one C-C binary confirmed to date in each galaxy (OGLE-LMC-CEP-1718 and CE Cas), it is impossible to favor/disfavor any particular variant of synthetic populations, based on this prediction alone.

The estimated binarity fraction of LMC Cepheids
Different detection methods allow for the discovery of binary Cepheids of different physical and orbital characteristics, but also have their own limitations. For example, binary Cepheids discovered due to eclipses in their light curves have companions of similar size, orbital inclination close to i = 90 • , and relatively short orbital periods, so that the eclipses can be observed over a finite time span. Within 29 yr of the Optical Gravitational Lensing Experiment project, five LMC eclipsing binaries, consisting of a classical Cepheid and an evolved companion (RG or another Cepheid) have been discovered (Soszynski et al. 2008;Pilecki et al. 2018). They share similarly high values of inclination, i 83 • , and short orbital periods, log(P/d) < 4. By comparing these observational data and our synthetic populations, we make a first estimate of the number of binary Cepheids in the LMC.
Let us assume that six binary Cepheids (in five systems, because one consists of two Cepheids) set a lower limit for the number of eclipsing binary Cepheids that can be observed given the detection conditions de-  , then six detected binary Cepheids with i 83 • constitute 8% of all 75 binary Cepheids, out of which 92% are undetectable due to their unfavorable inclination angles i < 83 • . We assume that binary Cepheids with longer orbital periods cannot be detected due to an insufficient time base of observations, regardless of their inclination angles. Next, we recall that binary Cepheids with RG+HB companions constitute 3−5% of the entire population for sets A, C, D and 0.8% for set B (Section 4.2 and Figures  10 and B3). We extract the fraction of these systems with log(P/d) < 4 and we make this value equal to 75 systems derived from the previous step. Now we can calculate the number of binary Cepheids with RG+HB companions over the entire log P range and compare it with the percentage of RG+HB companions in the whole population of binary Cepheids. This value, divided by the number of LMC Cepheids (4620, Soszyński et al. 2015), yields the binarity fraction of classical Cepheids in the LMC. Depending on the variant of our synthetic population, we get the lower limit for a binarity fraction that ranges from 55% to beyond 100% (Table 4). In particular, variants from set B produce nonphysically high values (above 700%), further confirming that set B should be disregarded altogether. Our crude estimation shows that the binarity fraction of classical Cepheids in the LMC should be at least 55% and likely much higher.

Mass ratios of MW Cepheid binaries with MS companions
In the MW, binary Cepheids with MS companions can be detected with interferometric observations; however, such detections are limited to large orbital separations ( 100 milliarcseconds) and relatively bright companions, i.e. the difference in the magnitude be-  tween the companion and the Cepheid (the so called contrast) must be smaller than 6 mag in H band (Gallenne et al. 2018(Gallenne et al. , 2019. In order to determine the physical properties of the binary components (especially masses), interferometric observations have to be supplemented with radial-velocity measurements, which is a time-consuming task, because the orbital period of a binary with a Cepheid is at least one year long (Neilson et al. 2015;Pilecki et al. 2018). Moreover, in the majority of cases, spectroscopic observations are performed in the visual domain, where hot MS companions are too faint to be detected. For such companions, more challenging UV spectroscopy has to be carried out from space (Gallenne et al. 2018). Our synthetic populations offer a unique insight into the relations between physical parameters of binary components and their observed properties. In particular, Figure 13 presents mass ratio as a function of contrast in the V and H bands for the two most divergent variants: set D, Anderson IS, SFH based on Cepheids' ages; and set A, parallel IS, uniform SFH. Complementing figures for B, I, J, and K bands are available in Appendix C. For comparison purposes, we overplot gray curves showing the contrast-mass ratio relation derived and averaged over second and third IS crossings by Anderson & Riess (2018, their Figure 2). We extend our analysis to IS1 crossers, and compare with observational data, summarized in Table 5. Data points with red bor-ders (V1334 Cyg, AX Cir, AW Per, U Aql) are binary Cepheids with contrasts in both the V and H bands, and are expected to occupy the same region of either IS1 or IS2+3 crossers on all plots in Figure 13. Indeed, V1334 Cyg and AX Cir lie on the IS2+3 trend, but AW Per and U Aql lie below it, in the area of IS1 crossers. While our result does not prove their evolutionary status, it entertains an uncommon idea that IS1 Cepheids might not be as rare as previously thought. Further investigation of the rates of pulsation period change would be required to determine their evolutionary status.
The contrast-mass ratio relation depends on the massluminosity relation, which was established by numerous authors based on models of stellar evolution (see Figure 3 and Section 3.1 for reference). Such models tend to overestimate Cepheid masses by 10 − 20% relative to pulsation models (the so called mass discrepancy problem, Keller 2008). The data point of V1334 Cyg, which was derived from the dynamical mass of the Cepheid and its companion (Gallenne et al. 2018), fits perfectly in the IS2+3 trends in both filters. However, other data points, for which the Cepheid masses were estimated from the evolutionary models via mass-luminosity and period-luminosity relations (Evans 1995), might be underestimated, meaning that they should be shifted upwards in the plot to fit the IS2+3 trend.  Table 5. Systems present on both panels are additionally marked in red. The underlying synthetic populations correspond to the following parameters: set D, Anderson IS, SFH based on Cepheids' ages (top); set A, parallel IS, uniform SFH (bottom). The rest of populations fit between these two most divergent variants. Complementing panels for the B, I, J, and K bands are presented in Figure C1.
Contrast between a Cepheid and its companion strongly depends on the pulsation phase of the Cepheid and is bigger when the Cepheid is brighter. If the Cepheid is observed at a random phase, an additional statistical error should be added to the error budget of the contrast value. V -band contrasts in Table 5 were calculated from Cepheid magnitudes averaged over their pulsation cycles. H-band contrasts were either calculated from random-phase observation or averaged over a couple of observations, and therefore their errors might be underestimated.
The above relations offer a promising opportunity to estimate the mass ratios using just a single interferometric image, instead of a number of radial-velocity measurements. This method could be applied to very wide binaries, for which spectroscopic observations could take many years to complete.

Detection of Cepheid Binaries above the Period-Luminosity Relation
Following the hypothesis of Pilecki et al. (2021), that Cepheids located well above the PLR are binaries, we reproduced their Figure 3 using our synthetic population of LMC Cepheids (set D, IS: Anderson, SFH: Cepheids' ages) with a binarity fraction of 100%, meaning that every Cepheid has a companion. The total brightness of a binary was calculated as follows: m tot = −2.5 log 10 −m A /2.5 + 10 −m B /2.5 .
Peak-to-peak amplitude of a binary was calculated as a difference between the minimum and maximum brightness of a binary, A tot = m tot, min − m tot, max . Next, the relative pulsation amplitude with respect to the original pulsation amplitude of a single Cepheid was calculated as 100% × A tot /A Cep . A relative amplitude of 100% means that the extra light from the companion was negligible, and the Cepheid retains its full amplitude, while a relative amplitude of 0% means that the Cepheid amplitude was completely diminished by the extra light from the companion. Figure 14 presents our synthetic population on the period-luminosity plane. Black lines mark thresholds for stars that are 50%, 100%, and 200% brighter than an average Cepheid of a given period (gray line). All outliers (above 50% threshold) are binary Cepheids with giant or supergiant companions. Colors encode relative amplitudes; particularly interesting are magenta data points, representing binary Cepheids with RG+HB companions, whose light contribution places the Cepheids above the dashed line while preserving about 30-60% of their pulsation amplitudes. Such stars present an unmatched observational opportunity to find binaries among Cepheids.
A handful of binary Cepheids located above the dotted 200% line, marked as yellow or bright orange, have AGB companions. Such systems evolved with Cepheid progenitors as secondary components (less massive at ZAMS). AGB companions dominate the brightness of the system, and also heavily diminish the apparent pulsation amplitude of Cepheids. Consequently, the Cepheid is found well above the PLR, but its pulsation amplitude is so low that it might be mistaken for another type of small-amplitude pulsator or overlooked altogether and treated as a constant star.
MS and WD companions, on the other hand, contribute little or almost no light to their systems (navy blue data points) and lie very close to the main trend and well below the 50% threshold. However, they outnumber companions of other evolutionary stages and constitute the vast majority of Cepheid binaries. Because of their location on the period-luminosity plane, they are virtually undetectable via photometric methods, yet their cumulative effect on the zero-point and the slope of the PLR is not negligible, and will be thoroughly characterized in Paper II.

SUMMARY
We presented the synthetic populations of binary Cepheids for three environments of different metallicity: the SMC, LMC, and MW. For each metallicity, we crated 16 different variants of synthetic populations, testing two prescriptions for the shape of the IS, two prescriptions for the SFH, and four sets of initial parameter distributions. Our synthetic populations are free from selection bias, and the percentage of Cepheid binaries is controlled by us via the binarity parameter.
We compared all variants with the literature and concluded that the most realistic synthetic populations are the ones created from the IS prescription of Anderson et al. (2016a), and the SFH based on the Cepheid ages (Bono et al. 2005). We dissuade using set B of the initial parameters as it resulted in unrealistically high fractions of binary Cepheids in the LMC, and unrealistically low fractions of binary Cepheids with giant, evolved companions.
Hot MS stars constitute 20 − 40% of all companions, which agrees with the empirical study of Evans (1992). Such companions show narrow contrast-mass ratio relations, already suggested by Anderson & Riess (2018), and replicated by us. Comparison of our theoretical results with empirical values of mass ratios and V -, Hband contrasts shows satisfactory agreement, and encourages further investigation of the contrast-mass ratio relations as an efficient tool to estimate the mass ratios of binary Cepheids.
We reported that giant, evolved stars constitute 3 − 5% of all companions, and by comparing observational data with our synthetic populations, we estimated the number of binary Cepheids in the LMC, which is at least 50% and probably much higher (close to 100%). We confirmed that Cepheid binaries with giant companions can be easily detected as outliers above the PLR. MS companions lie well below the detection threshold in the period-luminosity plane, but their effect on the PLR is nonnegligible and will be the focus of Paper II.

ACKNOWLEDGMENTS
We thank the anonymous referee, whose pertinent comments helped to significantly improve this paper. P.K. thanks P.W. for insightful discussions on the topic of this study and beyond. The research leading to these results has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreements No. 695099 and No. 951549 Software: StarTrack (Belczynski et al. 2002(Belczynski et al. , 2008, SciPy (Virtanen et al. 2020), NumPy (Harris et al. 2020), Pandas (McKinney 2010), JuPyter (Pérez & Granger 2007). All images were proudly made with Matplotlib (Hunter 2007).

A. INITIAL DISTRIBUTIONS
We present triangle plots of initial distributions of sets A, B, C, and D. Differences with respect to metallicities are negligible. For sets C and D, we assumed, following Moe & Di Stefano (2017), that the distributions of initial primary masses and orbital periods are independent, and drew eccentricities and mass ratios from the distributions given in their Table 13. Distributions in sets B-D are based on observational data, and approximated by analytical functions. Such functions are of different types (e.g. power-law, log-normal, log-uniform) and/or have different parameter values (e.g. the exponent of a power-law function) in different ranges of one distribution (e.g. q, e). As a consequence, such distributions are not smooth, but remain continuous.

B. CHARACTERISTICS OF BINARIES AND CEPHEIDS' COMPANIONS
Figures 7, 8, and 9 presented three binary characteristics: orbital period log(P/d), eccentricity e, and mass ratio q = M B /M A , respectively, for set A of the initial parameters. This section complements the presented results with the plots for sets B, C, and D.

C. MULTIBAND CONTRAST-MASS RATIO RELATIONS FOR MW CEPHEIDS AND THEIR COMPANIONS
We present contrast-mass ratio relations for the B, I, J, and K bands, complementing the V -and H-band relations, presented in Figure 13. Note that empirical contrasts are available only in V and H band (Table 5), and therefore Figure C1 presents solely theoretical results.  Figure A1. Visualization of the four sets of initial distributions described in Table 1.  Figure B1. Distributions of orbital periods (top) and mass ratios (bottom) in 12 variants of synthetic populations, for sets B, C, D of initial parameters. IS2+3 Cepheids (navy blue) and IS1 Cepheids (gray) are presented separately. The values on all y-axes were scaled linearly from 0 to 1, and then omitted for a clearer comparison of the shapes of the distributions.  Figure B2. Distributions of eccentricities (top) and companion effective temperatures (bottom) in 12 variants of synthetic populations, for sets B, C, D of the initial parameters. IS2+3 Cepheids (navy blue) and IS1 Cepheids (gray) are presented separately. The values on all y-axes were scaled linearly from 0 to 1, and then omitted for a clearer comparison of the shapes of the distributions.   Figure C1. Magnitude difference (contrast) between MW Cepheids and their companions in the B, I, J, and K bands, versus the mass ratio from two populations variants: set D, Anderson IS, SFH based on Cepheids' ages (left); set A, parallel IS, uniform SFH (right). The rest of populations fit between these two most divergent variants. This figure complements Figure 13, which shows V and H-band contrast-mass ratio relations.