Stellar properties of an actively accreting Algol-type eclipsing binary 2M16212643+2136590

Interacting binary stars undergo evolution that is significantly different from single stars, thus, a larger sample of such systems with precisely determined stellar parameters is needed to understand the complexities of this process. We present an analysis of a hierarchical triple containing a spectroscopically double-lined eclipsing binary, 2M16+21. Our calculations show that this system has undergone significant mass transfer, with the current mass and radius of the donor of 0.33 Msun and 2.55 Rsun, as well as the accretor of 1.37 Msun and 2.20 Rsun, resulting in a mass ratio of 4.2. Despite the already significant mass loss from the donor, shedding well over half its initial gas, mass transfer remains active. The shock from the accretion has produced a spot on the surface of the accretor that is ~2 times hotter than the photosphere, reaching temperatures of ~10,000 K and producing significant UV excess. This shock temperature is comparable to what is seen in the pre-main sequence stars that undergo active accretion. The compactness of the hot spot of just ~2 deg is one of the smallest observed in systems exhibiting binary mass transfer, pointing to the recency of its formation, as such it can be used to explicitly trace the point of impact of the accretion stream. The donor of this system may be a sub-sub-giant; comparing it with systems with similar initial conditions may help with understanding the formation processes of such stars.


INTRODUCTION
Eclipsing binaries (EBs) are valuable laboratories for understanding stellar evolution, as they offer a direct method probe of masses and radii of the individual stars that is independent of any modeldependent assumption.However, evolution of such systems can also be fundamentally different than in single stars.In particular, close binary systems can undergo transfer of mass from one star onto another if one of the stars fills out its Roche lobe, the size of which depends on the proximity of the two stars and their mass ratio.The resulting distortion of these stars can be observed in the light curves of EBs.If only one star fills out its Roche lobe, such systems are typically referred to as semi-detatched or Algol-type eclipsing binaries.
Mass transfer in binary systems can occur at different states of stellar evolution.Case A includes systems consisting of main sequence stars, Case B involves systems that have begun following exhaustion of H in the core and transition of the star into the red giant phase, and Case C following exhaustion of He (Ivanova 2015).Similarly, it is possible to catch systems at different times from the onset of mass transfer.Depending on this, mass accretion rates could vary wildly, ★ E-mail: marina.kounkel@unf.eduranging from 10 −4 to 10 −12 M ⊙ year −1 (van Rensbergen et al. 2011;Manzoori et al. 2017).
In order to understand these complex processes, a sample of multiple systems consisting of stars of different masses and evolutionary stages with precisely determined stellar parameters is needed.
Previously, we have modelled 2M17091769+3127589, a doublelined eclipsing binary system (Miller et al. 2021).This system was notable for having the more luminous star to have mass almost 6 times smaller than its fainter companion, which can be a byproduct of mass transfer from an evolved giant onto a main sequence star.The initial characterization of this system was serendipitous: it was one of the few eclipsing binaries identified in the ASAS-SN survey at the time (Pawlak et al. 2019) for which there were a sufficient number of high resolution spectra taken at different dates in which the star appeared as a double lined spectroscopic binary (SB2) that enabled to measure radial velocities (RVs) of both stars.
In this paper, we perform modeling of a system that is a twin to 2M1709, originating with stars of similar masses and comparable separation, but caught in the earlier stage of its evolution.This system, 2M16212643+2136590 (hereafter 2M16+21), is a semi-detatched eclipsing binary, and it still appears to have a significant ongoing accretion.In Section 2 we describe the available photometric, spectroscopic, and astrometric data that are currently available towards 2M16+21.In Section 3 we model the physical and orbital properties of 2M16+21, including the masses and radii of the individual stars.In Section 4 we analyze the accretion signatures, and construct a model of its evolutionary history.Finally, in Section 5 we conclude our findings.

Radial velocity
APOGEE is a multi-object spectrograph with high resolution (R∼22,500) covering the wavelength range of 1.51-1.7  (Majewski et al. 2017).Kounkel et al. (2021) have recently performed a systematic search of SB2s in the APOGEE data.They were identified using Gaussian fitting to the resulting cross-correlation function to measure both the RVs of the individual stars.When multiple spectral visits were available, mass ratio () was also estimated using the resulting RVs.This makes it possible to identify unusual systems, such as those with  > 2, which signifies brighter star in the system being significantly less massive.Although there are a number of such candidates, only 3 appear to be eclipsing binaries in TESS data, and only one, 2M16+21, has a sufficient number of APOGEE epochs to fully sample the RV phase space.
We adopt RVs for 2M16+21 from Kounkel et al. (2021).A total of 5 epochs are available, observed between June 2016 and February 2017.The source can be resolved as an SB2 during 3 of these epochs.Both sources are characterized by a wide profile attributable to fast rotation, with the full width of half maximum (FWHM) of the crosscorrelation function of 58 km s −1 for the more massive star, and 68 km s −1 for the less massive star.

Light curves
2M16+21 has a light curve morphology of an Algol-type EB with the period of ∼2.6 days.It has been observed by TESS in Sectors 25 (May 2020), and 51 (April 2022).However, Sector 51 has a significant portion of the data flagged as anomalous.As during a single sector, 2M16+21 is able to complete a full orbit 10 times, the light curve is well-sampled even with just one month of observations, as such we restricted our analysis solely to Sector 25.
The light curve was generated from TESS Full Frame Images using eleanor (Feinstein et al. 2019).The aperture was manually defined across 4 pixels that have shown the strongest periodic signal consistent with the source.We used CORR_FLUX, as it produced the cleanest light curve without any systematic trends and only minimal aberrations that were manually removed.
However, 2M16+21 is located in a somewhat crowded field.A nearby source, TIC 458723453 is located 8" away, which is not resolved in TESS due to its large pixel scale of 21".Modeling the pixel response function of all of the field stars in the region using Gaia photometry, we estimate that as much as 10% of all flux in the aperture is contamination, which is important to consider in the modeling of the lightcurve.
We extract ASAS-SN light curve for 2M16+21.It has been observed from February 2013 through August 2018 in V band, and from September 2017 onward in g band.Both light curves show similar morphology to the one in TESS, without any long term trends.However, some isolated epochs (most prominently in V band) show excessive level of noise, they were manually excluded from the analysis.Contamination from TIC 458723453 may also be a concern for ASAS-SN, as its pixel scale is 8" (Kochanek et al. 2017).
Finally, epoch photometry for 2M16+21 is also available as a part of Gaia DR3 (Eyer et al. 2022), in , , and  bands.Although only 71 epochs are available, they provide a comprehensive phase coverage of the light curve folded by the period.We also include these data in the analysis, as they are able to provide unblended fluxes unaffected by the nearby field stars.

Stellar properties
Parallax of 0.7485±0.1031mas for 2M16+21 was reported in Gaia DR3, however it also has an extremely high renormalized unit weight error (RUWE) of 6.482.Indeed, in Gaia DR2, parallax was reported to be 0.0906±0.1198mas, showing significant evolution with more data, as such either of these distance measurements are unlikely to be reliable.RUWE> 1.4 are considered as spurious astrometric solutions, typically thought to be affected by the orbital motion of multiple systems (Gaia Collaboration et al. 2021).However, the eclipsing binary itself is unlikely to be responsible for high RUWE: given the compactness of the orbit their orbital motion is unlikely to cause significant astrometric jitter.Rather, we suspect that this system is likely to be a hierarchical triple, housing a fainter and more widely separated companion.
As the system is a spectroscopic binary, spectroscoically derived properties are likely to be imprecise, but, nonetheless, they offer a sense of likely parameter ranges.Different pipelines analyzing APOGEE data report  eff ranging between 4500 and 5400 K (e.g.Ting et al. 2019;Olney et al. 2020;Sprague et al. 2022;Abdurro'uf et al. 2022), with the large scatter driven by different pipelines being sensitive to various features associated with both the primary and the secondary.Optical spectra for 2M16+21 have also been obtained by LAMOST, providing somewhat more stable  eff of 6240-6280 K (Luo et al. 2019;Xiang et al. 2019;Luo et al. 2019), most likely primarily sensitive towards the hotter star.Typically reported log  range between 3.5-4 dex.There is a significant scatter in the reported metallicitities, -0.1-0 dex from APOGEE, and 0.3-0.5 dex from LAMOST spectra.
We can achieve stricter constraints on  eff of the stars through performing an spectral energy distribution (SED) fit to the system considering two stars 1 .We follow the procedure outlined in Stassun & Torres (2016).Using synthetic stellar atmospheres from Kurucz   1993), we pick two templates of a given temperature, with log =3.5 and solar metallicity.The combination of these atmospheres are evaluated against the SED of the system, consisting of low resolution XP spectrum from Gaia, which spans the range in wavelength of 336-1020 nm, as well as the fluxes in Johnson and Cousins filters, alongside with fluxes 2MASS, WISE, SDSS, Gaia, and GALEX surveys.However, we include only NUV flux from GALEX, as FUV shows excess relative to the model in all of the fits.
We solve for radius of individual stars, distance and the extinction along the line of sight, limited to the maximum of   =0.153 mag, which is the integrated   along this line of sight from the map of Schlegel et al. (1998).Each combination of temperatures of the atmospheres is evaluated independently, and the lowest  2 of the fit is recorded.The best fit was produced by  eff,1 =6200 K,  eff,2 =4300 K, and angular radii of 14.0 and 11.4 as respectively.We note that this  eff,1 well matches the spectroscopically derived  eff from the optical spectra, and  eff,2 is somewhat consisted with the near-IR spectra.A family of similar solutions is also possible (Figure 2).The fit suggests that the hotter star has to have a larger radius than the cooler star for all  eff,2 <6300 K.

MODEL
Given the availability of the light curves that show the system as eclipsing, as well as RVs for both stars, it is possible to perform a full fit of all of the stellar parameters.To do this, we use PHOEBE (Conroy et al. 2020), which enables a comprehensive modeling of eclipsing binaries.In this section we describe different approaching to constructing such a model.

Simple model
In the process of modeling 2M16+21, we consider 21 parameters in total.They include 9 stellar parameters: effective surface temperatures  eff,1 ,  eff,2 , semi-mejor axis , orbital period , inclination , mass ratio , center-of-mass velocity , time of passage of the superior conjunction  0 , as well as the equivalent radius  1 of the hotter star (accretor).The system is treated as semi-detached, with  2 being set to the radius of the Roche lobe of the cooler star (donor).
Due to the compactness of the orbit, circular orbit is assumed.
In addition to these stellar parameters we fit the third light contamination ( 3 ) for each of six light curves in different bandpasses that are available for the star, to accommodate flux contamination from the nearby field stars, as well as from the likely tertiary within the system.We also include passband luminosity for all six light curves, which is a nuisance variable recommended by PHOEBE.
We performed an initial manual exploration of the ranges of stellar parameters, as well as a preliminary optimization with the Nelder-Mead solver.Due to the combination of stellar parameters required by this system being outside of the set of the set of the atmospheres incorporated into PHOEBE, we treat atmospheres of both of the stars as blackbodies, with limb darkening computed using LDTK (Parviainen & Aigrain 2015) from the larger set of PHOENIX atmospheres (Husser et al. 2013).
We then fit the system using the emcee sampler.It was initialized with 250 walkers, and a chain length of 5,000 steps, which was sufficient for emcee to converge after 2,000 steps.In the distribution, we adopted  and  and their uncertainties from Kounkel et al. (2021).
The resulting fitted parameters are also listed in Table 1, and the evaluation of the model against the data is shown in Figure 3.The fit well describes both the RV curve and the light curves, although it doesn't capture a slight asymmetry in the continuum flux in TESS.
The third light contribution in TESS is comparable but somewhat larger than the a-priori estimate of 10% from the neighboring field stars alone.Similarly, although Gaia has a significantly higher spatial resolution, it also requires a non-negligible  3 of ∼ 2%.This supports a presence of a significantly fainter third star in the system.
However, adopting these stellar parameter for the SED fit creates a significant discrepancy: the model significantly overestimates IR flux, and underestimates blue flux.This is partially due to the fact that temperatures of the two stars converge to a lower  eff than the ranges that are able to fit the SED.Forcing a tighter prior on  eff (such as at  eff,1 >6000 K) results in PHOEBE converging its  eff right at the edge of this prior.The light curve modeling often has greater sensitivity to the ratio of the temperatures than the precise temperatures themselves, but even scaling  eff to the range in Figure 2 cannot produce a good SED fit.The reason for this is because PHOEBE requires  2 >  1 in all permutations, inconsistent with the SED derived ratio of the radii at the point of intersection of  eff ratios between the two approaches.
As such, the model requires an additional hotter source of light that would be able to reconcile the two approaches.This system is likely a tertiary, but the third star cannot be responsible for producing excess blue light: fitting the SED with three stars, holding two of them fixed at PHOEBE solutions results the third star being at least a factor of 3 times brighter than either of the stars in the eclipsing binary, which is inconsistent with  3 estimates, and is unlikely to not have been reflected in the spectra of the system.3, but for a model that considers a presence of a hot spot; green curve in the SED panel is the spot contribution.Note that the fit for T band is improved over Figure 3, at a cost of significant overestimation of UV flux.
There is one additional discrepancy between the derived model and the data.Regardless of the adopted  eff,1 , the donor star is expected to have comparable flux to the accretor in NIR, but it is expected to be fainter at all wavelength.In contrast, APOGEE spectra suggest that the donor star to be more luminous in the H band.Even with the uncertainty driven by the third light, individual radii, as well as the ratio of  eff,2 / eff,1 are strongly constrained by the available light curves, which in turn constrain the luminosity of stars at a given wavelength.The mostly likely explanation for this is that in performing cross-correlation, lower  eff template of the donor star was favored as at lower  eff various spectral lines are stronger.Because of  eff difference of >1000 K, the template mismatch resulted in a weaker cross-correlation strength for the hotter accretor star, despite the latter being marginally more bright in this wavelength range.

Spotted model
In order to reconcile parameters from SED and light curve fitting, we consider the effect of spots.To increase the fraction of the blue light in the system, the spot has to be hotter than the photosphere.As the system is semi-detatched, likely with active mass transfer, the shock from the accretion stream from the donor is able to excite a small fraction of the photosphere of the accretor (Kouzuma 2019).Because of this, the accretor is more likely to contain a necessary hot spot.
The hint of such a hot spot is present in the data.As previously mentioned, TESS light curve exhibits a slight asymmetry in its flux continuum that cannot be reproduced by a simple model.Furthermore, there is FUV excess in the SED fit that cannot be accounted by any combination of  eff and radii of the two stars: both of these factors can be explained by a presence of a spot.
Given that there is no additional periodicity other than the orbital period, it is likely that the rotation of the stars is tidally locked to their orbit.Both compactness of the orbit and the mass transfer can be responsible for this.With the orbital period of just 2.7 days and the radii of ∼ 2.5 R ⊙ , this results in rapid rotation, with velocities on the order of 40-50 km s −1 , consistent with the velocity spread seen in the CCF taking instrumental broadening into the account (Figure 1).
To incorporate a spot in the SED fitting, we create a mesh in PHOEBE, which is then able to independently specify the temperature, surface gravity, and area within each of the facets of the mesh.SED is then constructed for all of the facets, and they are integrated together to produce the total flux for each star and the system as a whole.
We first explore the range of likely spot properties in the SED fitting.We randomly sample  eff,1 between 5500 and 6500 K, forcing  eff,1 =1.53 eff,2 , to match the previously derived temperature ratio, preserving PHOEBE derived radii from the spotless model.Spot sizes are drawn from 0 to 25 degrees, and temperature contrast from 1 to 2.5.Once again,  eff,1 <6000 K is strongly ruled out by this fit.
Larger spots require spot temperature closer to  eff , smaller spots have to be significantly hotter (Figure 4).
Similarly, keeping  eff fixed, we explore the dependence of the spot's size and contrast relative its position on the photosphere on the goodness of fit of both the light curves and the SED simultaneously by computing  2 for each.The most likely solution is a compact equatorial spot ∼140 • away from the point facing the donor (Figure 5. To better fit the UV excess, small spots are preferred, with temperature > 1.5 hotter than the photosphere; this is consistent both for the SED only and the light curve only fit (Figure 6).
We then initialize emcee with 250 walkers, adding temperature contrast, radius, and the position of the spot as parameters in addition to 21 parameters from the spotless model.The resulting parameters are listed in Table 1, and the fit is shown in Figure 7.
In large, the quality of the fit relative to the spotless model is very similar -adding a spot has significantly improved the fit to the TESS light curve, however, it did not substantively alter any of the derived parameters compared to a spotless model, thus there isn't  a substantive differences in the overall shape of the modelled light curves or the RV curves.As such, the quality of the resulting SED fit did not substantively improve.The model still overestimated IR flux because the derived  eff,1 was too small.And while the NUV flux is better fit, FUV flux is significantly overestimated, because the contrast ratio and/or the spot size appear to be too large to well fit the available photometry.
We subsequently attempted several permutations of the fitting procedure, placing different priors on  eff,1 and on the spot properties.Forcing hotter  eff,1 fixes discrepancy in IR, but not in UV, as the spot continued to be too luminous.Eventually we fixed these parameters to the values most likely required by the SED, namely  eff,1 =6150 K, spot size of 2.22•, and the temperature contrast of the spot of 1.8 (i.e.,  spot ∼11,000 K).The parameters are listed in Table 1, and the fit is shown in Figure 8.The quality of the resulting fit to the light curve is slightly impacted in comparison to the version of the model in which  eff is unconstrained, but this allows a more self-consistent solution between both approaches.The distance to the system based on the SED fitting is 782±13 pc.
The favored spot location is different between the models with unconstrained and with fixed  eff .Regardless of  eff,1 model with unconstrained spot parameters favors a larger equatorial hot spot.Forcing the spot to be less luminous forces it to migrate north.The smaller it is, the more likely model would place it near the polar region.The final model with fully constrained spot temperature and size places it only 6.5 • from the north pole. .SED fit to the system from Figure 8, with the green curve shows the included teritary component consisting of a main sequence star with  eff =4600 K.

Third Star
As previously mentioned, very large RUWE, as well as nonnegligible third light contribution in Gaia passbands is indicative of the presence of a tertiary star in the system.Since it amounts to only 2.6% of the total flux at RP band, it does not significantly alter the SED fit, but it is likely to be a low mass main sequence star with  eff ∼4600 K with the corresponding radius of 0.66 R ⊙ , which matches the RP flux contribution (Figure 10).However, such a star significantly overestimates third light at G and BP passbands.It is possible that the model underestimates this contribution.Indeed, in all of the models, the fitted light curve somewhat underestimates the depth of the secondary eclipse in these bands, despite matching well in all of the other passbands.
This may be what has fundamentally led to the difficulty of deriving accurate  eff within PHOEBE, and requiring constraining it primarily from the SED fit.For example, constructing a PHOEBE model fit using only T, V, and g light curves, which are all contaminated by TIC 458723453 (and disregarding Gaia lighcurves) produces third lights of  3,T =8%,  3,V =0.1%, and  3,g ∼0 (much smaller than  3,T =12% and  3,V,g =6% in the model with the Gaia data).This model also had incorrect  eff , favoring  eff,1 >7000 K, which could not be reconciled with the SED fit.
It appears that PHOEBE struggles to accurately represent  3 in cases where all of the passbands have some degree of contamination from a third body, favoring the lightcurve with the smallest contamination to be set at  3 ∼0, regardless of the absolute fraction, and this does impact the fitted  eff .
As such,  3 in all passbands may still be underestimated by a few %.In this case, the tertiary companion would likely be a main sequence 4600-5000 K star.

Accretion status
The current masses of the donor and the accretor star are 0.33 and 1.368 M ⊙ .In the event of the conservative mass transfer and assuming that the masses two stars may have been initially similar, the minimum initial mas of the donor would be >0.85M ⊙ .As such, the donor has already lost well over half of its initial mass.We note that non-conservative mass transfer is common, including in the lower mass stars, with mass frequently lost due to winds, jet formation, outflows from the disk, etc (e.g., Woods et al. 2012;Chen et al. 2020;Sun et al. 2021;Sun & Mathieu 2023), thus, this is a lower limit.Given the proximity of the two stars, the onset of it would have been rapid following the expansion of the donor into a red giant.What is surprising, however, is that despite the history of substantial mass transfer in this system, it appears to still be active, as indicative through UV excess and the presence of a hot spot.
In young stellar objects that still actively accrete gas from protoplanetary disks, gas flows along the magnetic fields of a (proto)star, shocking the photosphere upon impact (e.g., Hartmann et al. 2016).Such shock has been observed to create hot spot with the temperature contrast of >2 relative to photosphere ( spot ∼ 9000 K in a young star with Espaillat et al. 2021). spot can reach upwards of ∼ 12,000 K (de Sá et al. 2019).This is not dissimilar to  spot ∼ 11,000 K that we observe in 2M16+21.
In young stars, high density of the accretion flow tends to also produce spots with low filling factor of 0.1-1% of the total size of the photosphere (Calvet & Gullbring 1998), which can correspond to the spot sizes of a few degrees.On the other hand, lower density accretion columns will produce spots with larger filing factor (Ingleby et al. 2013).
Hot spots have been previously observed in a number of close binaries (Kouzuma 2019).However, in systems with a similar configurations as 2M16+21, i.e., a semi-detached binary in which a less massive star is overflowing its Roche lobe (SD2), contrast ratio of the spot to the photosphere is typically 1.1, and the spots themselves tend to be large, 10-30 • .As such, high  spot relative to the temperature of the photosphere, as well as the compactness of the hot spot make 2M16+21 relatively unusual in its class.Most likely, this difference is caused by significantly higher  of 2M16+21 in contrast to a number of other semi-detached eclipsing binaries.This may be due to this system being caught very early on its evolution: after the bulk of mass has already been transferred between the two stars, but not yet at a point where it reaches a relative equilibrium in their respective masses.
2M16+21 appears to be very similar to 2M17091769+3127589 (Miller et al. 2021) in terms of the total mass of the system as well as the mass ratio.However, while donor of 2M17+31 is an evolved red giant with  ∼4 R ⊙ , 2M16+21 appears to have only recently begun climbing up the red giant branch.

Origin of discrepancy in the resulting PHOEBE & SED models
As previously discussed in Section 3.3, the difficulty in determining accurate  eff in PHOEBE has most likely been caused by third light contamination, namely by underestimation of third light in G and BP bandpasses.However, even when holding  eff fixed in the model, there is a difference in the derived parameters for the spot between SED and PHOEBE fit.PHOEBE favors significantly more luminous spot, and, constraining the spot size and contrast results in a very different spot position.
In semi-detatched eclipsing binary systems, hot spots are most commonly equatorial (Kouzuma 2019), which is the location favored by the light curve in 2M16+21.This gives greater credence to the PHOEBE derived spot properties, at least concerning its position.The luminosity of the spot may be skewed for the similar reason why  eff of the photosphere appears to be incorrect.
There is also a possibility of both SED-only and PHOEBE-only fit of the spot luminosity being accurate.The primary effect on photometry of a more luminous spot is in the FUV bandpass from GALEX.The only light curve with significant sensitivity to the spot is from TESS, as ASAS-SN has significantly lower signal-to-noise, and Gaia light curve is very sparse.GALEX and TESS observations have been separated in time by approximately a decade.During this time, there could have been an evolution in its accretion/spot properties, which may have increased the FUV luminosity, even though it is not captured in the available data.

Possible evolutionary history of this system
To explore the past and future of this system, we utilized the Modules for Experiments in Stellar Astrophysics (MESA; Version 12115; Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 2019;;Jermyn et al. 2023) binary module.The simulation adopts an initial metallicity of  = 0.019.The model incorporates the same physics as described in Section 4 of Miller et al. (2021).Our initial goal is to fit the MESA model to the fixed  eff result from Table 1, on the mass versus orbital period plane.
Considering 2M16+21 has very similar features to 2M17091769+3127589, the initial parameters of the model were searched near the best-fit model of 2M1709.Figure 11 shows the best fit model for 2M16+21 with an initial donor mass of   = 1.2  ⊙ , an accretor mass of   = 1.15  ⊙ , and an initial orbital period of  orb = 3.32 days (comparable with initial donor mass of   = 1.2  ⊙ ,   = 1.11  ⊙ , and an initial orbital period of  orb = 3.43 days in 2M1709).To match the observed mass of the accretor star, the mass transfer efficiency is fixed at 23% during the mass transfer (lower than 49% in 2M1709).The system undergoes initial shrinkage during its evolution.The onset of mass transfer occurs at  orb = 1 day, when the donor star has already exhausted its central hydrogen.Around  orb = 0.9 day, with  d = 0.97  ⊙ and  a = 1.2  ⊙ , the system begins to expand in  orb and eventually reaches its current binary configuration, with a mass loss rate from the donor star of  ∼ 10 −9.5  ⊙ /yr.At the same time, the accretor star is in the subgiant phase, as indicated by both the modeled and observed radius and  eff .At the end of the simulation, the donor star evolves into a 0.23  ⊙ helium pre-white dwarf, while the accretor star becomes a 1.36  ⊙ red giant star.The simulation terminated when the accretor fills the Roche lobe radius, resulting in a contact binary.Most likely, after that point the system will evolve into a planetary nebula, as there is not enough mass to produce a Type 1a supernova.
This model explains the current masses, orbital period of the system, and the  eff and radius of the accretor star.Figure 12a illustrates the evolutionary track of the accretor star in terms of radius and  eff .The model begins at the zero-age main-sequence stage of the stars, characterized by ( eff , ) = (6070K, 1.1 ⊙ ).Mass accretion commences during the accretor's post-main-sequence evolution at ( eff , ) = (6000K, 1.6 ⊙ ).At the observed  orb = 2.66 days, the evolutionary track passes through the box indicating the error bar of the measurements, corresponding to a mass of 1.35  ⊙ .

𝑇 eff of the donor star
The temperature of the donor star is significantly hotter than observed, with a difference of approximately 800 K, as depicted in Figure 12b.The evolutionary track of the donor star commences at ( eff , ) = (6200K, 1.12 ⊙ ).The onset of mass transfer takes place at ( eff , ) = (5550K, 2.12 ⊙ ).The model indicates a case B mass transfer scenario.The discrepancy of 800 K between the model and data may not be solely attributed to model fine-tuning issues.We explored alternative parameter spaces, such as using a lower mass donor, but the  eff during mass transfer remains hotter than observed.None of the MESA models we have explored were able to more optimally fit the parameters of the system.It is unclear why there is such a significant discrepancy in the donor's  eff .One of the possibilities might have been due to the fast rotation of the donor produce gravitational darkening that caused the star to appear cooler in comparison to a non-rotating model.This, however, is unlikely, as the difference of ∼1000 K would require rotation rate of 0.9 times the critical (Paxton et al. 2019, see Figure 38); the rotation velocity of ∼60 km s −1 observed for this system is significantly smaller.
The donor star may be a sub-sub-giant (SSG), which are often identified though being luminous but cooler than the red giant branch  1.The size of the box correspond to the size of the error bars, while the color represents the mass of the accretor throughout its evolution.Bottom: same as above but for the donnor star.
stars (Mathieu et al. 2003).Much about their formation is still uncertain, and evolutionary models still struggle to fully reproduce their evolutionary history.SSGs may often have large cool spots on their surface which may partially explain  eff discrepancy (Stassun et al. 2023) -the light curves for 2M16+21 show no evidence of these cool spots, as a single hot spot is sufficient for accounting for the out-of-eclipse variability.Nonetheless, the data may not necessarily rule out a possibility of them being present.Given the similarity of 2M16+21 to 2M1709, it is possible that the donor star being significantly cooler than what is predicted by the model is a transitory stage, and as it expands further, it may return to its expected position on the red giant branch.Further resolving the temperature discrepancy of the donor star would be the focus of the future work.We note, however, given the strong presence of the accretion signatures that were discussed in Section 4.1, combined with an imprecise fit, it is possible that  in the model is underestimated due to the limitations of the evolutionary model being only 1D.

CONCLUSIONS
Through analyzing light curves from TESS, Gaia, and ASAS-SN, as well as radial velocities from APOGEE, we derive stellar and orbital properties of a semi-detatched eclipsing binary system, 2M16+21.This system has a significant history of mass transfer, with the donor star transferring well over half of its mass onto the accretor, resulting in the present day masses of 0.33 M ⊙ and 1.37 M ⊙ respectively.Despite his, the system still exhibits active accretion, as is evident by a very compact and very luminous hot spot on the surface of the accretor.
There have been some challenges in deriving accurate  eff for the individual stars due to third light contamination, both from a low mass main sequence companion orbiting around the binary, as well as a nearby field star TIC 458723453.Deriving accurate  eff for both stars required an iterative process of both the light curve and the SED modelling, as the latter can place firmer constraints on the allowable range of  eff .Incorporating SED fitting into the light curve fitting routines such as PHOEBE could yield more self-consistent models.
We note that the donor star in 2M16+21 appears to be a sub-sub-giant.However, this system appears to be a twin to 2M17091769+3127589 (Miller et al. 2021), potentially starting out with near-identical initial conditions, but significantly less evolved.Further analysis of these two systems may shed more light on the formation and evolution of SSGs.

Figure 1 .
Figure 1.APOGEE cross-correlation function of 2M16+21 showing both stars.The more pronounced blueshifted peak corresponds to the lower mass star.FWHM of the two components are 68 and 58 km s −1 .

Figure 2 .
Figure 2. Left: Goodness of fit of the SED using a combination of two model atmospheres of different  eff .The white line shows the best family of fits.The purple line shows the  eff ratio of 1.53 derived from the light curve analysis.Right: radius ratio of the two stars produced by the best fit for a given  eff combination.

Figure 3 .
Figure 3. Fit of the spotless model.Left panel shows the 6 different light curves used in the analysis, with the red line being the model, and black points being the phase-folded data.Top right panel shows the radial velocity curves for the primary and the secondary.Bottom right panel shows the SED fit, with black points corresponding to the broad-band photometry, black curve showing the Gaia XP spectrum, yellow curves showing the best-fitted SEDs of individual stars, red is the coadded SED, and blue dots show the interpolated fluxes of the model at the bandpasses corresponding to the data.Residuals between the model and the data are shown below; for light curves residuals are scaled proportionally to their respective uncertainties.

Figure 4 .Figure 5 .
Figure 4. Goodness of fit of the SED of the spotted model based on  eff,1 .

Figure 6 .
Figure 6.Goodness of fit of the spottded model with variable hot spot size and temperature contrast.Top: SED-only fit.Middle: Light curve only fit.Bottom: Joint solution.

Figure 7 .
Figure7.Same as Figure3, but for a model that considers a presence of a hot spot; green curve in the SED panel is the spot contribution.Note that the fit for T band is improved over Figure3, at a cost of significant overestimation of UV flux.

Figure 8 .
Figure 8. Same as Figure 7, but with fixed temperatures of both the spot and the photospheres.

Figure 9 .
Figure 9. Mesh showing the temperature distribution along the photospheres, including a spot.Top panel shows the unconstrained spotted model.Bottom panel shows the model with fixed  eff of the photosphere and of the spot.Note that in the latter case the spot is close to polar, highlighted by a black rectangle.
Figure10.SED fit to the system from Figure8, with the green curve shows the included teritary component consisting of a main sequence star with  eff =4600 K.

Figure 11 .
Figure11.The MESA model with initial donor mass  d = 1.2  ⊙ , accretor mass  a = 1.15  ⊙ and initial orbital period  orb = 3.32 days, with a mass transfer efficiency of 23%.In the upper subplot, the blue line represents the mass loss rate from the donor star, while the orange line represents the mass accretion rate of the accretor star.The lower panel illustrates the evolution of the donor (blue line) and accretor masses as a function of  orb during the binary evolution.The error bar of the donor star mass measurement is smaller than the size of the data point, so it is not displayed on the plot.

Figure 12 .
Figure12.Top: The evolutionary track of the fitted MESA model for the accretor star on the radius versus  eff plane.The data represented as red dot is from the third column of Table1.The size of the box correspond to the size of the error bars, while the color represents the mass of the accretor throughout its evolution.Bottom: same as above but for the donnor star.