Y Gem: A White Dwarf Symbiotic Star?

In this work we conduct a thorough investigation of the X-ray and ultraviolet (UV) properties of Y Gem based on six archival XMM-Newton and Chandra observations to explore the nature of the system. The results show that Y Gem has strong (1032–34 erg s−1) X-ray emission, including a hard (with a maximum emission temperature of 8–16 keV) and a soft (with emission temperatures of 0.02–0.2 and 0.2–0.9 keV) component. The integrated UV luminosity of Y Gem reaches ∼1035 erg s−1. We show that the previous asymptotic giant branch-main-sequence (AGB-MS) Roche-lobe overflow (RLOF) scenario is dynamically unstable and can hardly explain the ∼10 keV X-ray emission temperature. We propose Y Gem as a symbiotic star, where a white dwarf (WD) accretes from its AGB companion based on its X-ray and UV properties. We make numerical simulations to examine the evolutionary history of this system. The simulations can produce the observed properties of Y Gem in the wind WRLOF scenario. An ∼0.8M ⊙ WD with a ∼1.0–1.8M ⊙ companion in a ∼2000–32,000 day initial orbit may evolve to a Y Gem-like system. Our finding implies a potential population of symbiotic stars that may have been misclassified as AGB-MS binaries. What is more, their high mass accretion rates may enable mass accumulation to the WD and makes them candidates of Type Ia supernovae progenitors.


Introduction
Asymptotic giant branch (AGB) stars are medium-mass (∼1-8M e ) stars in their late evolutionary stage. The helium supply in their cores is exhausted so that the helium burning only continues in the inner shells. Binaries containing AGBs are particularly interesting, because they are not only important laboratories to test the stellar evolution theory (e.g., the massloss history of giants and supergiants), but also are related to interesting astrophysical objects or phenomena such as planetary nebulae and Type Ia supernovae (SNe Ia).
The population and evolution of AGB binaries have not been adequately explored. For example, the binary fraction of AGB stars is still unknown (see Sahai et al. 2008, and reference therein), and the mass-loss rate of AGB stars due to the presence of a companion is highly uncertain (Tout & Eggleton 1988;Meng & Han 2016). The population of AGBs with compact companions in the Milky Way is also not known, with predicted numbers from 10 3 -10 5 (Mukai et al. 2016;Mukai 2017;Munari 2019). To improve our knowledge of these issues, a good sample of AGB binaries is crucial. However, the search for AGB binaries via photometry or radial velocity measurements is difficult because AGB stars show strong intrinsic variability, and are so luminous (10 3 -10 4 L e ) that they often dominate their main-sequence (MS) companions in the optical band (van Winckel et al. 1999;Bond 2000).
Ultraviolet (UV) observations have been incorporated to search for AGB binaries. The basic idea is that AGB stars have too low surface temperatures to emit significant UV radiation beyond 2800 Å, so the detection of strong UV emission from an AGB most likely indicates the existence of a hot companion star and/or accretion process (e.g., Ortiz & Guerrero 2016, also see Montez et al. 2017 for a different interpretation). In the past decade or so, systematical investigations on AGB binaries have been made based on UV observations with the Far-Ultraviolet Spectroscopic Explorer, the Galaxy Evolution Explorer (GALEX), and the Hubble Space Telescope, which have greatly increased the number of AGB binaries and candidates (e.g., Sahai et al. 2008Sahai et al. , 2011Ortiz & Guerrero 2016;Montez et al. 2017). For example, Sahai et al. (2008) detected significant GALEX FUV emission from nine AGBs. Afterwards, Sahai et al. (2015) showed that three out of six far-UV (FUV) AGBs have strong X-ray emission with luminosities of 10 30-33 erg s −1 and plasma temperatures ranging from (35-160)×10 6 K, which were proposed as evidence of mass accretion to the companion stars (Sahai et al. 2015(Sahai et al. , 2018. Recently, Ortiz & Guerrero (2016) found 34 AGBs within 500 pc of the Sun that exhibit FUV emission and/or a near-UV (NUV) observed-to-predicted flux ratio >20, which were proposed to be related to the presence of hot companions or accretion flows.
The nature of the companion stars in these AGB binaries/ candidates is an interesting, but still not clear astrophysical question. Either an MS or a compact companion can accrete from the AGB star, heat the accreted matter, and emit UV photons. Specifically, unlike the MS accreting case, if the accretor is a white dwarf (WD), strong X-ray emission with temperature above several kiloelectronvolts is expected, making the system a symbiotic star (WD SySt; e.g., Luna et al. 2013). WD SySts are rare and interesting objects because some of them may contain massive WDs that could be progenitors of SNe Ia (for a recent review of WD SySts, see Mukai 2017 and reference therein). The total number of WD SySts in the Milky Way is highly uncertain, ranging from 10,000 to a few million (Munari & Renzini 1992;Lü et al. 2006;Akras et al. 2019), which may partly be due to the bias brought by the traditional spectroscopic detection method (Mukai et al. 2016). Some of the UV-emitting AGB binaries, especially those with strong X-ray emission, are good candidates for WD SySts. Thus the investigation of these binaries/candidates is also important to constrain the population of WD SySts and potentially the progenitors of SNe Ia.
In this work we focus on Y Gem, one of the three X-ray emitting FUV AGBs. Y Gem is significantly different from the other FUV AGBs in several aspects. Its X-ray luminosity (8.5 × 10 32 erg s −1 ) is about two orders of magnitude higher than the others (Sahai et al. 2015), and its plasma temperature (∼10 keV) is more than twice that of the others (Sahai et al. 2015). The mass accretion rate of the companion star inferred from the combined UV and X-ray luminosities reaches  > M 5 × 10 −7 M e yr −1 (Sahai et al. 2011(Sahai et al. , 2018. Theoretically, some of these results can hardly be explained with the MS companion-Roche-lobe overflow (RLOF) AGB scenario. First, the ∼10 keV emission temperature is too high for an accreting MS star (see Section 4.1.1 for details). Second, the inferred accretion rate (>5 × 10 −7 M e yr −1 ) is so high that the accretion is most likely in the form of RLOF (Sahai et al. 2018). However, for an AGB star, the RLOF mass transfer to a 0.35M e MS companion is dynamically unstable (see Section 4.1.1 for details). Third, the photometric flickering found in Y Gem (Snaid et al. 2018) indicates disk accretion onto a WD instead of an MS star. These inconsistencies deserve further investigations to improve our understanding of the nature of Y Gem, and more importantly, shed light on the origin and evolution of similar systems.
In this work, we make a thorough investigation on six currently available X-ray observations of Y Gem (including four XMM-Newton and two Chandra ones, which triples the number of X-ray observations in previous works). Based on the results, we propose a WD-SySt scenario to explain the observed X-ray and UV emissions of Y Gem. This scenario can naturally explain the high emission temperature detected in X-ray and high UV luminosity. We also make numerical simulations to explore possible progenitors of Y Gem. Our results implies a potentially hidden population of WD SySts that may have been misclassified as AGB-MS binaries.
The rest of the paper is organized as follows. In Section 2, we present the observations and the spectral and timing analysis methods. In Section 3, we present the results. We present a brief discussion on the nature of the source and perform numerical simulations in Section 4 and summarize our main findings in Section 5. In this work, the distance to Y Gem is adopted as

Observations and Data Analysis
Y Gem has been observed four times by XMM-Newton and twice by Chandra, as listed in Table 1. One of the XMM-Newton observations (ObsID 720340201) and one of the Chandra observations (ObsID 15714) have been analyzed in previous works (Sahai et al. 2015). In this work we make a spectral and timing analysis for all six observations, and focus on the XMM-Newton observations because they provide better counting statistics and simultaneous UV observations.

Spectral Analysis
We start by reprocessing the XMM-Newton data using the emproc (MOS1, MOS2) and epproc (PN) tools in the Science Analysis System (SAS v16.1.0) software. We exclude the time intervals when the count rate above 10 keV is higher than a threshold specific for each observation (see the last column in Table 1). 3 On-source events are extracted from a 40″ circle centered at the source, and background events from a 2′ circle on the same chip, beside the source region.
We then perform spectral fitting for the 0.3-10 keV XMM-Newton MOS1, MOS2, and PN spectra for each observation and for the combined spectra. We find that the spectra could not be well fitted by a simple absorbed one-temperature (apec) or multi-temperature (mkcflow) thermal plasma model, with large residuals between the 0.3 and 2 keV range, which is similar to the findings by Sahai et al. (2015). We then use the following two models, phabs 1 (pcfabs (apec 1 +Gaussian))+phabs 2 (apec 2 +apec 3 ) and phabs 1 (pcfabs(mkcflow+Gaussian))+phabs 2 (apec 2 +apec 3 ) to fit the spectra following Luna et al. (2013). In these models, the apec 1 or mkcflow component describes the 2-10 keV (hard X-ray hereafter) emission with respective foreground (phabs 1 ) and possible intrinsic (pcfabs) absorptions. The apec 2 and apec 3 components mainly describe the 0.3-2 keV (soft X-ray)  Sahai et al. (2015). b Calculated from EPIC-PN data for XMM-Newton observations. c Count rate threshold of flaring particle background. emission absorbed by phabs 2 , and the Gaussian component represents the Fe I-Kα line centered at ∼6.4 keV. The fitting results are presented in Table 2 and the physical meaning of the main parameters listed in the table is described as follows: N H,1 and N H, pcfabs represent the equivalent absorption column densities of the foreground material and those that are intrinsic to the binary, respectively; CF pcfabs represents the covering fraction of the intrinsic absorption; N H,2 represents the absorption caused by the stellar wind; T apec,1 or T max represents the characteristic or maximum temperature of plasma in the boundary layer; T apec,2 and T apec,3 represent the additional temperatures of the radiating plasma produced by the collided wind. We also fit the 5-8 keV XMM-Newton spectra to constrain the equivalence widths (EWs) and the flux ratios (I 6.4 /I 6.7 and I 7.0 /I 6.7 ) of the iron K shell lines. The results are presented in Table 2. The fitting is performed with the model phabs(apec +threeGaussian), where the metal abundance of the apec model is set to zero and the threeGaussian model was specifically built to measure I 7.0 /I 6.7 by Xu et al. (2016).
For Chandra observations, we use the chandra_repro tool in CIAO 4.11 to reprocess the ACIS-S data, and extract the spectra with specextract. This extraction uses an on-source 5″ radius circle centered on the source, and a background region as an annulus of the inner and outer radii equal to 10″ and 50″. The 0.5-8 keV Chandra spectra are then fitted in the same way as the XMM-Newton spectra, and the results are presented in Table 2. Due to the poor counting statistics, however, the Fe line EWs and the I 7.0 /I 6.7 cannot be constrained, and thus will not be discussed further.

Timing Analysis
Various methods (e.g., Muno et al. 2003;Hong et al. 2012) have been suggested to perform X-ray timing analysis. In this work, we follow the method used by Hong et al. (2012), applying the Lomb-Scargle method to the backgroundsubtracted light curves to search for periodic signals. First, the background-subtracted light curves are generated for each observation, where on-source and background are extracted from the same region as the spectral analysis. For XMM-Newton observations, we utilize epiclccorr to obtain the background-subtracted light curves for PN with a bin size of 10 s, where photon arrival times are barycenter corrected by the tool barycen in SAS. For Chandra observations, the background-subtracted light curves are generated for time bins of 12.8 s, and the photon arrival times are corrected by the tool  Notes. Model 1: phabs 1 (pcfabs(apec 1 +Gaussian))+phabs 2 (apec 2 +apec 3 ), Model 2: phabs 1 (pcfabs(mkcflow+Gaussian))+phabs 2 (apec 2 +apec 3 ). Errors show 90% confidence level. No error means the parameter cannot be constrained by observations. a Upper limit of temperature is set to 2 keV. b Unabsorbed luminosity from Model 2. c T max derived from I 7.0 /I 6.7 measurements here and the -I I T axbary in CIAO. Then, we utilize the Lomb-Scargle routine provided by the AstroPy (Astropy Collaboration et al. 2018) to search for period signals, which are further examined with the epoch folding method (Leahy et al. 1983). The candidate frequencies are determined with the false alarm probabilities P < 0.01, where Baluev (2008)ʼs approximation is used.

UV Data Analysis
Simultaneous UV observations in UVW2 band 4 (effective wavelength λ = 212 nm) were obtained for all four XMM-Newton X-ray observations using the Optical Monitor (OM) onboard XMM-Newton. The data were downloaded and reprocessed with the standard pipeline omfchain. The UV light curves and the mean count rates (N UVW2 ) in the UVW2 energy bands are also extracted and presented in Figure 1 and Table 3. We further estimate the unabsorbed UV flux density in UVW2 band (F UVW2 ) from the UVW2 count rate, using a count rate conversion factor of 5.71 × 10 −15 erg s −1 cm −2 Å −1 /(counts s −1 ) 5 , accounting for the color excess (E(B −V) = 0.051, (Ortiz & Guerrero 2016) and adopting R V = 3.1 and A UV /A V = 3.0 (Gordon et al. 2003). Assuming that the UV spectra follow Planck's law, we can estimate the integrated UV luminosities from F UVW2 . To do so, we adopt 30,000or 17,000 K as the temperature, which was given by Sahai et al. (2011) when Y Gem was at the high or low state. Two UV luminosities are calculated using the two temperatures, as listed in Table 3. Figure 2 presents the best-fitted XMM-Newton and Chandra spectra, and the fitting results are presented in Table 2. Judged by the χ 2 values, both models can well describe the 0.3-10 keV XMM-Newton spectra. The best-fitted foreground absorption column densities for the hard component (n H1 = 2 − 13 × 10 22 cm −2 ) are comparable to previous measurements (∼8 × 10 22 cm −2 , Sahai et al. 2015), but are higher than the one (2.8 × 10 20 cm −2 ) derived from the color excess (E(B −V) = 0.051 (Ortiz & Guerrero 2016) adopting a conversion factor of n H = 1.79 × 10 21 × A V cm −2 and (Predehl & Schmitt 1995;Gordon et al. 2003). The high foreground absorption value seems quite common for SySts and is likely to be due to the extra X-ray absorption produced by the accretion inflow/outflow and the AGB stellar wind (Luna et al. 2013). The measured emission temperatures (T apec ∼5-8 keV and T max ∼8-16 keV) are consistent with those derived from the I 7.0 /I 6.7 values, and are similar to those of SySts and nonmagnetic CVs (non-mCVs) in the solar vicinity (Luna et al. 2013;Xu et al. 2016). The soft component shows an unabsorbed luminosity comparable to the hard one (3-10 × 10 32 erg s −1 ), with emission temperatures of 0.02-0.9 keV and a small foreground column density (below several 10 21 cm −2 ). These values are again consistent with the so-called β/δ-type SySts with WDs as accretors (Luna et al. 2013).

Spectral Analysis
On the other hand, the T max and T apec of the hard component of the Chandra spectra are all above 10 keV. Both are significantly higher than the XMM-Newton measured values.
What is more, the unabsorbed X-ray luminosity of Chandra ObsID 16683 (∼2 × 10 34 erg s −1 ) are one order of magnitude higher than all the other five observations, indicating that Y Gem may have experienced an outburst due to an increase of the accretion rate (see Section 4.1.2 for details) during that observation.

Timing Analysis
For the period search, we find four candidate periods (488, 3005, and 6158 s in XMM-Newton and 4283 s in Chandra observations) with false alarm probability below 1% assuming white noise conditions. 6 We notice that the 4283 s candidate period is similar to the previous ∼1.35 hr signal by eye fit (Sahai et al. 2015). However, none of them are detected in all observations (Figure 3), and they are all too short to be the orbital period of a binary system containing an AGB star. What is more, these candidate periods are unlikely the spin periods, which are not supposed to change significantly in timescales of several days (e.g., between the second and third XMM-Newton observations).
We conclude that the spin or orbital period of Y Gem is not found in our analysis.

UV Flux and X-Ray-to-UV Luminosity Ratio
As shown in Figure 1 and Table 3, Y Gem showed both short-term and long-term variations in the UV band, which is similar to the findings of Sahai et al. (2015). The mean UVW2 count rates varied within an order of magnitude, from 51.8 counts s −1 in ObsID 0720340201 to 7-8 counts s −1 in the other three observations. The high count rate in ObsID 0720340201 is again consistent with Sahai et al. (2015). The integrated UV luminosity of Y Gem ranges from several 10 33 to ∼10 35 erg s −1 , consistent with the estimation (∼5 × 10 34 erg s −1 ) from Hubble observations by Sahai et al. (2018). Such high integrated UV luminosity results in rather low X-ray-to-UV luminosity ratios, from ∼0.02 to 0.18 (Table 3).

Discussion
In all the six observations spanning over 2 yr, Y Gem showed strong UV (L UV ∼ 10 33 erg s −1 to ∼10 35 erg s −1 ) and X-ray emissions (L X ∼ 10 32 -10 34 erg s −1 ) and high emission temperature (T max ∼8-16 keV). Based on these results, we discuss the nature of Y Gem, and perform numerical simulations to search the possible progenitors in this section.  Sahai et al. (2015Sahai et al. ( , 2018 proposed a 0.35M e , Roche-lobe accreting MS star as the companion to explain the observed X-ray and UV properties of Y Gem. However, both the high X-ray emission temperature (T max ∼8-16 keV) and the inferred high accretion rate (> 5× 10 −7 M e yr −1 , Sahai et al. 2015) can hardly be explained by this proposal.
First, the gravitational potential of a 0.35M e MS star could not heat the accreted matter to ∼10 keV. As a rough estimation, the maximum temperature that the accreted matter can reach 4 UVM2 data with effective λ = 231 nm was obtained only for ObsID 0720340201, and therefore is not used in this work. 5 https://www.cosmos.esa.int/web/xmm-newton/sas-watchout-uvflux 6 If both the white and red noise are taken into consideration, the false alarm probability of the four candidate periods are all above 50% judged by the Bayesian test introduced by Vaughan (2010).
(adopting freefall assumption) would be where μ is the mean molecular weight of the accreted matter, m H is the mass of the H atom, k is the Boltzmann constant, G is the gravitational constant, M and R are the mass and radius of the accreting star, respectively. For a 0.35M e MS star, T max ∼ 1 keV, 7 which is about an order of magnitude lower than the measured value. Figure 1. X-ray (0.3-10 keV for XMM-Newton and 0.3-8 keV for Chandra data, respectively) and UV light curves (UVW2 band for XMM-Newton data) of Y Gem. Black and red points represent UV and X-ray light curves, respectively. X-ray data have been rescaled with factor 60, 30, 40, and 20 in the first four panels. 7 The temperature in the disk accretion case would be even lower, about half of this value.
Second, the RLOF mass transfer from an AGB star to a 0.35M e MS companion is dynamically unstable and lasts for a maximum of 10 4 yr before the systems enters the so-called common envelope stage where the MS spirals in the envelope of the AGB star (e.g., Taam & Sandquist 2000). As a result, the chance of finding such an object would be extremely low. To illustrate this, we make numerical simulations with the stellar evolution code MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2017(Paxton et al. , 2019. We start with a zero-age MS binary with M 1 = 0.35M e , and the mass of the companion is set to M 2 = 1 or 2M e (the mass transfer of more massive companions would be even more unstable). The metallicity of both stars are set to 0.4Z e to match the observation. The initial orbital period is set to 500, 1000, and 2000 days to allow the companion to evolve to the AGB phase. The M 1 = 0.35M e companion is assumed to accumulate all accreted matter. Without exception, the mass transfer in all the simulations becomes unstable and the system enters the common envelope phase after 10 3-4 yr of high rate (> 5 × 10 −7 M e yr −1 ) mass transfer. This time interval is only ∼10 −5 to 10 −7 of the evolution timescale of a 1-2M e star (∼10 9-10 yr), thus the chance of observing such a system is also low. So we conclude that the MS companion could not reproduce the observed properties of Y Gem.

WD-SySt Scenario
The WD companion scenario was also introduced in Sahai et al. (2015), where the cool WD accretes from the AGB star via either RLOF or wind capture, making Y Gem a WD SySt. However, this scenario was later rejected by Sahai et al. (2018) because the outflow velocity (∼500 km s −1 ) derived from the absorption line profiles is about an order of magnitude lower than the escape velocity of a WD (∼5000 km s −1 ).
However, we argue that the low outflow velocity may not be used as a definitive evidence to reject the WD scenario because WD SySts could also have low outflow velocities, although the exact origin of the outflow has not been well established. For example, Sion et al. (2002) found an outflow velocity of 170 km s −1 for RW Hydrae, which is a WD SySt. Moreover, the recent discovery of the photometric flickering in Y Gem also favors the existence of a WD companion (Snaid et al. 2018). We then discuss if the accreting WD scenario can account for various observed properties of Y Gem as follows.
First, the X-ray and the UV emissions of Y Gem can be naturally explained by the WD scenario. The deep gravitational potential of a WD can easily heat the accreted matter to tens of kiloelectronvolts, as found in accreting WDs in the solar vicinity (e.g., Luna et al. 2013;Yu et al. 2018). What is more, the WD scenario only requires an accretion rate on the order of ∼10 −9 to 10 −8 M e yr −1 to explain the UV and X-ray luminosity, as shown in Section 4.2. In fact, WD SySts in the solar vicinity do show strong UV (up to several ∼10 34 erg s −1 with variability on timescales of minutes to hours) and X-ray (∼10 31-33 erg s −1 ) emissions that are comparable to those of Y Gem (Luna et al. 2013;Mukai et al. 2016;Mukai 2017).
Second, the X-ray spectra of Y Gem resemble those of the β/δ-type WD-SySt binaries. First classified by Muerset et al. (1997) and extended by Luna et al. (2013), WD SySts can be divided into several subclasses: the α type, where all the X-ray photons are detected below 1 keV, which originates from quasi-steady shell burning on the surface of the WD; the β type where the soft (<2.4 keV) X-ray emission is from the collision of winds from the WD and the companion; the δ type where the highly absorbed hard X-ray emission is suggested to come from the boundary layer; and the β/δ type, where both the soft (from wind collision) and hard (from boundary layer) X-ray emissions are discovered. Showing both the soft and highly absorbed hard X-ray components, Y Gem is most likely a member of the β/δ-type WD-SySt binaries.
The δ component is supposed to be from the boundary layer of the WD, where the luminosity could be estimated by in the standard picture. It implies L X-ray ≈ L UV , assuming the X-ray and UV emissions are mainly from the boundary layer and the accretion disk, respectively, which was not the case in Y Gem. Similar discrepancies have been observed in diskaccreting dwarf novae (DNe), a subtype of non-mCVs (e.g., in VW Hyi and SS Cyg, where L BL /L disk < 0.1 Mauche et al. 1991Mauche et al. , 1995 and almost all the β/δ subclass of WD-SySt binaries (see Table 2 in Luna et al. 2013 for details). It is currently an unresolved problem and may be ascribed to the following scenarios: (i) the WD spins at high speed so the gravitational energy released will be much less than predicted (e.g., Mauche 2004); (ii) The high accretion rate leads to the presence of the optically thick component (T ∼ 15-25 eV) in the boundary layer, so the boundary layer emits most photons in the UV band rather than in the X-ray band (e.g., see the review by Balman 2020, and references therein). Finally, the previously observed change in UV and X-ray luminosities of Y Gem could be explained with the change of accretion rates caused by the disk thermal-viscous instability, the pulsation of the AGB star, or the possible eccentric orbit.  Note. Columns: ObsID-XMM-Newton observation IDs; N UVW2 -UVW2 mean count rate; F UVW2 -UVW2 flux density; L UVW2 -the luminosity for the UVW2 filter; L UV -the integrated UV luminosity estimated from F UVW2 assuming a blackbody temperature of 30,000/17,000 K (Sahai et al. 2011); L X /L UVW2 -X-ray-to-UVW2 luminosity ratio; L X /L UV -X-ray-to-UV luminosity ratio.
First, outbursts occur in disk-fed WD accretors (e.g., DNe) due to the thermal-viscous instability of the disk (see Hameury 2020, and references therein). Here we perform a tentative test to see whether outbursts can occur in Y Gem. The thermalviscous instability theory (Lasota et al. 2008) predicts the critical local accretion rate as yr −1 . Obviously, the accretion rate of Y Gem fits in this zone through its evolution (see Section 4.3 for details); thus, Y Gem would experience the disk instability. As a result, the accretion rate of the WD undergoes changes. Second, the pulsation of the AGB star would certainly alter the stellar wind and thus the mass transfer rate (see Figure 5 for details). Third, we could not rule out the possibility that the binary is in an eccentric orbit. In this case, the mass transfer rate will increase when the WD reaches the periastron and decrease when the WD reaches the apoastron. As a result, the accretion luminosity would also be varying from one observation to another. The above possibilities should be verified with further observations.

AGB-NS Scenario
Besides WD-SySt binaries, there is a special subclass of SySts, namely, γ-type SySts, where the accretor is a neutron star (NS). Some γ-type SySts have X-ray luminosities comparable to those observed in Y Gem (e.g., 4U 1700+24 and 4U 1954+31; Masetti et al. 2006;Enoto et al. 2014). However, their X-ray spectral properties are different from SySts with WD accretors. Their X-ray spectrum commonly exhibits a soft blackbody component and a hard Comptonized continuum in an NS accretor (Shapiro & Teukolsky 1983). The Comptonized component is a nonthermal emission mechanism and would not generate emission lines like Fe XXV-Kα or Fe XXVI-Kα, which is not consistent with the results of Y Gem. Thus, the AGB-NS scenario is not discussed further in this work.

Estimating the System Properties of Y Gem
Adopting the WD-SySt scenario, we can estimate the system properties of Y Gem based on its observational properties.
First, the L X /L UV values suggest that Y Gem is likely a diskaccreting system. For WDs with magnetic field strength B = 10 6 G, the mass accretion is usually in the form of a disk. The presence of a disk causes strong UV emission from the inner disk and/or the boundary layer, which results in small X-ray-to-UV luminosity ratios (L X /L UV < 0.1), as found in nonmagnetic CVs in the solar vicinity (Patterson & Raymond 1985;Giovannelli 2008). On the contrary, for WDs with strong enough magnetic fields (B > 10 6 G), the magnetic field controls the accretion flow and prevent the formation of the accretion disks. The accreted matter are instead channeled directly to the magnetic poles. As a result, these binaries usually have L X /L UV > 0.1, as found in magnetic CVs in the solar vicinity (Patterson & Raymond 1985;Giovannelli 2008). In three out of four XMM-Newton observations (except in ObsID 0763050401), Y Gem shows L X /L UV < 0.1, and thus implies the existence of an accretion disk in Y Gem. The disk scenario is also supported by the photometric flickering in Y Gem (Snaid et al. 2018).
Second, we try to constrain the (lower limit) mass of the WD from the X-ray spectra. For a disk-accreting WD, the mass of the WD (M WD ) determines its gravitational potential, and therefore the maximum shock temperature (T max ). The T max is commonly obtained from the spectral fitting of the X-ray continuum (e.g., Kimura et al. 2021), but it is difficult to constrain the WD mass without 10 keV X-ray spectra. Fortunately, the iron line flux ratio could be an indicator on T max , and therefore the WD mass. The higher temperature means more Fe XXV ions would be ionized into Fe XXVI ions, leading to the higher iron line flux ratio (I 7.0 /I 6.7 ). As a result, a positive correlation between M WD and I 7.0 /I 6.7 exists if the boundary layer is optically thin (e.g., Yu et al. 2018). On the other hand, if the boundary layer is optically thick, the boundary layer can still emit some residual optically thin X-ray emission where the measured T max is lower than in the optically thin case (Wheatley et al. 2003). Assuming Y Gem is a diskaccreting WD and follows their M WD -I 7.0 /I 6.7 correlation, the WD mass can be estimated as . This WD mass should be considered as the lower limit because Y Gem may maintain a (partially) thick boundary layer like some of the β/δ-type WD SySts (Luna et al. 2013) due to its high mass transfer rate.
Third, we estimate the mass accretion rate in Y Gem from its peak X-ray plus UV luminosities, assuming both of which are powered by accretion. 8 With a peak accretion luminosity of L UV + L X ≈ 1 × 10 35 erg s −1 , the mass accretion rate to the WD can be inferred as ∼1 × 10 −8 M e yr −1 assuming a 0.8M e accreting WD. Certainly, the required mass accretion rate decreases with the increasing WD mass, e.g., to 6.5 × 10 −9 and 3.8 × 10 −9 M e yr −1 for a 1 and 1.2M e WD, respectively. On the other hand, the real mass transfer rate of Y Gem may be higher than the above value, as what was found in long-period (P orb > 5 −6 hr) CVs in the solar vicinity (Baskill et al. 2005;Suleimanov et al. 2019). This phenomena could be due to the mass loss from the binary in the form of disk/WD wind (e.g., Drew & Proga 2000;Luna et al. 2013). The matter lost in the wind may then collide with the AGB stellar wind and generate the β-type 0.3-2 keV soft emission in addition to the δ-type 2-10 keV hard X-ray emission (Luna et al. 2013). Additionally, the UV and X-ray luminosities imply no quasi-steady shell burning on the surface of the WD, which places a constraint to the accretion rate of the WD. We adopt the critical mass accretion rate above which the quasi-steady shell burning would initiate from that in Shen & Bildsten (2007), and obtain the upper limit of the allowed mass accretion rate as ∼10 −7 M e yr −1 by assuming a 0.8M e WD with the observed metallicity.

Simulation
To evaluate the reliability of the WD-SySt scenario, we perform numerical simulations with MESA to see if it can reproduce the observational results. The Dutch wind scheme is adopted in the simulations, and magnetic braking is assumed to only work for stars less massive than 1.4M e which have convective envelopes (Verbunt & Zwaan 1981;Rappaport et al. 1983;Lin et al. 2011;Kalomeni et al. 2016;Shao & Li 2020). The metallicity is set to Z = 0.4Z e to match the observation. All other parameters are set to the default values in the code.
The binaries are assumed to contain a 0.8M e WD, which accretes from the companion via either stellar wind or RLOF. The initial companion masses are set from 0.8-2.4M e with a step of 0.2M e ; less massive stars are not included since the time required for them to evolve to AGB stars would be longer than the age of the universe. The initial orbital periods are set from 10 3.0 -10 5.0 days with an interval of 0.1 in logarithmic space.
To compare with observations, we adopt the present system parameters as follows: L b = 5800L e (Sahai et al. 2011) and T eff = 2800 K (Sahai et al. 2011) for the AGB companion 9 , and E(B − V ) = 0.051 (Ortiz & Guerrero 2016). An extra 20% uncertainty is assumed for L b and T eff , since no error ranges were provided in previous works. As required by the observed X-ray and UV luminosities, only systems with mass transfer rates of 1 × 10 −8 M e yr −1   -M M 1 10 7   yr −1 are considered as candidates.
In the simulation, we find that neither RLOF nor wind accreting binaries could fulfill the observables of Y Gem. The RLOF systems with L b = 5800L e would experience unstable mass transfer which leads to common envelope evolution. On the other hand, the default Dutch wind loss rate adopted in the code are too low (with maximum wind-accretion rate ∼10 −9 M e yr −1 ) to power the UV and X-ray emissions (required   -M M 10 8  yr −1 ). Thus, a new wind-accretion scenario is required to power the observed X-ray and UV luminosities. Mohamed & Podsiadlowski (2007) introduced a physical mechanism named the wind Roche-lobe overflow (WRLOF), which leads higher mass transfer rates than the Bondi-Hoyle mechanism (Bondi & Hoyle 1944) in SySts with long orbital periods and have been applied in many works (e.g., Iłkiewicz et al. 2019;Belloni et al. 2020). We then adopt the wind-accretion efficiency The soft X-ray component was suggested to have originated from the shocked wind (Luna et al. 2013), and the inclusion of this component only makes a <2% difference of the peak accretion luminosity in ObsID 0720340201, and thus does not affect the validity of the following discussions. 9 We do not adopt the surface gravity ( g log ) value obtained by Ortiz & Guerrero (2016) because it is obviously too high for a 2800K AGB star, as shown in the same paper (Ortiz & Guerrero 2016).
given by Abate et al. (2013), where q, R d , and R RL represent the mass ratio, the dust formation radius, and the effective Rochelobe radius, respectively. In the above equation, R d follows as suggested by Höfner (2007), where R * , T eff , and T cond represent the radius and the effective temperature of the companion and the condensation temperature of the dust near the companion, respectively. We then perform simulations again and show the allowed parameter spaces in Figure 4, where the allowed systems are marked as rectangles and systems shown with triangles evolve to Y Gem-like binaries only when the companions are in the thermal pulse AGB (TPAGB) phase. Other systems that could not evolve to Y Gem-like systems are shown with crosses with various colors: not fulfilling  M (purple), too high L b (orange), not enough time to evolve in Hubble time (red), and the luminosity and the effective temperature are consistent with the observed only when the RLOF is present during the AGB phase (cyan). We also show the end products of the companions with different background colors in Figure 4: CO WD (blue) and common envelope evolution (purple). We note that some systems in the WRLOF scenario will undergo RLOF mass transfer within 10 −4 to 10 −2 M e yr −1 , during which the companions expand to large radii in AGB phases, but durations are only several decades and they cannot evolve to the common envelope stage. We find the WRLOF scenario can reproduce the observed properties of Y Gem, where the initial companion masses are between 1.0 and 1.8M e and the initial orbital periods are between 10 3.3 and 10 4.5 days. When the progenitors evolve to Y Gem-like systems, they are in ∼2000-220,000 day orbits, and their companions are ∼0.61-1.77M e AGB stars, whose bolometric luminosities are ∼4640-6960L e and effective temperatures are ∼3080-3360 K. The WD is now accreting at 10 −8 to ∼3 × 10 −8 M e yr −1 , enough to power the observed X-ray and UV emissions. The WD will continue accreting with a rate of 10 −8 to ∼3 × 10 −8 M e yr −1 for ∼5 × 10 4 to ∼7 × 10 5 yr (typically ∼4 × 10 5 yr). We present an evolutionary path of Y Gem-like in the WRLOF scenario as an example in Figure 5. The simulation starts with a 0.8M e WD with a 1.2M e companion in an ∼3162 day orbit. After ∼ 4.846 × 10 9 yr, the companion evolves to a TPAGB star with L b ∼ 4700L e , T ∼ 3300 K, and M ∼ 1.18M e . After several thermal pulses, the binary evolves to the stage at ∼4.847 × 10 9 yr, where its observables best match those of Y Gem: L b ∼ 6200L e , T eff ∼ 3200 K,  ´-M M 2.4 10 8 yr −1 , log(g) ∼ −0.31. During the Y Gem-like stage, the windaccretion rate is between ∼ 1.1 × 10 −8 and ∼3 × 10 −8 M e yr −1 (which is well below the critical mass accretion rate leading to the quasi-steady shell burning), the corresponding bolometric luminosity is between ∼4640 and ∼6960L e , and the effective temperature is between ∼3110 and ∼3320 K. The total time that the system is in the Y Gem-like stage is ∼ 6 × 10 5 yr. Finally, the companion evolves to a ∼0.67M e WD, making the system a binary WD in a ∼4720 day orbit at ∼4.848 × 10 9 yr.
The above simulations certainly suffer from various uncertainties and limitations. First, the 0.8M e WD mass is derived from the maximum emission temperature, which represents the temperature in the optically thick boundary layer case and should be more or less lower than the value in the optically thin case. As a result, this WD mass should be considered as the lower limit of the real WD mass. In reality, more massive WDs would stabilize the mass transfer and result in larger allowed parameter spaces. Second, we do not include the evolution of the progenitor of the primary WD in the simulation, which may affect the evolutionary path. For example, the progenitor star of the current WD could have evolved to late stage and have lost mass to the current AGB star and caused an orbital shrinkage, which should be considered in future works. Third, we have not taken the possible orbital eccentricity into consideration in the simulations, which may have a minor influence on the results. Lastly, the effective temperatures of the AGB star in simulations (∼3080-3360 K) are within 20% of the assumed uncertainty range, but seems a little bit higher than the best observed value (2800 K). Better matches could be reached for higher metallicity values (e.g., for Z = Z e ), since higher metallicities will lead to larger stellar envelopes, and thus lower effective temperatures.
Accounting for the above uncertainties and limitations, we propose further photometric and spectroscopic observations to put tighter constraints on Y Gem's bolometric luminosity, effective temperature, and orbital period. Such results would be essential for better measurements on the masses and the orbital separation of Y Gem, and to distinguish its accretion form (e.g., disk-like or accretion column). Further X-ray and UV observations would also be useful to follow the mass transfer process in Y Gem.

Potential Hidden Population of Symbiotic Stars
The finding that Y Gem could be a WD-AGB binary implies a potentially undiscovered population of WD SySts and related astrophysical objects. For example, there could be more UVemitting WD-SySt binaries that were misclassified as AGB-MS binaries due to the lack of X-ray observations and/or the absence of strong emission lines (Mukai et al. 2016;Munari et al. 2021). As shown in this work, the main differences between these WD-SySts and AGB-MS binaries are that the former shows higher X-ray emission temperature and stronger UV and X-ray emissions. Thus, the most effective way to find WD SySts would be X-ray surveys such as eROSITA. Another possible method would be looking for giants and AGBs with high UV luminosities. To illustrate this point, we tentatively cross correlate the GALEX All-sky Imaging Survey GR6 + 7 catalog (Bianchi et al. 2017) and the Gaia DR2 catalog (Lindegren et al. 2018) with a 3″ matching radius. The matching results in 20 giants/AGB stars with FUV plus NUV luminosities comparable to Y Gem (>10 34 erg s −1 , see also Montez et al. 2017 for a full table of UV-emitting AGB stars detected by GALEX), which are good candidates of Y Gemlike WD SySts and should be discussed in further works. With possible contaminating sources, the number of these candidates is already higher than the 16 confirmed SySts within 500 pc (Akras et al. 2019), and may represent the hidden population of WD SySts that have been missed in low resolution spectroscopic surveys as suggested by Mukai et al. (2016). The confirmation/rejection of these binaries to be WD SySts would certainly place tighter constraints on the population of SySts, e.g., if there are 10 3 or 10 5 WD SySts in the Milky Way galaxy (Munari & Renzini 1992;Lü et al. 2006).
We strongly encourage further spectroscopic and X-ray (e.g., eROSITA) surveys on more FUV-emitting giant/AGB stars to search for Y Gem-like systems. Combined with Gaia data, a volume-limited sample of such binaries would be reachable, which would greatly improve our understanding of the population of WD-SySt binaries and progenitors of SNe Ia.

Summary
In this work we investigate the nature of Y Gem based on six archival XMM-Newton and Chandra observations. We propose that Y Gem is a symbiotic star, where a WD is accreting from its AGB companion. We also search for possible evolutionary paths of Y Gem by conducting numerical simulations. Our main findings are summarized as follows: 1. Y Gem appeared to be a highly variable X-ray source, with 0.3-10 keV luminosity in the range of 5 × 10 32 -2 × 10 34 erg s −1 . The system showed both hard and soft X-ray components. The former can be reasonably well characterized by either a single-temperature (T apec ∼ 5 −8 keV) or a multitemperature (T max ∼ 8-16 keV) optically thin thermal plasma model. The latter can be described by two optically thin thermal plasma model with T apec ∼ 0.02-0.2 keV and T apec ∼ 0.2-0.9 keV. The UV luminosity is ∼10 34 −35 erg s −1 . The derived accretion rate is on the order of ∼10 −9 to ∼10 −8 M e yr −1 assuming a 0.8M e WD as an accretor. The timing analysis shows no convincing periodic signal.
2. We propose that Y Gem is a WD-symbiotic binary, where the WD accretes from the AGB star. The X-ray and UV luminosities, as well as the X-ray spectra of Y Gem resemble those of β/δ-type WD-symbiotic stars. The low X-ray-to-UV luminosity ratio suggests the accretion is likely in the form of an accretion disk.
3. We perform numerical simulations to investigate possible evolution paths of Y Gem. The results suggest that the mass transfer in Y Gem is most likely in the form of WRLOF accretion. The AGB star is now an M-type AGB star with a mass of ∼0.61−1.77M e , a bolometric luminosity of ∼4640−6960L e , and an effective temperature of Figure 5. A possible evolution path of Y Gem for the WRLOF scenario. The system initially contains a 0.8M e WD and a 1.2M e zero-age MS companion in a 10 3.5 (∼ 3162) day circular orbit. The top left panel shows the mass of the companion as a function of the stellar age. The remaining five panels show the effective temperature (top right), luminosity (middle left), surface gravity (middle right), wind-accreted rate (bottom left), and orbital period (bottom right) as a function of the companion mass. In all the panels, the systems begins its evolution from left and evolves to the right. Red curves represent the Y Gem-like stage and the purple dot represents the best match.
Our finding implies that there may be more such systems that have been classified as AGB-MS binaries. We propose further observations to put tighter constraints on the parameters of Y Gem, e.g., the bolometric luminosity and the orbital period to better understand the nature of Y Gem. Statistical investigations combining X-ray, UV, and parallax measurements (e.g., with Gaia) are also necessary to explore the WD-symbiotic star population and their relative contribution to SNe Ia.