# Detached Eclipsing Binaries in Compact Hierarchical Triples: Triple-lined
systems BD+44 2258 and KIC 06525196

Ayush Moharana,1 K. G. Hełminiak,1 F. Marcadon,1,2 T. Pawar,1 M. Konacki,3 N.
Ukita,4,5 E. Kambe,6 and H. Maehara4
1Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul.
Rabiańska 8, 87-100 Toruń, Poland
2Villanova University, Dept. of Astrophysics and Planetary Sciences, 800 East
Lancaster Avenue, Villanova, PA 19085, USA
3Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul.
Bartycka 18, 00-716 Warszawa, Poland
4Okayama Astrophysical Observatory, National Astronomical Observatory of
Japan, 3037-5 Honjo, Kamogata, Asakuchi, Okayama 719-0232, Japan
5The Graduate University for Advanced Studies, 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan
6Subaru Telescope, National Astronomical Observatory of Japan, 650 North
Aohoku Place, Hilo, HI 96720, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Compact Hierarchical Triples (CHT) are systems with the tertiary star orbiting
the inner binary in an orbit shorter than 1000 days. CHT with an eclipsing
binary as its inner binary can help us extract a multitude of information
about all three stars in the system. In this study, we use independent
observational techniques to estimate the orbital, stellar, and atmospheric
parameters of two triple-lined CHT: BD+44 2258 and KIC 06525196. We find that
the masses of stars in BD+44 2258 are $1.011\pm 0.029M_{\odot}$, $0.941\pm
0.033M_{\odot}$, and $0.907\pm 0.065M_{\odot}$ while in KIC 06525196 the
estimated masses are $1.0351\pm 0.0055M_{\odot}$, $0.9712\pm 0.0039M_{\odot}$,
and $0.777\pm 0.012M_{\odot}$. Using spectral disentangling, we obtained
individual spectra of all the stars and combined it with light curve modelling
to obtain radii, metallicities and temperatures. Using stellar evolution
models from MESA, we constrain the log(age) of BD+44 2258 to be 9.89 and 9.49
for KIC 06525196. Two stars in BD+44 2258 are found to be sub-giants while all
three stars in KIC 06525196 are main-sequence stars. We constrain the mutual
inclinations to certain angles for BD+44 2258 and KIC 06525196 using numerical
integration. Integrating with tidal interaction schemes and stellar evolution
models, we find that KIC 06525196 is a stable system. But the inner binary of
BD+44 2258 merges within 550 Myrs. The time of this merger is affected by the
orientation of the tertiary, even rushing the collapse by $\sim 100$ Myrs when
the mutual inclination is close to 90 degrees.

###### keywords:

binaries: eclipsing – binaries: spectroscopic – stars: fundamental parameters
– stars: evolution – stars: individual: BD+44 2258, KIC 06525196 – stars:
kinematics and dynamics

††pubyear: 2022††pagerange: Detached Eclipsing Binaries in Compact
Hierarchical Triples: Triple-lined systems BD+44 2258 and KIC 06525196–D

##  1 Introduction

The multiplicity of stars is a well-established phenomenon (Duchêne & Kraus,
2013). The incidence of multiplicity varies with the spectral type of the
stars. Multiplicity is 44% among the Solar-like stars, out of which $8\%$ are
triple stars (Raghavan et al., 2010). These numbers increase in O, B, and A
types (Mason et al., 2009; Shatsky & Tokovinin, 2002; Kobulnicky & Fryer,
2007). There have been a lot of studies towards understanding binaries which
has created an almost complete picture of their evolution and formation. The
next step in decoding the multiple-architecture is understanding triple
systems.

Binary formation channels can be largely classified into disk instability,
core-fragmentation, and N-body interactions. The complexity of formation
increases in a triple system with a combination of these formation channels
being responsible for their formation (Tokovinin, 2021). Most of the studies
to understand these scenarios using orbital architectures, metallicity
variation among stars, and mass distributions in triples, have been usually
restricted to wide triples (Tokovinin, 2017, 2022; Lee et al., 2019).

The evolution of triple stars also departs from the simplified evolution of a
single star. A triple-star system has the additional complexity of multiple
dynamical interactions. The outer companion can directly alter the formation
of the components of the inner binary, leading to their apparent difference in
properties (e.g.apparent age; Marcadon et al. 2020), and has been used in
explaining various evolutionary phenomena of single and binary stars.

Hierarchical triple systems have been intensely studied in order to understand
the formation of closest main-sequence binary systems (Eggleton & Kiseleva-
Eggleton, 2001; Naoz & Fabrycky, 2014; Moe & Kratter, 2018). One of the
probable theories of the formation of blue straggler stars involves
perturbations from a third body (Perets & Fabrycky, 2009). One in a thousand
high-mass x-ray binaries evolve through interaction with a third star to form
low-mass x-ray binaries (Eggleton & Verbunt, 1986). Asymmetry of Planetary
Nebulae has been linked to evolution in a triple system (Akashi & Soker, 2017;
Jones et al., 2019) and has even been suspected to play a role in driving
white dwarf mergers towards type Ia supernova explosions (Maoz et al., 2014).

Recent population synthesis studies have shown that 65-75 % of triples undergo
mass transfer (Toonen et al., 2020). Hamers & Dosopoulou (2019) have shown
that this can occur uniformly throughout the orbit or at certain points due to
eccentricity and inclination changes known as von Zeipel-Lidov-Kozai (ZLK)
oscillations (von Zeipel, 1910; Lidov, 1962; Kozai, 1962). Furthermore, these
systems perturb Roche-lobe potentials and also undergo Roche-Lobe Over Flow
(RLOF), which can occur for the three individual stars and can be circumbinary
too. Therefore, understanding the stellar evolution coupled with the dynamical
evolution is important when studying triples. Most of the known triple systems
have long tertiary periods and therefore their dynamical effects can have
timescales of decades or centuries. There is a subset of these triples, called
Compact Hierarchical Triples (CHT), which offer more potential for
observational astrophysics (Borkovits, 2022). The timescale of the changes due
to these interactions is short in CHT and can be observed easily over a few
years, e.g., vanishing eclipses of HS Hya (Zasche & Paschke, 2012) or re-
appearing eclipses of V907 Sco (Zasche et al., 2023). These are triples with
the outer orbit period shorter than 1000 days. Due to this, dynamical
processes in CHTs can be observed in timescales of years. CHTs were thought to
be rare (Tokovinin, 2004) but with new space-based photometic missions we are
discovering more of these systems (Rappaport et al., 2013; Borkovits et al.,
2016).

Detached Eclipsing Binaries (DEB) are known as the source of the most accurate
stellar parameters (e.g., mass, radius, etc.). Accuracy of less than 1% can be
attained by coupling high-precision photometry and high-resolution
spectroscopy (Torres et al., 2010). The accuracy is robust and independent of
different models and methods, even varying slightly due to different numerical
implementations (Maxted et al., 2020; Korth et al., 2021).

If a CHT has a DEB as its inner binary, there is an added advantage of
obtaining accurate stellar parameters of not only the binary but of the
tertiary as well (Hełminiak et al., 2017). Using light curve modelling,
eclipse timing Variations, spectral analysis, and RVs, we can obtain an
accurate picture of the orbits, geometry, stellar parameters, metallicity,
age, and evolutionary status.

There has been a surge in interest in CHT recently. Kepler (Howell et al.,
2014) and Transiting Exoplanet Survey Satellite (TESS) both have been crucial
in detection and analysis of these systems (Borkovits et al., 2015; Borkovits
et al., 2020). Tertiary stars in CHT have been found to host tidally induced
pulsations (Fuller et al., 2013). Ongoing projects are using triply-eclipsing
triples (TET) to characterise CHT (Rappaport et al., 2022). Studies even show
that CHT can produce exotic Thorne-Żytkow objects (Eisner et al., 2022). Last
but not the least, CHT have proven to be useful to study ZLK oscillations and
their effect on stellar evolution (Borkovits et al., 2022). Though most of
these studies have provided us with mass and radii of all the stars in a CHT,
the tertiary radii-space is dominated by TET or planar systems.

In this paper, we report the detection of a DEB in a CHT, BD+44 2258
($\alpha=$ 13:15:06.66, $\delta=$44:02:33.48; hereafter BD44). BD44 has been
previously observed in UV (GALEX; Bianchi et al. 2011) and X-Ray (ROSAT; Voges
et al. 1999) but as a single source. We obtain stellar, orbital, and
atmospheric parameters of all the stars in BD44 and a previously detected CHT,
KIC 06525196 ($\alpha=$ 19:30:52.32, $\delta=$41:55:20.81; hereafter KIC65)
with TESS photometry and HIDES spectroscopy. We explain the different
observations and methods used for extracting parameters in Sections 2 and 3.
We use these parameters to estimate the age and evolutionary stages of the
components. Further, using the orbital parameters, we study the evolution and
stability of the systems as explained in Section 4.

##  2 Observations

###  2.1 Photometry

We use photometry from TESS (Ricker et al., 2015) for our light curves111From
GI programmes: G022003, G022062, G04047, G04171, G04234.(LC). BD44 (TIC
284595199) has 2-minute cadence photometry obtained from Sectors 16, 22, and
49. KIC65 was a target from the main Kepler mission field, and these data have
been analysed in Hełminiak et al. (2017). In addition to this, the TESS
2-minute cadence photometry (TIC 137757776) is available from Sectors 40 and
41. In this work, we only model a segment of TESS data from Sector 41 for
KIC65 while for BD44 we use segments from Sectors 16, 22, and 49. These
segments were selected on the basis of spot variability in the LCs. This was
also the reason for selecting only Sector 41 for KIC65 because the overall
structure of the LC was similar in both sectors but Sector 41 has fewer
fluctuations than Sector 40. The structure of the LC for BD44 changes in the
sectors and therefore we model all the sectors to check for consistency.

We filtered out the points which had the best quality-flag for our purposes.
There seems to be no other star in the Full Frame Image of BD44 but there
seems to be some contamination for KIC65. For BD44, we compared normalised LCs
obtained using different apertures on the FFIs only to find they all are
identical when normalised. Therefore we consider that any third light that
would show up from LC modelling will be solely due to the third star. We
considered the Simple Added Photometry or SAP fluxes for modelling both our
systems.

###  2.2 Spectroscopy

We use spectra collected with the HIDES (Izumiura, 1999; Kambe et al., 2013),
attached to the 1.88-m telescope at the Okayama Astrophysical Observatory.
Observations were conducted in the fiber mode with an image slicer ($R\sim
50\,000$), without I2, and with ThAr lamp frames taken every 1-2 hours. The
spectra are composed of 62 rows covering 4080–7538 Å, of which we use 30
(4365–6440 Å) for radial velocity (RV) calculations. A detailed description of
the observing procedure, data reduction and calibrations is presented in
Hełminiak et al. (2016).

For KIC65 we used exactly the same set of 14 spectra and RV measurements as in
Hełminiak et al. (2017). For BD44 we took a total of 28 spectra.

##  3 Analysis

In the following sections, we use A-B notation to denote the CHT, where B is
the tertiary. A is the eclipsing binary with components Aa and Ab. Aa
corresponds to the primary classified according to temperature, usually the
deepest eclipse in the LC if not affected by spots. Since we would be talking
about six stars in total (three in each of the two CHT), we use the short form
for the star’s name along with the alphabetical notation to exclusively denote
each star, e.g., the secondary (cooler star in the binary) of KIC65 is
referred as KIC65Ab.

###  3.1 RV extraction and fitting

Since both targets were observed with the same spectrograph as part of the
same programme, the approach to RV calculations and fitting was essentially
identical. It is described in detail in Hełminiak et al. (2017), but we
summarise it here briefly.

The RVs were calculated with a TODCOR method (Zucker & Mazeh, 1994) with
synthetic spectra computed with ATLAS9 code as templates. Measurement errors
were calculated with a bootstrap approach and used for weighting the
measurements during the orbital fit, as they are sensitive to the signal-to-
noise ratio (SNR) of the spectra and rotational broadening of the lines.
Though this code is optimised for double-lined spectroscopic binaries and
provides velocities for two stars ($u_{1},u_{2}$), it can still be used in
triple-lined systems as well. The RVs of the eclipsing pair were found from
the global maximum of the TODCOR map since in both targets these components
contribute more to the total flux than the third star. The tertiary’s
velocities were found from a local maximum, where $u_{1}$ was set for the
tertiary, and $u_{2}$ for the brighter component of the eclipsing pair. This
scheme was used previously for KIC65 in Hełminiak et al. (2017), and the $rms$
of the resulting tertiary’s RVs was comparable to the stability of the
instrument estimated from RV standard stars. In some cases, the velocity
difference between two components was too small to securely extract their
individual RVs, and measurements were taken only for the remaining one. For
this reason, each component of a given triple may have a different number of
RV data.

The orbital solutions were found using our own procedure called v2fit (Konacki
et al., 2010). It applies a Levenberg-Marquardt minimisation scheme to find
orbital parameters of a double-Keplerian orbit, which can optionally be
perturbed by a number of effects, like a circumbinary body. The fitted
parameters are: orbital period $P$, zero-phase $T_{P}$222Defined as the moment
of passing the pericentre for eccentric orbits or quadrature for circular.,
systemic velocity $\gamma$, velocity semi-amplitudes $K_{1,2}$, eccentricity
$e$ and periastron longitude $\omega$, although in the final runs, the last
two parameters were usually kept fixed to zero. We also included the
difference between systemic velocities of two components,
$\gamma_{2}-\gamma_{1}$, and a circumbinary body on an outer orbit,
parameterised analogously by orbital parameters $P_{3}$, $T_{3}$, $K_{3}$,
$e_{3}$, and $\omega_{3}$. In such case, $\gamma$ is defined in the code as
the systemic velocity of the whole triple.

Systematic errors that come from fixing a certain parameter in the fit are
assessed by a Monte-Carlo procedure, and other possible systematics (like
coming from poor sampling, small number of measurements, pulsations, etc.) by
bootstrap analysis. All the uncertainties of orbital parameters given in this
work already include the systematics.

In addition, for each observation where three sets of lines were sufficiently
separated, we also calculated the systemic velocities $\gamma(t_{i})$ of the
inner pair, using the formula:

$\gamma(t_{i})=\frac{v_{1}(t_{i})+qv_{2}(t_{i})}{1+q},$ (1)

where $v_{1,2}(t_{i})$ are the measured RVs of the inner binary, and $q$ is
the mass ratio, found from the RV fit with a circumbinary perturbation. With
these values as the center-of-mass (COM) RVs of the binary and RVs of the
tertiary component, we can treat the long-period outer orbit as an SB2
(Fig.1), and independently look for its parameters. The final values of
$P_{3}$, $K_{3}$, $e_{3}$, etc., actually come from such fits.

Figure 1: Phased RV profiles of BD44, for the inner binary (top) and the whole
system (COM of the binary vs the tertiary). The green lines are the best v2fit
models. Though the tertiary has a good fit, the residuals from the binary fit
contribute to the final uncertainties in the tertiary parameters.

###  3.2 Broadening Functions

Broadening function (BF) is a representation of the spectral profiles in
velocity space. BF contains the signature of the RV shifts of different lines
and also intrinsic stellar effects like rotational broadening, spots,
pulsations, etc (Rucinski, 1999a). The BF was calculated using the algorithm
described in Rucinski (1999b). We modified a single-order BF code, bf-
rvplotter333https://github.com/mrawls/BF-rvplotter, to calculate multi-order
BF and also fit the function with multiple Gaussian or rotation functions. The
BF was calculated in a wavelength range of 5050-5600 Å. We used a synthetic
solar-type spectrum with zero projected rotational velocity
($v\mathrm{sin}(i)$) as our template. The final BF generated was smoothed with
a Gaussian smoother of a 3 $\mathrm{km}\,\mathrm{s^{-1}}$ rolling-window.
Three clear peaks were visible in all epochs of spectra which implied that we
did not observe the system spectra during eclipses of any of the three stars.
The peaks were fitted with the rotational profile from Gray (2005),

$G(v)=A\left[c_{1}\sqrt{1-\left(\frac{v}{v_{\mathrm{max}}}\right)^{2}}+c_{2}\left(1-\left(\frac{v}{v_{\mathrm{max}}}\right)^{2}\right)\right]+lv+k$
(2)

where $A$ is the area under the profile and $v_{\mathrm{max}}$ is the maximum
velocity shift that occurs at the equator. $c_{1}$ and $c_{2}$ are constants
that are a function of limb darkening themselves while $l$ and $k$ are
correction factors to the BF ”continuum". The fits revealed distortions (sharp
kinks) of the BF from the ideal rotational profiles (Fig.2). These are most
likely due to spots, and this gives us a qualitative idea about the relative
number of spots on the stars. These fits for different epochs were used to (i)
calculate light contributions (the parameter $A$) from different companions,
and (ii) make an initial estimate of $v\mathrm{sin}(i)$ (the parameter
$v_{\mathrm{max}}$) for spectral analysis. We also used BF to calculate RVs
and found them consistent with the TODCOR RVs. Therefore, for the sake of
consistency and familiarity, we used RVs from TODCOR.

Figure 2: Broadening functions for all the three components of BD44 (top) and
KIC65 (bottom). The green and orange dashed line represents best-fit
rotational profile for BD44 and KIC65 respectively. The sharp kinks or
distortions on the BF peaks are most likely due to spots.

###  3.3 Spectral Disentangling

A detailed study of stellar evolution needs a model-independent estimate of
stellar metallicity. To estimate atmospheric parameters and abundances, we
need spectra for all three stars in the CHT. We use the technique of Spectra
Disentangling (hereafter: spd; Simon & Sturm, 1994; Hadrava, 1995), for
separating individual spectra of the component stars from the composite
spectra. This method, though, takes in the assumption that the line profiles
are not intrinsically variable. This would mean that we should consider only
the out-of-eclipse spectra so as to avoid such variability during the eclipse
(e.g., Rossiter-McLaughlin effect). One of the advantages of this method is
that it can detect faint companions, i.e., $<3\%$ light contribution (Holmgren
et al., 1999; Mayer et al., 2013). This, though, requires good phase coverage
and high signal-to-noise ratio.

For our purpose, we used a python-based
wrapper444https://github.com/ayushmoharana/fd3_initiator around the
disentangling code FDbinary (Ilijic et al., 2004), which can disentangle up to
three components. The wrapper takes two inputs: (i) an estimate of orbital
parameters and (ii) RV corrections and light-ratios at each epoch of spectra
used from disentangling. We used the solution from RV fitting as starting
values in our optimisation. The light-ratios at every epoch were calculated
from the BFs. To make the computation easy and avoid any wavelength
dependency, we divided the total spectral range into four sections with
overlapping regions. The final disentangled spectra were stitched after
normalisation, and removal of the edges of the segmented-disentangled spectra.
The overlapping regions acted as check-points for normalisation as they helped
us choose the normalisation function which gave the same line-depths for a
particular, overlapping spectral line. The errors of the disentangled spectra
were taken as the sum of errors calculated from SNR, and flux-scaled residuals
from the disentangled routine.

###  3.4 Light curve fitting

We use the version 4 of phoebe 2 code555http://phoebe-project.org/ (Prša et
al., 2016; Horvat et al., 2018; Jones et al., 2020; Conroy et al., 2020) for
our LC modelling. phoebe 2 models eclipsing binaries (or single stars) by
discretising the surface of each star. It also distorts the stellar surfaces
according to their Roche potentials. The key feature that made us choose this
code is its ability to model spots and also solve the inverse problem with
spot parameters free for optimisation.

The first look on the LCs of BD44 and KIC65 suggests the presence of cold
spots in the stars. Comparing the light-curve of the different sectors reveals
that the spots are time-evolving. We, therefore, approach our modelling by
dividing the light-curve into segments with relatively-stable spot signatures.
While we model all available sectors for BD44 (Fig.3), we only model only
Sector-41 for KIC65 (Fig.4). Further, the distortions on the BF suggest the
secondary to be more active than the primary for both stars. Therefore, in our
models, we assume more spot(s) on the secondary than the primary. The third-
light ($l_{3}$) is expected in a triple-lined CHT unless the signatures are
removed with detrending methods. The $l_{3}$ is highly degenerate with the
inclination ($i$) too. Estimating $l_{3}$ for our systems, from LC, adds
another level of complexity with the cold spots affecting the depths of the
eclipse. We assume that the values of light-fractions obtained from BF of
optical spectra are similar to the light-fraction of the components in the
TESS band. Therefore, we start with an initial $l_{3}$ equal to the flux-
fraction of the tertiary from BFs.

Considering the above assumptions and fixing parameters obtained from RV
fitting, i.e., mass-ratio ($q$), semi-major axis ($a$), and period of the
binary($P_{b}$), we used lc_geometry estimator in phoebe 2 for initial
parameter estimates. We then added spots one by one and minimising overall
trends in the residuals manually. We added 1 spot on the secondary of KIC65
and stopped at three spots (2 on secondary, 1 on primary) for BD44. More spots
would have compromised our computational resources. We then optimised for
stellar parameters using Nelder-Mead (Nelder & Mead, 1965) optimisation module
in phoebe 2. The optimised parameters include radii of primary and secondary
($R_{\mathrm{Aa}}$ and $R_{\mathrm{Ab}}$ respectively), time of super-
conjugation ($T_{0}$), ratio of secondary temperature to the primary
temperature ($T_{\mathrm{ratio}}$), $P_{A}$, $l_{3}$, $i_{A}$, and passband
luminosity of the primary ($L_{\mathrm{Aa}}$). Since the RV fitting didn’t
show any substantial eccentricity for the binary ($e_{A}$), we kept it fixed
at zero which also reduced the computational cost. After minimising the
residuals below 1% of the total flux, we optimised for all spot parameters
(co-latitude666Co-latitude is measured along the spin axis with the North pole
as $0^{\circ}$. $c^{\mathrm{spot}}$, longitude $l^{\mathrm{spot}}$, relative-
temperature $T^{\mathrm{spot}}$, and radius $r^{\mathrm{spot}}$) for all the
three spots. We then randomly optimised, combinations of spot parameters and
stellar parameters, to check the robustness of the optimisation.

We calculated the errors through the Monte Carlo Markov Chain (MCMC) sampling,
implemented in emcee (Foreman-Mackey et al., 2013; Foreman-Mackey et al.,
2019) and available as a sampler in phoebe 2. We decided on 8 parameters to be
sampled including the relative temperature of the biggest spot in both the
systems ($T^{\mathrm{spot}}_{Ab}$). Due to the lack of computational
resources, our initial sampling consisted of 40 walkers and was sampled till
we got stable chains for 1000 iterations. Then we started from this sample-
space and ran another sampling for 3500 stable chains with 80 walkers. To
further check the convergence, we use auto-correlation plots (Box & Jenkins,
1976) with Bartlett’s formula (Francq & Zakoïan, 2009) to find the minimum
limit for auto-correlated chains. The uncertainties represent 68.27%
confidence interval. To validate our results, we translated the errors on the
light curves and also the residuals (Fig. 3 and 4).

Figure 3: Phased TESS LC (for stable segments from all sectors) for BD44. The
green lines represent best-fit phoebe 2 model. The shaded areas are translated
from the adopted errors of the MCMC solution from Sector-16. Figure 4: Phased
TESS LC for KIC65. The orange lines represent best-fit phoebe 2 model. The
shaded areas are translated from the adopted errors of the MCMC solution from
Sector-41. Figure 5: Corner plot for MCMC parameters from LC fitting of BD44.
The contours on the maps correspond to 1$\sigma$, 2$\sigma$, and 3$\sigma$
errors on the distribution. The effect of spots is visible in the distorted
maps of the eccentricity parameter-space. Figure 6: Corner plot for MCMC
parameters from LC fitting of KIC65. The contours on the maps correspond to
1$\sigma$, 2$\sigma$, and 3$\sigma$ errors on the distribution. The effect of
spots is visible in the distorted maps of the eccentricity parameter-space.
Also note that the distortion of the radii parameter-space is smaller than
that of BD44.

###  3.5 Spectroscopic Analysis

We used iSpec (Blanco-Cuaresma et al., 2014; Blanco-Cuaresma, 2019) for
spectroscopic analysis. Before the analysis, we prepared the spectra by
correcting the RV offsets that was present in the spd products. We fitted and
corrected the continuum for the spd spectra using a spline function of degree
3, consisting of 30-150 splines, depending on the component. We used the
synthetic spectral fitting (SSF) method where iSpec generates synthetic
spectra on the go for the fitting. We used line masks for specific regions,
generated by the spectral synthesis code
spectrum777http://www.appstate.edu/~grayro/spectrum/spectrum.html to calculate
the $\chi^{2}$. The regions were generated on the line lists from Gaia-ESO
Survey (GES; Gilmore et al. 2012; Randich et al. 2013) version 5.0 which
covers the wavelength range of 420nm-920nm. There are two sets of GES line
lists available in iSpec, which we use: (i) list for abundances estimates
($\mathrm{LL_{a}}$) and (ii) list for parameter estimate ($\mathrm{LL_{p}}$).
The spectral synthesis was again done using spectrum with model atmospheres
from Gustafsson et al. (2008). The solar abundances were chosen from Grevesse
et al. (2007). In our initial test runs, with all parameters free, we found
that the $v\sin{i}$ values were consistent with the values obtained from BFs.
Therefore in our further runs, we fixed the respective $v\sin{i}$ for the
components. Using $\mathrm{LL_{a}}$, we did a search for metallicity, [M/H],
on all three stars of the two CHT. Since the [M/H] of the stars of the
respective CHT were consistent within the uncertainty range, we averaged the
[M/H] estimates for the system and fixed it for further runs. In the final
runs, we fixed the $\log(g)$ (obtained from LC models) for primary and
secondary. The limb darkening coefficient (values were taken from Claret &
Bloemen 2011) and resolution were kept fixed while marco-turbulence velocity
was calculated automatically from an empirical relation established by GES and
built in the code. For the estimation of temperatures (and $\log(g)$ for the
tertiary), we selected only the part of the spectra with SNR of more than 40
(more than 18 for the tertiary). We chose $\mathrm{LL_{p}}$ for this fitting.
The best fit synthetic spectra for all the stars in the study are shown in
Fig.7. We also calculated the radii $R$ (in R⊙) of the tertiary stars using
$\log(g)$ (in dex) from spectra, by applying the formula:

$R=A_{c}\sqrt{\frac{M}{10^{\log(g)}}},$ (3)

where $M$ is the mass of the tertiary calculated from LC and RV fitting (in
M⊙) and $A_{c}\equiv\sqrt{GM_{\odot}}/R_{\odot}(=168.589888477)$ is a constant
necessary for transformation to solar units.

We calculated the $\alpha$-enhancement with $\mathrm{LL_{a}}$ and fixed the
rest of the parameters as given in Table.1. Using this setup, the
abundances888The abundances were obtained in the 12-scale as $A(X)$, where,
$A(X)=\log\left(\frac{n_{X}}{n_{H}}\right)+12$, in which $n_{X}$ and $n_{H}$
are the number of atoms of the element X and of hydrogen, respectively. were
calculated using the SSF method but with free abundance and [M/H] for a
particular element. We did not consider the abundances of the elements where
we got large errors and/or where [M/H] was out of the pre-calculated bounds.

Figure 7: A section of the disentangled spectra and the best-fit synthetic
template of all the components of BD44 (top) and KIC65 (bottom). The fit was
obtained through the SSF method available in iSpec. While the SNR of the
primary and secondary vary between 50-60, the tertiary component has a SNR
$\sim$18\. The section contains one of the wings of the Mg triplet and several
Fe lines.

###  3.6 Isochrone Fitting

The age of each system was estimated with a grid of isochrones generated using
a dedicated web interface,999http://waps.cfa.harvard.edu/MIST/ based on the
Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al., 2011;
Paxton et al., 2013, 2015, 2018), and developed as part of the MESA Isochrones
and Stellar Tracks project (MIST v1.2; Choi et al., 2016; Dotter, 2016). The
grid was prepared for iron abundance, [Fe/H]101010It is reasonable to assume
that without significant deviations from solar amounts of $\alpha$-elements,
the iron abundance [Fe/H] sufficiently approximates the amount of metals
[M/H]., values from -4.0 to 0.50 dex with 0.05 dex steps, as well as for ages
108.6 to 1010.2 Gyr, in logarithmic scale, every $\log(\mathrm{age})$=0.01.

On each isochrone we were looking for a triplet of points that simultaneously
best reproduce the observed values of any selected parameters from the
following: masses $M_{1,2,3}$, radii $R_{1,2,3}$, and effective temperatures
$T^{\rm eff}_{\rm 1,2,3}$ of three components, flux ratio of the inner
binary’s components $l_{2}/l_{1}$ in the given band, as well as the [Fe/H],
and distance $d$. The method of simultaneously obtaining distances and
reddening $E(B-V)$ reproduced by a given model (triplet of points on an
isochrone) is described in Appendix A, as well as in Hełminiak et al. (2021).
In the recent 3rd Gaia Data Release (GDR3; Gaia Collaboration et al., 2022),
solutions from the Part 1. Main Source have high value of the RUWE parameter
($\sim$2.6), and better solutions, with significantly different distances, are
presented in Part 3. Non-single stars. We used the latter ones as the
constraints in our isochrone fitting.

It is worth noting that the ispec value of [M/H] was also used as a
constraint, and the “best-fitting” [Fe/H] was searched for in the fitting
process. For this reason, the values of [Fe/H] (assumed equivalent to [M/H],
since [$\alpha$/Fe]=0) may not be the same as [M/H] found from spectra. In the
literature, isochrone or evolutionary track fitting is often made under the
assumption of fixed [M/H]. We find this approach incorrect.

###  3.7 Numerical Integration of Orbital Dynamics

The accurate orbital parameters obtained enable us to probe the significant
dynamical changes in these compact systems. The major orbital parameters that
drive the dynamics are usually the masses, semi-major axes or periods, and all
inclinations, including mutual inclination ($i_{m}$). While we get almost all
orbital parameters from combined LC and RV analysis, we still lack the
information about the longitude of ascending nodes ($\Omega_{\mathrm{A}}$ and
$\Omega_{\mathrm{AB}}$ for inner and outer orbit respectively) and $i_{m}$. An
estimate of the range of the $i_{m}$ can be calculated from the constraints
arising due to geometry, using $i_{\mathrm{A}}$ and $i_{\mathrm{AB}}$.
Simplifying the calculations from Gronchi & Tommei (2007) we get,

$\cos{i_{m}}=\cos{(\Omega_{\mathrm{A}}-\Omega_{\mathrm{AB}})}\times\sin{i_{\mathrm{A}}}\sin{i_{\mathrm{AB}}}+\cos{i_{\mathrm{A}}}\cos{i_{\mathrm{AB}}}$
(4)

Since the value of $\cos{(\Omega_{\mathrm{A}}-\Omega_{\mathrm{AB}})}$ can vary
between -1 and 1, we estimate the range of $\cos{i_{m}}$ to be,

$\displaystyle\cos{i_{m}}\geq\cos{i_{\mathrm{A}}}\cos{i_{\mathrm{AB}}}-\sin{i_{\mathrm{A}}}\sin{i_{\mathrm{AB}}}$
(5)
$\displaystyle\cos{i_{m}}\leq\cos{i_{\mathrm{A}}}\cos{i_{\mathrm{AB}}}+\sin{i_{\mathrm{A}}}\sin{i_{\mathrm{AB}}}$
(6)

Using trigonometric identities, we get the range of $i_{m}$ to be,

$i_{\mathrm{A}}-i_{\mathrm{AB}}\leq i_{m}\leq i_{\mathrm{A}}+i_{\mathrm{AB}}$
(7)

This is the solution when the cosine of $i_{m}$ is positive (say configuration
A). There exists another set of solution for a negative cosine configuration
(say configuration B),

$180-(i_{\mathrm{A}}+i_{\mathrm{AB}})\leq i_{m}\leq
180-(i_{\mathrm{A}}-i_{\mathrm{AB}})$ (8)

This gives us a range of $0^{\circ}-169.4^{\circ}$ (for Config.A) or
$10.6^{\circ}-180^{\circ}$ (Config.B) for the $i_{m}$ of KIC65. While for
BD44, we get a range of $3.98^{\circ}-156.27^{\circ}$ (Config.A) or
$23.73^{\circ}-176.02^{\circ}$ (Config.B). To further constrain this range, we
rule out the unrealistic $i_{m}$ by looking at average $i_{\mathrm{A}}$
variations in the numerical integration of orbital parameters, and comparing
with the $i_{\mathrm{A}}$ from the observations. For our work, we use
rebound111111https://github.com/hannorein/rebound, an open-source collisional
N-body code (Rein & Liu, 2012). rebound can also be used to simulate
collision-less problems such as the three-body hierarchical orbit in our case.
We use the symplectic integrator whfast which is designed for long-term
integration of gravitational orbits (Rein & Tamayo, 2015). whfast uses mixed-
variables (Jacobi and Cartesian) and also a symplectic corrector which ensures
accurate and fast integration of the N-body dynamical equations. The general
setup is defined by masses and orbital parameters obtained with our
observations. The values of inclinations, $\omega$ and other orientation
parameters are also taken from our observations (Table.1). Further, we use the
reboundx (Tamayo et al., 2020) library for adding in tidal forces and
dissipation to this setup. The implementation uses a constant time-lag model
from (Hut, 1981) to raise tides on the larger (and more massive) stars. The
constant time-lag parameter ($\tau$) is given by,

$\centering\tau=\frac{2R^{3}}{GMt_{f}}\@add@centering$ (9)

where $R$, $M$ are the radius and mass of the body with tides (BD44Ab and
KIC65Aa in our case) and $t_{f}$ is convective friction time which is assumed
as 1yr as described in reboundx setup121212We followed the process of adding
tides as explained in:
https://github.com/dtamayo/reboundx/blob/master/ipython_examples/TidesConstantTimeLag.ipynb.
$G$ is the gravitational constant. The other inputs from observations include
the rotation frequency of the secondary and its radius. The tidal Love number
of degree 2 (tctlk2) is assumed to be 0.01. Using this setup, we simulate the
systems for a time equal to the time between the observation of the first and
last LC of the corresponding systems, i.e., 2.7yrs for BD44 and 6.9yrs for
KIC65. We track the inclination changes of the inner binary for different
values of $\Omega_{\mathrm{AB}}$ (and assuming
$\Omega_{\mathrm{A}}=0^{\circ}$) and then compare it with the observations.
A long-term look at such close systems ($P_{A}\sim 3d$) would require a
coupled stellar evolution and dynamical evolution treatment. For this, we use
the parameter-interpolation module in reboundx. We update the radius and
masses in the numerical simulation from MIST evolutionary tracks for the tidal
stars (BD44Ab and KIC65Aa) in the inner binary. We start with an integration
timestep of 1/30 times the inner binary period but update the integration
timestep, if the orbit shrinks, so that we can track close encounters or
collisions. With this setup, we simulate our systems for 600 Myrs.

##  4 Results and Discussion

###  4.1 Physical parameters of the Binary

The binary components of the two systems are found to be of similar orbital
configurations. This gives us an interesting case of comparing the effect of
the tertiary on the binary system. The inner mass-ratio of BD44 is 0.931 ( or
1/1.074) while for KIC65 it is 0.9383. The corresponding inner orbit periods
of BD44 and KIC65 are 3.4726217d and 3.4205977d, respectively. The differences
in the binaries of the two systems appear when we look at the LC solutions.
The radius of the most massive star in BD44 is larger
($R_{\mathrm{Ab}}=2.4387_{-0.0048}^{+0.0035}R_{\odot}$) but it corresponds to
the secondary (shallower) eclipse in the LC. This is also evident from the
light-fractions obtained from the BF and also the fast rotation as reflected
in the high $v\mathrm{sin}(i)$. Going by the convention, we will address this
star as the secondary (BD44Ab; thus the secondary-to-primary mass ratio is
$>$1). The primary of BD44 (BD44Aa) is inflated compared to the solar radius
($R_{Aa}=1.3118_{-0.0053}^{+0.0056}R_{\odot}$). Both the stars in KIC65 are
very much solar-like ($R_{Aa}=1.0719_{-0.0020}^{+0.0019}R_{\odot}$ and
$R_{Ab}=1.0101_{-0.0014}^{+0.0017}R_{\odot}$). A quick look at the MCMC corner
plots (Fig. 5) tells us that the uncertainty on radii measurements is mostly
affected by degeneracy in $L_{\mathrm{Aa}}$ which is itself degenerate with
$l_{3}$. The $e_{\mathrm{A}}$ is adopted to be zero for both systems. But the
MCMC maps show shifts of $2\sigma$ and $3\sigma$ from zero, for BD44 (Fig.5)
and KIC65 (Fig.6) respectively. This small eccentricity can possibly be
induced in the LC due to spots or be due to perturbations from the tertiary
star. Unfortunately with current observations, it is difficult to decouple
these effects.

The temperatures of the BD44 binary stars are lower
($T^{\mathrm{eff}}_{\mathrm{Aa}}=5822\pm 202K$ and
$T^{\mathrm{eff}}_{\mathrm{Ab}}=5449\pm 100K$ ) compared to those of KIC65
($T^{\mathrm{eff}}_{\mathrm{Aa}}=6490\pm 129K$ and
$T^{\mathrm{eff}}_{\mathrm{Ab}}=6397\pm 123K$). While the similar temperatures
of KIC65Aa and KIC65Ab can explain the similar eclipse depths, we expect a
temperature ratio of 0.772 (compared to 0.936 from the spectral analysis) from
the LC fitting of BD44. This discrepancy can arise due to the cold spots in
BD44, which affect the spectroscopic temperature measurement.

The stars in BD44 have quite different $v\mathrm{sin}(i)$ with the BD44Aa
having a value of 23.21$\,\mathrm{km\,s}^{-1}$ compared to
39.60$\,\mathrm{km\,s}^{-1}$ for BD44Ab. KIC65 has similar $v\mathrm{sin}(i)$
for the stars in the inner binary as expected from stars with similar radii.
But when comparing with calculated synchronised velocities, we find that the
observed velocities are larger than expected. The $v\mathrm{sin}(i)$ of the
tertiary of KIC65 and BD44 are similar.

Both KIC65 and BD44 are metal-poor systems having [M/H] $-0.28\pm 0.22$ and
$-0.24\pm 0.21$ respectively. While most of the stars are found to be
$\alpha$-enhanced systems, BD44Ab and BD44B were the only stars with negative
$\alpha$. We calculated the abundances of some elements whose lines showed up
on the spectra. We compare the abundances (see Tables 2 and 3) for all three
stars, in both the systems, in Fig.8. While all the stars in KIC65 have
similar abundances (within error bars), BD44Ab and BD44B have a rise in
$A(Si)$ and a dip in $A(Mg)$ compared to BD44Aa (Fig.8). The total set of all
parameters is given in Table. 1 for comparison.

Table 1: All adopted parameters for BD44 and KIC65 (except $i_{m}$).
Unsymmetrical errors correspond to uncertainties estimated using MCMC
sampling.

| BD +44 2258 | KIC 06525196
---|---|---
Orbital Parameters
| Aa–Ab | A–B | Aa–Ab | A–B
$t_{0}$ [BJD - 2450000] | $8740.8268426$ | $2421.6563579$
$P$ [days] | $3.4726217\pm 0.0000015$ | $254.84\pm 0.05$ | $3.42059774\pm 0.00000014$ | $418.0\pm 0.4$
$a$ [R⊙] | $12.06\pm 0.12$ | $240.2\pm 3.3$ | $12.053\pm 0.018$ | $330.8\pm 2.1$
$e$ | $0.000^{a}$ | $0.598\pm 0.002$ | $0.000^{a}$ | $0.301\pm 0.003$
$i$ [deg] | $80.1271_{-0.0057}^{+0.0084}$ | $76.15_{-1.61}^{+1.91}$ | $84.7083_{-0.0063}^{+0.0051}$ | $84.7_{-2.6}^{+5.3}$
$\omega$ [deg] | $-$ | $203.0\pm 0.2$ | $-$ | $276\pm 1$
$q$ | $1.074\pm 0.021$ | $0.465\pm 0.019$ | $0.9383\pm 0.0026$ | $0.3875\pm 0.0045$
$K_{1}$ [km s-1] | $89.66\pm 0.98$ | $18.30\pm 0.75$ | $85.96\pm 0.12$ | $11.67\pm 0.12$
$K_{2}$ [km s-1] | $83.47\pm 1.37$ | $39.39\pm 0.22$ | $91.62\pm 0.12$ | $30.15\pm 0.17$
Stellar and atmospheric parameters
| Aa | Ab | B | Aa | Ab | B
Flux fraction (from Spectroscopy) | $0.3758\pm 0.0268$ | $0.5139\pm 0.0285$ | $0.1104\pm 0.0088$ | $0.4851\pm 0.0225$ | $0.4010\pm 0.0154$ | $0.1139\pm 0.0181$
Flux fraction (from Photometry) | $0.4021_{-0.0030}^{+0.0039}$ | $0.5101_{-0.0031}^{+0.0021}$ | $0.0878_{-0.0015}^{+0.0016}$ | $0.4990_{-0.0084}^{+0.0101}$ | $0.4060_{-0.0068}^{+0.0051}$ | $0.0950_{-0.0034}^{+0.0033}$
$M$ [M⊙] | $0.941\pm 0.033$ | $1.011\pm 0.029$ | $0.907\pm 0.065$ | $1.0351\pm 0.0055$ | $0.9712\pm 0.0039$ | $0.777\pm 0.012$
$R$ [R⊙] | $1.3118_{-0.0053}^{+0.0056}$ | $2.4387_{-0.0048}^{+0.0035}$ | $1.67_{-0.59}^{+0.70}$ | $1.0719_{-0.0019}^{+0.0020}$ | $1.0101_{-0.0017}^{+0.0014}$ | $0.74_{-0.38}^{+0.80}$
$T_{\mathrm{eff}}$ [K] | $5822\pm 202$ | $5449\pm 100$ | $5261\pm 181$ | $6490\pm 129$ | $6397\pm 123$ | $5393\pm 319$
$\log(g)$ [dex] | $4.19^{b}$ | $3.69^{b}$ | $3.95\pm 0.39$ | $4.39^{b}$ | $4.41^{b}$ | $4.58\pm 0.63$
$v_{\mathrm{mic}}$ [km s-1] | $2.36\pm 0.75$ | $1.78^{c}$ | $1.19\pm 0.5$ | $1.39\pm 0.54$ | $1.18\pm 0.53$ | $1.5^{c}$
$v_{\mathrm{mac}}$ [km s-1]c | $4.69$ | $3.76$ | $3.43$ | $9.28$ | $8.43$ | $5.94$
$v\mathrm{sin}(i)$ [km s-1] | $23.21\pm 2.19$ | $39.60\pm 3.26$ | $5.34\pm 1.41$ | $20.69\pm 1.45$ | $18.39\pm 1.15$ | $5.23\pm 1.34$
$\alpha$ [dex] | $0.21\pm 0.12$ | $-0.22\pm 0.08$ | $-0.21\pm 0.11$ | $0.07\pm 0.11$ | $0.16\pm 0.11$ | $0.29\pm 0.24$
System parameters
$\log$(age) [dex] | $9.89_{-0.05}^{+0.03}$ | $9.49\pm 0.06$
$\mathrm{[M/H]}_{\mathrm{iSpec}}$ [dex] | $-0.24\pm 0.21$ | $-0.28\pm 0.22$
$\mathrm{[Fe/H]}_{\mathrm{isoc}}$ [dex] | $-0.40^{+0.15}_{-0.10}$ | $-0.45^{+0.15}_{-0.35}$
$E(B-V)^{d}$ [mag] | $0.176^{+0.019}_{-0.049}$ | $0.053^{+0.040}_{-0.030}$
Distanced [pc] | $194.5_{-10.3}^{+3.3}$ | $220.01^{+2.26}_{-1.68}$

a Fixed while optimisation. b Fixed from LC fitting solutions. c Obtained from
empirical tables. d Based on isochrone fitting.

Table 2: Abundances of individual elements of the components of BD44. Solar abundances are presented for comparison (Asplund et al. 2009). Elements | Primary | Secondary | Tertiary | Solar
---|---|---|---|---
12Mg | 7.59 $\pm$ 0.03 | 6.94 $\pm$ 0.06 | 7.16 $\pm$ 0.03 | 7.60 $\pm$ 0.04
14Si | 7.20 $\pm$ 0.39 | 8.26 $\pm$ 0.07 | 8.63 $\pm$ 0.06 | 7.51 $\pm$ 0.03
22Ti | 4.90 $\pm$ 0.06 | 5.25 $\pm$ 0.14 | 5.01 $\pm$ 0.03 | 4.95 $\pm$ 0.05
24Cr | 5.21 $\pm$ 0.09 | 5.31 $\pm$ 0.31 | 5.49 $\pm$ 0.03 | 5.64 $\pm$ 0.04
25Mn | 5.13 $\pm$ 0.10 | 5.39 $\pm$ 0.19 | 5.20 $\pm$ 0.03 | 5.43 $\pm$ 0.05
26Fe | 7.18 $\pm$ 0.03 | 7.20 $\pm$ 0.05 | 7.30 $\pm$ 0.01 | 7.50 $\pm$ 0.04
28Ni | 5.77 $\pm$ 0.10 | 6.19 $\pm$ 0.20 | 6.07 $\pm$ 0.05 | 6.22 $\pm$ 0.04
Table 3: Abundances of individual elements of the components of KIC65. Elements | Primary | Secondary | Tertiary | Solar
---|---|---|---|---
12Mg | 7.36 $\pm$ 0.11 | 7.33 $\pm$ 0.12 | 7.41 $\pm$ 0.17 | 7.60 $\pm$ 0.04
22Ti | 4.67 $\pm$ 0.28 | 4.67 $\pm$ 0.30 | 5.16 $\pm$ 0.49 | 4.95 $\pm$ 0.05
24Cr | 5.33 $\pm$ 0.30 | 5.44 $\pm$ 0.29 | 5.17 $\pm$ 0.43 | 5.64 $\pm$ 0.04
26Fe | 7.12 $\pm$ 0.14 | 7.21 $\pm$ 0.13 | 7.13 $\pm$ 0.41 | 7.50 $\pm$ 0.04

Figure 8: Abundance variation for all three stars in BD44 (top) and KIC65
(bottom). The plotted abundances are relative to solar abundances, which is
denoted by the dotted line.

###  4.2 Physical parameters of the Tertiary

The tertiary stars orbit the inner binary with the periods of $254.84\pm 0.05$
and $418.0\pm 0.4$ d for BD44 and KIC65, respectively. They are different in
both mass and radius. The mass of the tertiary in BD44 is close to solar
($0.907\pm 0.065$ M⊙) while tertiary of KIC65 is less massive ($0.777\pm
0.012$ M⊙). The estimated radii have large errors, but a comprehensive look at
all the signatures suggests that the tertiary in BD44 has an inflated
atmosphere and therefore a radius larger than one expected for a main-sequence
star of this mass. The tertiary of KIC65 is most likely to have a radius of
0.74 $R_{\odot}$. The stars themselves are orbiting with different periods
around the inner binary system. The tertiary of BD44 is in an orbit with
higher eccentricity ($e_{\mathrm{AB}}=0.598$) than KIC65 ($e_{\mathrm{AB}}\sim
0.3$; Fig.9). These different configurations will mostly affect the timescale
of secular perturbations which depend on $e_{\mathrm{AB}}$ (Ford et al.,
2000).

Figure 9: Orbital path of all the stars of BD44 (top) and KIC65 (bottom). The
orbits are viewed perpendicular to the orbital plane of the inner binary (for
$i_{m}=4^{\circ}$) and are integrated over one outer-orbit period.

###  4.3 Possible Mutual Inclinations

Short-term numerical integration gave us an estimate of inclination changes of
the inner binary ($\Delta i_{\mathrm{A}}$) for different values of
$\Omega_{\mathrm{AB}}$. We then obtained observed $\Delta i_{\mathrm{A}}$ by
using inclination values obtained by Hełminiak et al. (2017) for KIC65,
subtracted from the values obtained in this work. While for BD44 we used the
inclinations observed in Sector-16 and Sector-49 of TESS LCs (Table.4). This
gave us possible values of $\Omega_{\mathrm{AB}}$ for the observed $\Delta
i_{\mathrm{A}}$ (Fig.10). We then used Eq.(4) to translate the possible
$\Omega_{\mathrm{AB}}$ values to possible $i_{m}$ values for all possible
configurations (Table.5).

Table 4: Variation of $i_{\mathrm{A}}$ over time. System | Initial $i_{\mathrm{A}}$ [deg] | Final $i_{\mathrm{A}}$ [deg] | Time (yrs)
---|---|---|---
BD44 | $80.1271\pm 0.008$ | $80.1521\pm 0.017$ | 2.7
KIC65 | $85.15\pm 0.34$ | $84.7083\pm 0.005$ | 6.9

Figure 10: Variation of $|\Delta i_{\mathrm{A}}|$ from numerical simulations for different values of $\Omega_{\mathrm{AB}}$, for BD44 (top) and KIC65 (bottom). The timescale of the simulation corresponds to the length of the available LCs. The grey lines show the observed $|\Delta i_{\mathrm{A}}|$ with the shaded regions representing the errors in measurements. The overlapping regions represent the possible values of $\Omega_{\mathrm{AB}}$. The lower panel of the top figure represents a zoomed-in view of the observed variations. Table 5: Estimates of mutual inclination for possible values of $\Omega_{\mathrm{AB}}$. System | $\Omega_{\mathrm{AB}}$ | $i_{m}$
---|---|---
|  | Config.A | Config.B
BD44 | $359.46\pm 0.76$ | $4.03\pm 0.10$ | $175.97\pm 0.01$
$266.97\pm 0.79$ | $90.56\pm 0.76$ | $89.45\pm 0.76$
$180.23\pm 0.83$ | $156.27\pm 0.01$ | $23.73\pm 0.01$
$92.16\pm 0.80$ | $89.71\pm 0.77$ | $90.29\pm 0.77$
KIC65 | $353.80\pm 4.85$ | $6.18\pm 4.83$ | $173.82\pm 4.83$
$269.41\pm 11.12$ | $90.09\pm 11.02$ | $89.91\pm 11.02$
$180.42\pm 11.18$ | $164.60\pm 0.30$ | $15.40\pm 0.30$
$90.67\pm 11.10$ | $90.16\pm 11.00$ | $89.83\pm 11.00$
| $6.19\pm 4.85$ | $6.17\pm 4.83$ | $173.83\pm 4.83$

###  4.4 Age and Evolution

Isochrone fitting puts the log(age) for BD44 between 9.84 and 9.92 (95%
confidence level), for metallicity range of -0.25 to -0.50 dex. The formally
best fit was found for log(age) = 9.89 (7.8 Gyr) and [Fe/H] = -0.40 dex. While
the primary of BD44, BD44Aa, is a main-sequence star, BD44Ab is a sub-giant.
With the large uncertainties in the parameters of BD44B, it is hard to
determine its evolutionary state. But the simultaneous mass-radius and mass-
temperature isochrone fit depicts it as a sub-giant star (Fig.11). The other
signatures of BD44B being a sub-giant are found in (i) large amplitude of BF
(Fig.2), and (ii) similar abundances (Fig.8) as that of the BD44Ab (which is a
sub-giant itself). However, with the available data, we can not completely
rule out the possibility that it is a main-sequence star, less massive and
smaller than the primary.

The isochrone-based, reddening-free distance was found to be 194.5 pc, which
is significantly larger than the GDR3 Part 3 value of 165.8$\pm$0.6 pc. The
distribution of acceptable models is very skewed, and none of the acceptable
models reached a distance lower than 184 pc. The tension is probably caused by
the tertiary, which was formally found to be less massive but seemingly more
evolved than the primary. Its parameters could probably be better determined
with additional observations around the outer orbit’s pericenter, where the
tertiary’s RVs reach their minimum. It should also be noted that the GDR3
solution is of worse quality than for KIC65.

Even though having the binary stars with masses similar to BD44, both the
primary and secondary of KIC65 are main-sequence stars along with its
tertiary. This immediately suggests that KIC65 is significantly younger than
BD44. The fitting procedure for KIC65 resulted in log(age) = 9.49 (3.1 Gyr)
and [Fe/H] = $-0.45$ dex, with the 95% confidence level ranges of 9.43 to
9.55, and $-0.80$ to $-0.30$ dex, respectively. The isochrone-based reddening-
free distance ($\sim$220$\pm$2 pc) is in excellent agreement with, and of
comparable precision to, the GDR3 solution for an astrometric binary model
(222.3$\pm$1.7 pc), even when it was not used as a constraint.

Figure 11: Mass-Radius (left column) and Mass-Temperature (right column)
isochrone fits for BD44 (top; in green) and KIC65 (bottom; in orange). The
shaded area represents the parameter-space corresponding to the errors of the
fitting. The respective measurements from the stars are plotted in grey. The
size of the points are relative to their radii.

###  4.5 Dynamical Evolution

Long-term evolution of KIC65 shows that the system is stable for 600 Myrs
more. But BD44 becomes a binary system within 550 Myrs due to the
collision/merger of the inner binary. This collision will be driven by
increasing tidal forces due to the increasing radius of the sub-giant BD44Ab
(Fig.12: upper-panel). The radius of the star will exceed the Roche-limit
(Eggleton, 1983) at around 450 Myr and will drive the merger process unless
the formation of a contact binary stabilises this system. This merger is
mostly due to the tides in the inner binary as a lack of tertiary companion
would have only delayed the merger by a few Myrs (Fig.12: lower-panel) for
most of the estimated $i_{m}$. But a $i_{m}$ near $90^{\circ}$ will make the
system merge faster ($<$400 Myr) than the lower values of $i_{m}$ (Fig.13).

Figure 12: Variation of the radius of BD44Ab and the semi-major axis of the
inner binary of BD44 as simulated in reboundx for $i_{m}=4^{\circ}$. The
radius variation is interpolated from MIST grids. The grey dotted line in the
top panel shows Roche-limit for the system. The dashed line in the lower panel
shows the change in the collision time if no tertiary star was present. Figure
13: The simulated collision time of the inner binary of BD44 for different
mutual inclinations. The green dot-dashed line shows the time for collision
without a tertiary companion. The vertical grey lines represent the possible
values of $i_{m}$. Solid lines represent possible $i_{m}$ for Config.A while
dashed lines are for Config.B.

###  4.6 Spot Evolution

The light curve of BD44 is highly varying over different sectors owing to the
migrating and evolving cold spots. The activity of BD44 is corroborated by its
ultraviolet and X-ray emissions. The distortions on the BFs of the stars in
BD44 indicate the secondary (BD44Ab) to have more spots. The occurrence of the
fast-evolving and migrating spots is common on sub-giant stars. The biggest
spot on the secondary is still visible in the newest Sector-49 of TESS
observations. This enabled us to study its migration. We used spot parameters
obtained from LC fitting using phoebe 2, from different sectors, to quantify
the migration (see Table 6). Using the values of longitudes and the mid-times
of the first primary eclipse ($T_{\mathrm{seg}}$) for each segment, we
calculated the rate of change of the longitude. We found that the spot moves
$\sim$ $0.696^{\circ}$ per orbital cycle of the inner binary. The spot moves
from a longitude of $0^{\circ}-180^{\circ}$ in 2.5 years. Extrapolating this,
we get a spot migration period of 5 years. The migration is most probably
caused by differential rotation because the spot in question is near the
poles, as seen in low $c^{\mathrm{spot}}$ values (Fig.14). Petit et al. (2004)
represents the differential rotation as,

$\Omega(c^{\mathrm{spot}})=\Omega_{\mathrm{eq}}-\Delta\Omega\cos^{2}(c^{\mathrm{spot}})$
(10)

where, $\Omega$ is the differential rotation as a function of co-latitude (or
latitude), $\Omega_{\mathrm{eq}}$ is the rotation at the equator, and
$\Delta\Omega$ is the difference in rotation rate between the pole and the
equator. Calculating an estimate of $\Delta\Omega$ in BD44Ab gives us
$\Delta\Omega=0.0044\pm 0.0004\,\mathrm{rad\,d^{-1}}$. This small differential
rotation has been seen in the K1 sub-giant primary in the RS CVn system HR
1099 (Petit et al., 2004) and also K-type main-sequence of V471 Tau (Hussain
et al., 2006) with $\Delta\Omega=0.0152\pm 0.0008\,\mathrm{rad\,d^{-1}}$ and
$\Delta\Omega=0.0016\pm 0.0060\,\mathrm{rad\,d^{-1}}$ respectively.

Table 6: Spot parameters for the coldest spot on the secondary of BD44. The parameters have been obtained for each sector using phoebe 2 modelling. Parameters | S16 | S22 | S49
---|---|---|---
$T_{\mathrm{seg}}$ [BJD-2457000] | 1740.8270 | 1900.5753 | 2640.2640
$r^{\mathrm{spot}}_{\mathrm{Ab}}$ [deg] | 42.846 | 37.311 | 39.611
$T^{\mathrm{spot}}_{\mathrm{Ab}}$ | 0.8988 | 0.9296 | 0.8500
$c^{\mathrm{spot}}_{\mathrm{Ab}}$ [deg] | 22.000 | 35.745 | 21.913
$l^{\mathrm{spot}}_{\mathrm{Ab}}$ [deg] | 0.000 | 29.975 | 180.233

Figure 14: Spot positions on the inner binary stars of BD44 (at the same
phase) for Sector-16 (top), Sector-22 (middle), and Sector-40 (bottom) as
obtained with phoebe 2 modelling. The stable-spot near the pole of the
secondary is the spot that affects the eclipse depths the most and is the one
that migrates with a timescale of $\sim$ 5 yr.

##  5 Conclusions

We obtained independent measurements of different parameters for two triple-
lined CHT. Using LC modelling, RV modelling, and spd followed by spectral
analysis, we obtained stellar, orbital and atmospheric parameters of all the
six stars in the two CHT. A multi-parameter isochrone fitting constrained the
ages of the two systems to be in the order of Gyrs. Isochrones, along with
abundances obtained from the disentangled spectra, helped us classify the
evolutionary state of the tertiary in the two systems. Furthermore, we
gathered the following information about the two systems:

  * •

KIC65: The period-ratio of the CHT is $\approx$ 122\. All stars in the system
are main-sequence stars. The system is a metal-poor one and is
$\alpha$-enhanced. The tertiary has 5 possible configurations of mutual
inclination with a possibility of near co-planar orbit. Due to the
comparatively smaller mass of the tertiary and a wider orbit, the system is
stable in the long-term for all values of mutual inclinations. The distance
estimated in our study is consistent with the distance obtained in Gaia-DR3.

  * •

BD44: This system is a relatively tighter CHT with a period-ratio $\approx$
73\. But still, this system is well above the dynamical instability131313This
limit is derived for prograde co-planar motion. limit as defined in Mardling &
Aarseth (2001). The system consists of one main-sequence star (almost at the
turn-off) and two sub-giant stars. The abundance patterns of the two sub-
giants are similar while the main-sequence differs in Mg and Si abundances.
The system has large and cold spots which affect the measurement of some of
our parameters. But using the spots we were able to calculate a differential
rotation in the sub-giant of the inner binary. This sub-giant component of the
inner pair also contributes significantly to tidal forces. Numerical
simulations with tidal interactions show that the inner binary will
collide/merge in a few hundred of Myrs due to the radius of the sub-giant
exceeding the Roche limit. This leaves behind a wide binary unless the
formation of a contact binary stabilises the system. The tertiary does hasten
this merger but the effects are drastic if the tertiary is orbiting in an
orbit perpendicular to the inner binary orbit.

Both the targets can benefit from further photometric and spectroscopic
observations which will improve the estimate of obtained parameters. The
photometric observations themselves will be quite crucial to check for
inclination variations of the inner binary and therefore will give better
constraints on the mutual inclination. This will be helpful in constraining
evolution scenarios of the CHT as well as star-formation scenarios of these
close triples. Nevertheless, this study shows that the use of parameters
obtained using independent observations is crucial for the realistic modelling
of CHT.

## Acknowledgements

The authors thank the referee for the invaluable comments and suggestions. The
authors thank Dr. Kyle Conroy and Dr. Andrej Prša for their valuable
suggestions and help with phoebe 2. The authors also acknowledge Dr. Hanno
Rein for his help in setting up rebound simulations. This work is funded by
the Polish National Science Centre (NCN) through grant 2021/41/N/ST9/02746.
A.M., F.M., T.P., and M.K. are supported by NCN through grant no.
2017/27/B/ST9/02727. F.M. gratefully acknowledges support from the NASA TESS
Guest Investigator grant 80NSSC22K0180 (PI A. Prša). This paper includes data
collected with the TESS mission, obtained from the Mikulski Archive for Space
Telescopes (MAST) data archive at the Space Telescope Science Institute
(STScI). Funding for the TESS mission is provided by the NASA Explorer
Program. STScI is operated by the Association of Universities for Research in
Astronomy, Inc., under NASA contract NAS 5–26555. This work also presents
results from the European Space Agency (ESA) space mission Gaia. Gaia data are
being processed by the Gaia Data Processing and Analysis Consortium (DPAC).
Funding for the DPAC is provided by national institutions, in particular, the
institutions participating in the Gaia MultiLateral Agreement (MLA).

## Data Availability

The TESS data used in this article are public and hosted on MAST. The data can
be accessed from http://dx.doi.org/10.17909/fp8v-5705. The spectroscopic data
will be made available upon request.

## References

  * Akashi & Soker (2017) Akashi M., Soker N., 2017, MNRAS, 469, 3296
  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  * Bianchi et al. (2011) Bianchi L., Herald J., Efremova B., Girardi L., Zabot A., Marigo P., Conti A., Shiao B., 2011, Ap&SS, 335, 161
  * Blanco-Cuaresma (2019) Blanco-Cuaresma S., 2019, MNRAS, 486, 2075
  * Blanco-Cuaresma et al. (2014) Blanco-Cuaresma S., Soubiran C., Heiter U., Jofré P., 2014, A&A, 569, A111
  * Borkovits (2022) Borkovits T., 2022, Galaxies, 10, 9
  * Borkovits et al. (2015) Borkovits T., Rappaport S., Hajdu T., Sztakovics J., 2015, MNRAS, 448, 946
  * Borkovits et al. (2016) Borkovits T., Hajdu T., Sztakovics J., Rappaport S., Levine A., Bíró I. B., Klagyivik P., 2016, MNRAS, 455, 4136
  * Borkovits et al. (2020) Borkovits T., Rappaport S. A., Hajdu T., Maxted P. F. L., Pál A., Forgács-Dajka E., Klagyivik P., Mitnyan T., 2020, MNRAS, 493, 5005
  * Borkovits et al. (2022) Borkovits T., Rappaport S. A., Toonen S., Moe M., Mitnyan T., Csányi I., 2022, Monthly Notices of the Royal Astronomical Society
  * Box & Jenkins (1976) Box G. E. P., Jenkins G. M., 1976, in Holden-Day Series in Time Series Analysis.
  * Cardelli et al. (1989) Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Claret & Bloemen (2011) Claret A., Bloemen S., 2011, A&A, 529, A75
  * Conroy et al. (2020) Conroy K. E., et al., 2020, ApJS, 250, 34
  * Dotter (2016) Dotter A., 2016, ApJS, 222, 8
  * Duchêne & Kraus (2013) Duchêne G., Kraus A., 2013, ARA&A, 51, 269
  * Eggleton (1983) Eggleton P. P., 1983, ApJ, 268, 368
  * Eggleton & Kiseleva-Eggleton (2001) Eggleton P. P., Kiseleva-Eggleton L., 2001, ApJ, 562, 1012
  * Eggleton & Verbunt (1986) Eggleton P. P., Verbunt F., 1986, MNRAS, 220, 13P
  * Eisner et al. (2022) Eisner N. L., et al., 2022, MNRAS, 511, 4710
  * Ford et al. (2000) Ford E. B., Kozinsky B., Rasio F. A., 2000, ApJ, 535, 385
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Foreman-Mackey et al. (2019) Foreman-Mackey D., et al., 2019, The Journal of Open Source Software, 4, 1864
  * Francq & Zakoïan (2009) Francq C., Zakoïan J.-M., 2009, Journal of Time Series Analysis, 30, 449
  * Fuller et al. (2013) Fuller J., Derekas A., Borkovits T., Huber D., Bedding T. R., Kiss L. L., 2013, MNRAS, 429, 2425
  * Gaia Collaboration et al. (2022) Gaia Collaboration et al., 2022, arXiv e-prints, p. arXiv:2208.00211
  * Gilmore et al. (2012) Gilmore G., et al., 2012, The Messenger, 147, 25
  * Gray (2005) Gray D. F., 2005, The Observation and Analysis of Stellar Photospheres. Cambridge University Press
  * Grevesse et al. (2007) Grevesse N., Asplund M., Sauval A. J., 2007, Space Sci. Rev., 130, 105
  * Gronchi & Tommei (2007) Gronchi G. F., Tommei G., 2007, Discrete and Continuous Dynamical Systems - B, 7, 755
  * Gustafsson et al. (2008) Gustafsson B., Edvardsson B., Eriksson K., Jørgensen U. G., Nordlund Å., Plez B., 2008, A&A, 486, 951
  * Hadrava (1995) Hadrava P., 1995, A&AS, 114, 393
  * Hamers & Dosopoulou (2019) Hamers A. S., Dosopoulou F., 2019, ApJ, 872, 119
  * Hełminiak et al. (2016) Hełminiak K. G., Ukita N., Kambe E., Kozłowski S. K., Sybilski P., Ratajczak M., Maehara H., Konacki M., 2016, MNRAS, 461, 2896
  * Hełminiak et al. (2017) Hełminiak K. G., et al., 2017, MNRAS, 468, 1726
  * Hełminiak et al. (2021) Hełminiak K. G., et al., 2021, MNRAS, 508, 5687
  * Holmgren et al. (1999) Holmgren D. E., Hadrava P., Harmanec P., Eenens P., Corral L. J., Yang S., Ak H., Bozić H., 1999, A&A, 345, 855
  * Horvat et al. (2018) Horvat M., Conroy K. E., Pablo H., Hambleton K. M., Kochoska A., Giammarco J., Prša A., 2018, ApJS, 237, 26
  * Howell et al. (2014) Howell S. B., et al., 2014, PASP, 126, 398
  * Hussain et al. (2006) Hussain G. A. J., Allende Prieto C., Saar S. H., Still M., 2006, MNRAS, 367, 1699
  * Hut (1981) Hut P., 1981, A&A, 99, 126
  * Ilijic et al. (2004) Ilijic S., Hensberge H., Pavlovski K., Freyhammer L. M., 2004, in Hilditch R. W., Hensberge H., Pavlovski K., eds, Astronomical Society of the Pacific Conference Series Vol. 318, Spectroscopically and Spatially Resolving the Components of the Close Binary Stars. pp 111–113
  * Izumiura (1999) Izumiura H., 1999, in Chen P. S., ed., Proc. 4th East Asian Meeting on Astronomy, Observational Astrophysics in Asia and its Future. Kunming Yunnan Observatory, p. 77
  * Jones et al. (2019) Jones D., Pejcha O., Corradi R. L. M., 2019, MNRAS, 489, 2195
  * Jones et al. (2020) Jones D., et al., 2020, ApJS, 247, 63
  * Kambe et al. (2013) Kambe E., et al., 2013, PASJ, 65, 15
  * Kervella et al. (2004) Kervella P., Thévenin F., Di Folco E., Ségransan D., 2004, A&A, 426, 297
  * Kobulnicky & Fryer (2007) Kobulnicky H. A., Fryer C. L., 2007, ApJ, 670, 747
  * Konacki et al. (2010) Konacki M., Muterspaugh M. W., Kulkarni S. R., Hełminiak K. G., 2010, ApJ, 719, 1293
  * Korth et al. (2021) Korth J., Moharana A., Pešta M., Czavalinga D. R., Conroy K. E., 2021, Contributions of the Astronomical Observatory Skalnate Pleso, 51, 58
  * Kozai (1962) Kozai Y., 1962, AJ, 67, 591
  * Lee et al. (2019) Lee A. T., Offner S. S. R., Kratter K. M., Smullen R. A., Li P. S., 2019, ApJ, 887, 232
  * Lidov (1962) Lidov M. L., 1962, Planet. Space Sci., 9, 719
  * Maoz et al. (2014) Maoz D., Mannucci F., Nelemans G., 2014, ARA&A, 52, 107
  * Marcadon et al. (2020) Marcadon F., Hełminiak K. G., Marques J. P., Pawłaszek R., Sybilski P., Kozłowski S. K., Ratajczak M., Konacki M., 2020, MNRAS, 499, 3019
  * Mardling & Aarseth (2001) Mardling R. A., Aarseth S. J., 2001, MNRAS, 321, 398
  * Mason et al. (2009) Mason B. D., Hartkopf W. I., Gies D. R., Henry T. J., Helsel J. W., 2009, AJ, 137, 3358
  * Maxted et al. (2020) Maxted P. F. L., et al., 2020, MNRAS, 498, 332
  * Mayer et al. (2013) Mayer P., Harmanec P., Pavlovski K., 2013, A&A, 550, A2
  * Moe & Kratter (2018) Moe M., Kratter K. M., 2018, ApJ, 854, 44
  * Naoz & Fabrycky (2014) Naoz S., Fabrycky D. C., 2014, ApJ, 793, 137
  * Nelder & Mead (1965) Nelder J. A., Mead R., 1965, The Computer Journal, 7, 308
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton B., et al., 2018, ApJS, 234, 34
  * Perets & Fabrycky (2009) Perets H. B., Fabrycky D. C., 2009, ApJ, 697, 1048
  * Petit et al. (2004) Petit P., et al., 2004, MNRAS, 348, 1175
  * Prša et al. (2016) Prša A., et al., 2016, ApJS, 227, 29
  * Raghavan et al. (2010) Raghavan D., et al., 2010, ApJS, 190, 1
  * Randich et al. (2013) Randich S., Gilmore G., Gaia-ESO Consortium 2013, The Messenger, 154, 47
  * Rappaport et al. (2013) Rappaport S., Deck K., Levine A., Borkovits T., Carter J., El Mellah I., Sanchis-Ojeda R., Kalomeni B., 2013, ApJ, 768, 33
  * Rappaport et al. (2022) Rappaport S. A., et al., 2022, MNRAS, 513, 4341
  * Rein & Liu (2012) Rein H., Liu S. F., 2012, A&A, 537, A128
  * Rein & Tamayo (2015) Rein H., Tamayo D., 2015, MNRAS, 452, 376
  * Ricker et al. (2015) Ricker G. R., et al., 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003
  * Rucinski (1999a) Rucinski S., 1999a, Turkish Journal of Physics, 23, 271
  * Rucinski (1999b) Rucinski S., 1999b, in Hearnshaw J. B., Scarfe C. D., eds, Astronomical Society of the Pacific Conference Series Vol. 185, IAU Colloq. 170: Precise Stellar Radial Velocities. p. 82 (arXiv:astro-ph/9807327)
  * Shatsky & Tokovinin (2002) Shatsky N., Tokovinin A., 2002, A&A, 382, 92
  * Simon & Sturm (1994) Simon K. P., Sturm E., 1994, A&A, 281, 286
  * Tamayo et al. (2020) Tamayo D., Rein H., Shi P., Hernandez D. M., 2020, MNRAS, 491, 2885
  * Tokovinin (2004) Tokovinin A., 2004, in Allen C., Scarfe C., eds, Revista Mexicana de Astronomia y Astrofisica Conference Series Vol. 21, Revista Mexicana de Astronomia y Astrofisica Conference Series. pp 7–14
  * Tokovinin (2017) Tokovinin A., 2017, ApJ, 844, 103
  * Tokovinin (2021) Tokovinin A., 2021, Universe, 7, 352
  * Tokovinin (2022) Tokovinin A., 2022, ApJ, 926, 1
  * Toonen et al. (2020) Toonen S., Portegies Zwart S., Hamers A. S., Bandopadhyay D., 2020, A&A, 640, A16
  * Torres et al. (2010) Torres G., Andersen J., Giménez A., 2010, A&ARv, 18, 67
  * Voges et al. (1999) Voges W., et al., 1999, A&A, 349, 389
  * Zasche & Paschke (2012) Zasche P., Paschke A., 2012, A&A, 542, L23
  * Zasche et al. (2023) Zasche P., Vokrouhlický D., Barlow B. N., Mašek M., 2023, The Astronomical Journal, 165, 81
  * Zucker & Mazeh (1994) Zucker S., Mazeh T., 1994, ApJ, 420, 806
  * von Zeipel (1910) von Zeipel H., 1910, Astronomische Nachrichten, 183, 345

##  Appendix A Distance and reddening estimates from the isochrones

The reddening-free distance $d_{0}$ are estimated simultaneously with the
reddening $E(B-V)$, using the available observed total magnitudes in different
filters, and the predicted total brightness of the system in the same filters,
for a given triplet of points (= stellar masses) on the same isochrone. We use
the $T_{\rm{}_{eff}}$–surface brightness relations from Kervella et al.
(2004), to calculate the distances $d_{\lambda}$ and distance moduli
$(m-M)_{\lambda}=5\log(d_{\lambda})+5$ in each band. To obtain the extinction-
free modulus $(m-M)_{0}$ we fit a straight line on the $A_{\lambda}$ vs.
$(m-M)_{\lambda}$ plane, where the $A_{\lambda}$ are extinction coefficients
in each band. We followed the extinction law of Cardelli et al. (1989)
$A_{U}:A_{B}:A_{V}:A_{R}:A_{I}:A_{J}:A_{H}:A_{K}=4.855:4.064:3.1:2.545:1.801:0.88:0.558:0.36$,
which assumes $R_{V}=3.1$. The slope of the fitted line in this approach is
the reddening $E(B-V)$, while the intercept is the extinction-free modulus
$(m-M)_{0}$, which can be translated into the distance $d_{0}$. In the
isochrone fitting process, where the distance is used as one of the
constraints, the reproduced $d_{0}$ value is the one that is being compared
with $d$.

##  Appendix B Radial Velocities

The RVs for BD44, extracted using TODCOR method, are given in Table.7.

Table 7: RV measurements of BD44. BJD-2450000 | $v_{1}\,(km\,s^{-1})$ | $\epsilon_{1}\,(km\,s^{-1})$ | $v_{2}\,(km\,s^{-1})$ | $\epsilon_{2}\,(km\,s^{-1})$ | $\gamma\,(km\,s^{-1})$ | $\epsilon_{\gamma}\,(km\,s^{-1})$ | $v_{3}\,(km\,s^{-1})$ | $\epsilon_{3}\,(km\,s^{-1})$
---|---|---|---|---|---|---|---|---
7022.350624 | -21.499 | 0.697 | -146.291 | 0.505 | -81.664 | 0.754 | -78.406 | 0.079
7061.170892 | 13.564 | 0.428 | -149.393 | 0.973 | -65.001 | 0.956 | -106.205 | 0.100
7062.263807 | — | — | -3.669 | 0.863 | — | — | — | —
7110.041913 | -22.493 | 0.471 | -139.493 | 0.394 | -78.901 | 0.654 | -78.143 | 0.075
7111.161334 | -160.274 | 0.578 | 3.675 | 0.355 | -81.230 | 0.878 | -76.942 | 0.086
7142.999934 | -148.156 | 0.714 | -23.570 | 0.426 | -88.090 | 0.746 | -63.215 | 0.083
7144.008005 | -13.377 | 0.565 | -161.962 | 1.019 | -85.013 | 0.929 | -62.731 | 0.135
7148.097089 | -16.490 | 1.381 | -160.184 | 0.379 | -85.768 | 1.023 | -62.282 | 0.077
7490.046763 | -164.541 | 1.803 | -2.964 | 0.496 | -86.641 | 1.250 | -68.064 | 0.227
7526.015228 | -10.156 | 2.387 | -159.806 | 0.391 | -82.306 | 1.452 | -76.503 | 0.095
7528.026002 | -164.479 | 1.040 | 7.485 | 0.497 | -81.571 | 1.032 | -77.027 | 0.140
7530.053976 | -13.474 | 1.193 | -154.691 | 0.526 | -81.558 | 0.965 | -78.010 | 0.087
7539.052136 | — | — | -41.870 | 0.339 | — | — | — | —
7540.110845 | 3.205 | 0.606 | -166.138 | 0.514 | -78.439 | 0.925 | -82.127 | 0.123
7755.261393 | -6.755 | 0.688 | -164.869 | 0.429 | -82.985 | 0.881 | -69.878 | 0.103
7813.057543 | -146.164 | 2.782 | 2.688 | 0.778 | -74.399 | 1.660 | -92.769 | 0.105
7813.175027 | -137.988 | 1.525 | -8.703 | 0.446 | -75.657 | 1.037 | — | —
7814.125822 | -3.179 | 1.464 | -141.592 | 0.745 | -69.911 | 1.081 | -93.820 | 0.093
7816.093416 | -150.858 | 2.080 | 12.322 | 0.828 | -72.185 | 1.402 | -95.491 | 0.128
7816.322403 | -152.826 | 2.013 | 13.918 | 0.420 | -72.435 | 1.343 | -95.617 | 0.126
7846.107080 | 16.267 | 1.456 | -123.088 | 0.526 | -50.919 | 1.051 | -139.779 | 0.231
7891.177386 | -14.944 | 1.876 | -163.681 | 0.491 | -86.653 | 1.240 | -66.297 | 0.118
7892.950768 | -152.526 | 0.821 | -12.576 | 0.450 | -85.053 | 0.839 | -66.027 | 0.078
7894.027357 | — | — | -154.775 | 0.303 | — | — | — | —
7950.023399 | -3.531 | 0.670 | -173.935 | 0.750 | -85.687 | 0.978 | -62.241 | 0.097
7954.983416 | -164.296 | 0.700 | -3.259 | 0.266 | -86.656 | 0.882 | -62.650 | 0.101
8066.362164 | -156.728 | 2.378 | 14.806 | 0.152 | -74.028 | 1.496 | -91.838 | 0.117

##  Appendix C Spot Parameters

The phoebe 2 model for BD44 included one spot on the Aa star and two on the Ab
star. While we considered the most stable spot (spot-2 on Ab) for our
calculations of differential rotation, the other parameters are important
while modelling the LC and are given in Table.8. The spot parameters vary over
the three TESS sectors: Sector-16 (S16), Sector-22 (S22), and Sector-49 (S49).

Table 8: Spot parameters for spots on primary (Aa) and secondary (Ab) of BD44. Parameters | S16 | S22 | S49
---|---|---|---
$T^{\mathrm{spot}}_{\mathrm{Aa}}$ | 0.9535 | 0.9362 | 1.04848
$T^{\mathrm{spot1}}_{\mathrm{Ab}}$ | 0.9922 | 0.9806 | 0.9545
$T^{\mathrm{spot2}}_{\mathrm{Ab}}$ | 0.8988 | 0.9296 | 0.8499
$r^{\mathrm{spot}}_{\mathrm{Aa}}$ [deg] | 29.449 | 29.846 | 33.943
$r^{\mathrm{spot1}}_{\mathrm{Ab}}$ [deg] | 31.9971 | 30.446 | 11.390
$r^{\mathrm{spot2}}_{\mathrm{Ab}}$ [deg] | 42.846 | 37.311 | 39.611
$c^{\mathrm{spot}}_{\mathrm{Aa}}$ [deg] | 80.735 | 69.893 | 80.735
$c^{\mathrm{spot1}}_{\mathrm{Ab}}$ [deg] | 90.977 | 89.047 | 90.977
$c^{\mathrm{spot2}}_{\mathrm{Ab}}$ [deg] | 22.000 | 35.745 | 21.913
$l^{\mathrm{spot}}_{\mathrm{Aa}}$ [deg] | 180.000 | 194.891 | 84.0973
$l^{\mathrm{spot1}}_{\mathrm{Ab}}$ [deg] | 162.488 | 196.044 | 262.305
$l^{\mathrm{spot2}}_{\mathrm{Ab}}$ [deg] | 0.000 | 29.975 | 180.233

##  Appendix D BF Fitting Tables

The BF fitting was done on multiple spectra of different epochs. Though the
profile of every fit was similar, the flux fraction of each component varied
slightly, which was necessary to consider while spectral disentangling. This
variation was noticed for $v\mathrm{sin}(i)$ too but their variations were not
significant to consider during spectral analysis. The complete tables for
these parameters are available in the online version.