Investigating Period Variability Mechanisms in Eclipsing Binary Stars through Eclipsing Time Variation Analysis: A Case Study of TZ Bootis

We present an effective strategy for extensive analysis of eclipsing time variations (ETVs) using modern and sophisticated optimization methods that comprise a set of tools to investigate period variability mechanisms in eclipsing binary stars such as the light-time effect, the Applegate mechanism, and mass transfer. We implement these methods for the first time assuming that the above mechanisms can act simultaneously in the puzzling W UMa–type binary star TZ Bootis by using archival and new TESS data spanning 75 yr and reexamining the up-to-date ETVs. Preliminary analysis of the TESS data revealed for the first time the presence of a second binary in agreement with previous spectroscopic data and astrometric results from Gaia DR3. We consider the most credible scenario for the ETV: two stellar circumbinary companions of minimum masses M 3 = 0.5 M ☉ and M 4 = 0.14 M ☉ in highly eccentric orbits e 3 = 0.70 and e 4 = 0.82 with periods P 3 = 38 yr and P 4 = 20 yr along with a 24 yr magnetic activity of the secondary component and a long-term period increase (dP/dt = 1.2 × 10−8 days yr−1), interpreted as a conservative mass transfer from the secondary to the primary component at a rate of dM 1/dt = 3.7 × 10−9 days yr−1. Further spectroscopic observations, analytical modeling of the second pair, and ETV analysis of both pairs are needed to investigate the quadruple nature of the system.


Introduction
O − C (observed minus calculated) eclipse timing residual diagrams can be a useful tool for investigating period variability mechanisms in eclipsing binary (EB) stars since observed eclipsing time variation (ETV) may reveal periodic, irregular, or secular behavior.Each of these trend behaviors can be attributed to a physical or apparent mechanism such as mass and/or angular momentum transfer or loss, magnetic activity, apsidal motion, third or more bodies in the system, and other intrinsic evolutionary processes of the binary star.
Long-term cyclic period changes (years to decades) are common in many classes of close binary systems (Mayor & Mazeh 1987;Chambliss 1992;Tokovinin 1997).The lighttime effect (LITE) due to a third body around the binary (Irwin 1952;Pribulla & Rucinski 2006), the Applegate mechanism caused by the magnetic activity of one or both components (Applegate 1992), and apsidal motion in eccentric systems have been proposed as the most plausible explanations for these variations.The latter unlike the other mentioned mechanisms affects the primary and secondary eclipses in different ways (Lacy 1992;Wolf et al. 2013).There may also be additional mechanisms, such as modulation through spots (Kalimeris et al. 2002;Tran et al. 2013) and pulsation (Hełminiak et al. 2017;Gaulme & Guzik 2019), that result in erroneous ETVs.
In the new era of space-based photometry, the importance of multiplicity, especially the interest in exoplanets in multiplestar systems, has triggered many efforts in the detection of LITE (e.g., Conroy et al. 2014;Borkovits et al. 2016;Hajdu et al. 2022b;Rappaport et al. 2022).In tight triples, LITE is supplemented by the effects of dynamical perturbations and other phenomena (see Borkovits et al. 2016).Recently Borkovits et al. (2020Borkovits et al. ( , 2021) ) presented spectrophotodynamical analyses of triply eclipsing triples discovered with the Transiting Exoplanet Survey Satellite (TESS), while Tokovinin (2022) compiled a sample of 392 low-mass hierarchical triple stellar systems within 100 pc based on Gaia data.All of these discoveries highlight the necessity of combining different data sets in order to confirm the existence of three or more companions in a binary-star system.Such confirmation is of great usefulness mainly to filling the gaps in our understanding of the origin and evolution of stellar and planetary systems.
On the other hand, the detection of ETVs through O − C diagram analysis has many caveats due mainly to problems arising in the determination of the times of minima (TOMs; for a detailed review see Nelson et al. 2014).Equally important is, however, the nature of the selected optimization techniques to solve the inverse mathematical problem by finding the bestfitting curve to the O − C diagram, due to nonlinear equations and a strong correlation of parameters responsible for a parameter space topology full of valleys and peaks (Kallrath & Milone 2009).This is why several studies of ETV analysis were not supported by the most recent data.Papageorgiou (2015) and Papageorgiou & Christopoulou (2015) suggested a fitting strategy for O − C analysis that incorporates optimization methods to explore the parameter space such as genetic algorithms (GAs;Holland 1975;Charbonneau 1995), heuristic scanning with parameter perturbation, the combination of Nelder-Mead downhill simplex and Levenberg-Marquardt (NMS-LM) algorithms, and Markov Chain Monte Carlo (MCMC)/Monte Carlo-based methods for error budget and parameter determination (Papageorgiou et al. 2023).Parimucha et al. (2018) and Gajdoš & Parimucha (2019) also proposed a fitting strategy that uses GA optimization followed by MCMC fitting for the O − C analysis of extrasolar planets and Algoltype EBs.We follow here the same fitting strategy enriched by new methods for global optimization.
In this paper, we present an effective strategy for ETV analysis using modern and sophisticated optimization methods in order to overcome many problems in the best available manner.It is also the first ETV analysis where both mechanisms, LITE and Applegate, are assumed to work simultaneously.In the case of the latter mechanism we implement for the first time the analytical approach of original models (Völschow et al. 2016) using the radial density profile of a polytropic star with index n = 1.5.We therefore choose an EB with a strong period change that may indicate its triple or quadruple nature as confirmed by the published spectroscopic data.This is the TZ Boo system, for which we present an extensive analysis of the ETV updated by all previously published TOMs and including new minima from the TESS space telescope.
We formulate the fundamental mathematical framework in Section 2 for the LITE-and Applegate-modified ETV analysis.In Section 3 we outline the steps of the proposed method of optimization algorithms implemented here for the case of TZ Boo.The compilation of TOMs and the derivation of the new TOMs of TZ Boo are described in Section 4, together with the results of the ETV analysis for three cases.Finally, we discuss conclusions in Section 5.

LITE
One or more bodies orbiting a binary star result in the periodic variation of the eclipsing pair distance from the observer.Consequently, an apparent period variation of the EB can be observed.This phenomenon of a binary's oscillating motion around the common barycenter of a triple or multiple system is called LITE and was first formulated by Irwin (1952Irwin ( , 1959)).
The light time, τ LITE , is the observable minimum time variation defined as (Irwin 1952;Borkovits et al. 2015) where A LITE is the observed semiamplitude of LITE variation in days while c = 2.590 × 10 10 km day −1 is the velocity of light.The parameters υ 3 , ω 3 , e 3 , and i 3 (the true anomaly, argument of periastron, eccentricity, and inclination, respectively) describe the orbit of the third body and i sin 12 3 a is the projected semimajor axis of the eclipsing pair around the common barycenter with the third body.The true anomaly depends on the period P 3 and the time of periastron passage Tp 3 and hence they do not explicitly appear in Equation (1).
We can define the time of minimum light (Min I ) for an EB of period P bin with a suspected companion as where T 0 is a reference TOM, while cases of continuous period change with a constant dP/dt such as mass transfer, mass outflow, or magnetic braking are described by the quadratic coefficient Q: Consequently, we can define the merit function for our optimization problem as the sum of normalized square residuals over N data set values of observed minimum times taking into account the associated uncertainties σ i or the equivalent weights denoted by w i .The observed values are represented by O, while the computed TOM (Equation ( 3)) is denoted by C.
The set of eight parameters (T 0 , P bin , e 3 , A, ω 3 , P 3 , Tp 3 , Q) that minimizes the χ 2 (Equation ( 5)) describes the period variation of an EB depending on the quantity and quality of the data set.A good constraint of parameters can be achieved by a sufficient number of data points of high accuracy and density; however, this is not always the case.

The Applegate Mechanism
According to the Applegate model (Applegate 1992) a magnetically active star in a binary can transfer angular momentum toward the surface from its interior and backward during a complete magnetic cycle, leading to oscillations of the orbital period.The magnetic star's energy budget must be sufficient to adjust its quadrupole moment for this mechanism to work, resulting also in a variable luminosity (Equations ( 28) and (30) in Applegate 1992).
One of the model's observable predictions is the period change relative to the binary period given by Therefore, we can construct the relevant O − C equations and χ 2 merit function like the ones defined for the LITE case (Equation ( 5)) by replacing the LITE term with an Applegate one.
In addition, as the Applegate model predicts that the magnetically active star gets bluer as it brightens, orbital period changes should be consistent with the light and color variations of the system during the same cycle or with any other evidence of magnetic activity (e.g., starspots, emission cores in Ca II or Mg II lines, the X-ray emission, etc.).In most cases such luminosity and color variations are difficult to detect with certainty either because there are insufficient data in the same filter(s) or because the variations are smoothed in time due to the fact that the Kelvin-Helmholtz relaxation timescale in the convective envelope is much longer than the timescale for the resulting luminosity variations (Watson & Marsh 2010;Khaliullin & Khaliullina 2012).Furthermore, when asymmetries in the light curves (LCs) are present the problem becomes more complex.Völschow et al. (2016) proposed a more accurate and realistic two-zone model, where they treated the stellar density profile of the inner and outer regions for two different densities including the change of the quadrupole moment both in the shell and in the core, adding to previous improvements of the mechanism (Brinkworth et al. 2006;Lanza 2006;Tian et al. 2009).Their formula for the minimum energy threshold is given by where k 1 and k 2 are two numerically evaluated coefficients.
Here we implement for the first time an analytical approach to the original thin-shell model of Applegate (1992) and the two-zone model of Völschow et al. (2016) using the radial density profile ρ(r) of a polytropic star with index n = 1.5.Based on the assumption of an upper shell mass limit M shell = 0.1 M of the active star (Applegate 1992) we calculated the minimum energy thresholds (Equation ( 8)) for shell mass values in the range (0, 0.1 M).The distribution of mass values within this range was calculated by solving the Lane-Emden equation for polytropic index n = 1.5 in order to acquire a radial density profile, approximating the secondary star as a fully convective star.
In the case of the thin-shell model, we considered a core of fixed mass concentrated on a point at the center of the shell (Applegate 1992), whereas, in the two-zone model, we defined a core and a shell region of two different densities (Völschow et al. 2016) according to the shell mass range (0, 0.1 M) and the radial density profile of the polytrope.Therefore, we were able to calculate numerically the k 1 and k 2 coefficients of Equation (8) following the same approach of Völschow et al. (2016) but using the approximation of a polytropic radial profile instead of a stellar profile from an evolutionary code.

Optimization Algorithms
Briefly, the strategy that we used was comprised of the following steps: 1. First a solution was approximated using global optimization algorithms such as GAs, differential evolution There is extensive literature about the basic concepts and details of all the algorithms that we use in this research and therefore we present them only by their names and references.Their selection was based on their wide use in the field of natural sciences and their public availability through Python libraries (for details see the software section).

Period Variation Analysis of TZ Boo
TZ Boo is one of the most interesting and unusual W UMa binaries because it has exhibited several period changes and variations in both minima's depth since its discovery by Guthnick and Prager (1926).Qian & Liu (2000) pointed out several jumps in the orbital period.Awadalla et al. (2006) showed that the ETV residuals' variations may be caused by a magnetic activity cyclic mechanism.Christopoulou et al. (2011) found a cyclic variation of about 31.2 yr and suggested that the ETV (1948ETV ( -2010) ) may indicate that TZ Boo is either a triple system with M 3 = 0.99 M ☉ (coplanar case) and e 3 = 0.63 ± 0.06 or a quadruple system as confirmed from the published spectroscopic data (Pribulla et al. 2009).Nelson et al. (2016) in light of new data  proposed a twocompanion LITE model with e 3 = 0.83 ± 0.01, e 4 = 0.303 ± 0.032, P 3 = 45.5 ± 0.6 yr, P 4 = 22.4 ± 0.1 yr, and a secular increase of the period at a rate of dP/dt = (4.6 ± 0.1) × 10 −8 days yr −1 .However, a previous adaptive optics search (Rucinski et al. 2007) did not reveal a companion at a separation <5″ although a very close companion was suspected within 0 3-0 4. Gaia DR3 does not distinguish stars that are ∼ 0 5 apart and lists TZ Boo as a single object astrometrically and spectroscopically with a parallax of π = 5.3956 ± 0.0740 mas and g log 3.9921 0.0835 0.1064

= -+
. However, the so-called astrometric excess noise parameter (AEN) is as large as ∼11,776σ.AEN is a potential signature of the orbital wobble of individual components in a binary system or of the inner system in a triple.This is also supported by the large value of the renormalized unit weight error (RUWE) 6.135 as a RUWE value greater than 1.4 can typically be an indicator of systems with unresolved companions (Stassun & Torres 2021).
The LC modeling of a system with a dark spot on the more massive component (Christopoulou et al. 2011) resulted in absolute parameters of the components with values of

Data Acquisition and TOM Determination
We compiled all available TOMs of TZ Boo from Christopoulou et al. (2011) and Nelson et al. (2016) updated up to 2022 with 54 new TOMs from the literature.Out of a total of 530 data points, there are 36 visual, 158 photoelectric, and 336 CCD data covering the period from 1948 to 2022 (Table A1).
In addition, TZ Boo (TIC 20212631) was observed by TESS in 2020 (sectors 23 and 24) and 2022 (sectors 50 and 51) at a 2 minute cadence.The LCs are available from the Mikulski Archive for Space Telescopes (MAST). 1o determine a handful of TOMs instead of hundreds from the entire range of TESS time series, we phase-folded and binned each time-series LC into 1000 equally spaced phase bins of normalized flux using the Lightkurve Python package for Kepler and TESS data analysis (Lightkurve Collaboration et al. 2018).The steps of the computational process are as follows: 1.A seventh-order polynomial is fitted separately on the primary and secondary minima of the phase-folded and binned LC (FBLC) following Hajdu et al. (2022a).2. Its minimum is found by minimizing the best-fit polynomial using the NMS algorithm followed by an LM least-squares minimization to find the midphase of the primary and secondary eclipses.3. The resulting phase and flux values for each eclipse are used as query points for a KDTree nearest-neighbor searching algorithm (Maneewongvatana & Mount 1999) resulting in a unique pair of primary and secondary eclipses of the phase-folded but unbinned time-series LC. 4. Finally, we apply the same process of polynomial fitting and determine the midtime of the unique pair of eclipses.Reliability of errors is ensured using the affine-invariant MCMC sampler emcee (Foreman-Mackey et al. 2013) as shown by Papageorgiou et al. (2021).
In Figure 1(a) we present a segment of the time series of TESS sector 50 determining a primary and a secondary TOM.In total 13 segments of continuous observations were extracted from the four available TESS LC files and 26 TOMs were calculated from their FBLCs.These were converted from BJD to HJD by using the publicly available code of Eastman et al. (2010) and are listed in Table 1.For comparison, we also calculated the TOMs using the method of Kwee & van Woerden (1956;KvW) and found that 80% of the TOMs are in agreement within the limits of our MCMC error.As the error of  Note.TESS TOMs with their error estimation calculated using the emcee MCMC sampler and the KvW method.
KvW is usually underestimated as it reflects the assumption of a perfectly symmetrical LC, we used the calculated TOM values and errors of the emcee MCMC sampler.
We have to point out that Figure 1(b) clearly shows a marked LC asymmetry and our preliminary analysis of maximum flux level variations for all the sectors of TESS data clearly shows the O'Connell effect, which can be attributed to the presence of spots in one or both components of the system.To quantify the expected minimum shift Δt due to spots, we used the approximation of Pribulla et al. (2012, Equation (A5)) using the amplitude of the O'Connell effect (∼0.006-0.03).The resulting shift is 2.5-17 s, smaller than our error in Table 1.
In addition, the TESS LCs show brightness dips with a strict periodicity of 4.75 days, and thus we searched for light contamination due to the large pixel scale (21″ × 21″) of the TESS CCD from the star TIC 2021633 (T = 16.78 mag) accounting for its proximity (11 77) but its contribution was found to be negligible.Nevertheless, our preliminary analysis resulted in a disentanglement of the target signal (Figure 1(c)) with Silicups2 and revealed a secondary one with a period of 9.5 days (Figure 1(d)), which can be attributed to the presence of a detached single-line binary of a 9.5 day period.This result agrees with half of this period's brightness drops in the TESS LCs and with the spectroscopic results of Pribulla et al. (2009) of a single-line binary with a period of 9.48 days.Our tests of calculating the minima from the disentangled TESS LC of TZ Boo showed that the mean difference from the calculated minima of Table 1 is 0.0005 days and thus they are within the error limits.
By collecting all 556 TOMs in the period 1948-2022 and using a linear ephemeris (Christopoulou et al. 2011)

Min
HJD 2452500.1602170.2971599316.9 The epochs and the O − C values of all TOMs were calculated (Table A1) and are illustrated in Figure 2. Figure 2 shows the ETV nature of the system throughout the past 74 yr .To overcome the problem of old archival TOMs with low precision or without errors we adopted the weighting scheme of Kallrath & Milone (2009) and Wolf et al. (2016Wolf et al. ( , 2018) ) and assigned w v = 1 for visual observations and w ccd,pe = 10 for CCD + photoelectric observations to all TOMs of Table A1 accordingly (related to variances σ v = 0.0019 and σ ccd,pe = 0.0006, respectively).These are not affected by the additional signal, as they fall within the error limits.
We ensured through many optimization runs that the first visual and photoelectric data points that seem isolated at the left side of Figure 2 (1948)(1949)(1950)(1951)(1952)(1953)(1954)(1955)(1956)(1957)(1958)(1959)(1960)(1961) do not drive the optimization process to a different solution whether they are included or discarded; in the end, we chose to include them in our analysis.
Fitting a quadratic equation (a mass transfer term) together with a single sinusoidal fit (one circumbinary companion; Equation (3)) or with Applegate's mechanism resulted in a poor fit (not presented here) with residuals yielding the presence of an additive phenomenon.Therefore we applied a global and local vicinity search of the best solution as we described in Section 3 for three cases: 1. two circumbinary companions and mass transfer (model 1) 2. one circumbinary companion, magnetic activity, and mass transfer (model 2) 3. two circumbinary companions, magnetic activity, and mass transfer (model 3)

Model Fitting
Adding a second LITE term to Equation (3) leads to and increases the set of adjusted parameters to 13 in order to minimize the merit function (Equation ( 5)) of model 1.
Equivalently the equations for models 2 and 3 are defined as with a total of 10 and 15 parameters for adjustment, respectively.All minimization algorithms (global and local) were implemented as described in Section 3, and during the minimization process, all parameters were freely adjusted but bound to a wide range of physically accepted values.The following process was used for each model.
As a first step, we applied the global optimization algorithms (GA, DE, and GSA) separately to approach the global minimum as a first estimation of the best-fitting curve as described in Section 3. Afterward, this estimation formed the starting point of a local vicinity search implemented by the NMS-LM algorithms.
The stability of the solution was tested by a heuristic scanning of 5000 NMS-LM iterations with parameter kicking using a standard normal distribution as a random offset.This random effect essentially acted as a perturbation to each parameter before the start of each iteration with value k * N(0, 1) * p, where k is the size of the perturbation defined by the user, N(0, 1) is the standard normal distribution, and p is the value of each parameter.The purpose of this process was the exhaustive scan of the parameter space inside the local minimum while at the same time a large perturbation value in the range 5%-10% should be enough to escape the local minima and drive the solution to another area of parameter space (Prša & Zwitter 2005).
Finally, the resulting solution of heuristic scanning was used as a starting point for a Bayesian regression using the NUTS MCMC sampler.Normal prior distributions were selected for our model parameters while the sampling process was set with a 3000-sample burning period, 5000 samples as draw iterations, and 15 chains resulting in a total of 75,000 samples.
Trace plots offer a visual indicator to assess the convergence of MCMC runs and as shown in Figure 3(a) there is a good mixing of chains for all models suggesting that the sampler explored the parameter space efficiently.Another indicator for assessing the convergence and diagnosing poor exploration is to plot the energy transition distribution and marginal energy distribution to quantify how well momentum resampling matches the marginal energy distribution according to the BFMI (Betancourt 2016).
In Figure 3(b) the momentum resampling effectively generates exact draws from the marginal energy distribution indicated by the well-matched distributions and BFMI values close to unity for all Markov chains.This indicates a fast random walk that explores efficiently the target distribution with small autocorrelations between samples.
The WAIC and LOO indicators provide sophisticated measures of model comparison to assess the predictive performance of models and account for both model fit and complexity using penalization terms to avoid overfitting.Based on this cross-validated estimate of model performance we can The results of all optimization methods are presented in Table 2, and Figures 4(a)-(c) contain the best-fitting curve for each model together with the residuals of the synthetic curve from data points at the lower panels.
The posterior and joint distributions of the parameters are shown in Figure A1 and in the online version of the appendix.The rank-normalized split R-hat indicator (Vehtari et al. 2021) was estimated during the optimization process with value R 1 ˆ= suggesting a good convergence and mixing of MCMC chains.
Adopting the results of MCMC for the quadratic coefficient (Equation ( 4)) we estimated a long-term secular period increase (dP/dt) with its values presented in Table 2 for each model.Mass transfer is expected due to the binary's contact nature, and the resulting positive value of the quadratic term can be interpreted as conservative mass transfer from the secondary to the primary member.Using the equation of Hilditch (2001) we calculated a mass transfer rate dM 1 /dt = (2.6 ± 0.3) × 10 −9 M ☉ yr −1 for the case of two LITEs (model 1).Notes.Optimization results for the three composite models.The number in parentheses is the numerical value of the uncertainty referring to the corresponding last digits of the reported result.The global optimization algorithms (GA, DE, and GSA) do not provide an error value due to their fundamental computational principles.
Energy requirements for the feasibility of the Applegate mechanism were investigated for the secondary star using the analytical approach of the original thin-shell model of Applegate (1992) and the two-zone model of Völschow et al. (2016) as described in Section 2.2.
The comparison of the two approximations for model 2 is plotted in Figure 4(d) showing the strong feasibility of the mechanism (ΔE/E 2 < 1) for (0.009−0.1)M 2shell /M 2 = (0.002-0.021)M ☉ in the case of the thin-shell model and for (0.012−0.1)M 2shell /M 2 = (0.003-0.021)M ☉ in the case of the two-zone model.For the same models and M 2shell = 0.1 M 2 the relative luminosity variation (ΔL/L 2 ) was calculated to be 0.28 and 0.40, respectively.
In the case of model 3 we calculated the minimum energy threshold of the Applegate mechanism for the same range of shell mass values and the relative luminosity variation (ΔL/L 2 ) was found to be 0.30 and 0.41 in the limit M 2shell = 0.1 M 2 for the thin-shell model and two-zone model, respectively.The almost identical values to those of model 2 result from their very similar period modulation and amplitude values.

Conclusions
We present a strategy for extensive analysis of ETVs using modern and sophisticated optimization algorithms to investigate thoroughly the period variability in EBs.Our strategy is a first-solution approximation using global optimization algorithms such as GA, DE, and GSA followed by local optimization ones (NMS, LM, HS, and NUTS MCMC sampler).However one could interchangeably use these two global and local optimization routes.Applying the complete set of algorithms on the 75 yr ETV diagram of the W UMa-type EB TZ Boo, we investigated three models-(i) two LITEs, (ii) LITE and Applegate, and (iii) two LITEs and Applegate-considering (iii) as the most plausible for the following reasons.
Our preliminary analysis of TESS LCs resulted for the first time in the disentanglement of the target and revealed a detached binary of a ∼9.5 day period in accordance with the spectroscopic study of Pribulla et al. (2009); however, their spectroscopic ephemeris does not predict precisely the minima of the second pair.The same authors suggested that TZ Boo is more likely a hierarchic quadruple system consisting of a contact EB and a second, non-eclipsing single-line binary with a period of 9.48 days (although the possibility of eclipses was not ruled out), revolving in a 34 yr orbit.From our orbital parameters and the distance to the system from Gaia DR3 (d = 210.9629pc), the computed angular separations of the third and fourth components are about 60 and 42 mas, respectively.The late spectral type of TZ Boo, the asymmetry of the TESS LCs, and the spotted (although not unique) photometric solution of Christopoulou et al. (2011) support the theoretical considerations about a magnetic dynamo mechanism and a magnetic activity that exists simultaneously with two LITEs, and lead to the observed ETV diagram.This model is also supported by the most recent TOMs (see Table A1) when they are added to our ETV diagram (not presented here).
We investigated the feasibility of the Applegate mechanism for the first time simultaneously with LITE, by approximating the secondary star of TZ Boo as a polytrope (n = 1.5) and calculating its radial density profile.Our numerical solutions revealed adequate energy requirements for M 2shell ∼ (0.003-0.021)M ☉ in the case of the thin-shell and the most accurate two-zone model.
According to the Applegate (1992) theory, if either of the stars in TZ Boo is magnetically active, the mean bolometric luminosity changes globally, in response to the magnetic cycle, under the assumption that no local photospheric inhomogeneities exist.The luminosity changes published by Awadalla et al. (2006) and Elkhateeb & Nouh (2013) using the available LCs through the years 1927-1990 and 1990-2011, respectively, are in agreement with our calculations of L L 0.28 2 min ( ) D ~per activity cycle (23.7 yr).Their reported variation of the amplitude (depth) of the secondary minimum is in the range 0.25-0.28mag over a period of 24 yr  in the V and B filters, respectively, corresponding to a change ΔL/L 2 ∼ ±0.26-0.29.
The corresponding change of the mean magnitude difference relative to the mean light level comes out as Δm ∼ ±0.06-0.09mag.In order to check for maximum light variations, we performed a search on archival photometric data for long-term TZ Boo time series.Unfortunately, the data of TZ Boo in most massive astronomical surveys were saturated.However, we found ∼12 yr time span photometric data series of TZ Boo in the V band and I c band in the Kamogata/Kiso/Kyoto Wide-field Survey (KWS).3Using the method described in Papageorgiou et al. (2018), we attempted to fit linear, parabolic, and sinusoidal models on maximum light data points.Although there is a spurious variation it is not statistically significant in the given time span and photometric accuracy (ΔBIC 1, Verr ∼ 0.07 mag).No significant color variations were observed at maximum light (V -I c ), given the time span and the photometric accuracy.So, we could not exclude the possibility that maximum light variations are present and thus support/reject Applegateʼs model for O − C variations.The obvious O'Connell effect in the historical and TESS LCs and the presence of a second binary make the ETV changes more complex.
TZ Boo is unusual in its depths of minimum change, and the LC shows very large variations of its shape.The presence of circumbinary components or of the revealed second binary may be the cause of the large perturbations observed in the LC of the eclipsing pair.Our next step is to model the high-quality TESS LCs of the eclipsing pair of TZ Boo by studying the O'Connell effect to improve the modeling of the system and shed more light on the surface activity of its components.To investigate the quadruple nature, we would have to model the second pair and study the ETV analysis signals of both pairs on the basis of archival photometry and new data, to confirm that both of them actually orbit a common barycenter (Zasche et al. 2022(Zasche et al. , 2023)).In addition, supplementary methods based on astrometric, spectroscopic, or even interferometric observations and stability analysis should be used for cross-validation.
(DE;Storn & Price 1997), and generalized simulated annealing (GSA;Tsallis & Stariolo 1996;Xiang et al. 1997).2. The resulting parameters of the global optimization algorithms were used as the initial values of the local optimization routine: (a) A first model fit was performed using NMS.(b) The NMS output was narrowed down to more precise parameter values and uncertainties using the LM algorithm.3. The solution stability was tested by using these two methods (NMS-LM) in a heuristic scheme with parameter kicking (HS) for a large number of iterations and an initial grid of solutions (Prša & Zwitter 2005; Christopoulou & Papageorgiou 2013).4. Bayesian regression was performed using the No-U-Turn Sampler (NUTS; Hoffman & Gelman 2014) to explore even more extensively the parameter space and take advantage of various diagnostic convergence tools (trace plots, energy plots, and R-hat (Vehtari et al. 2021); Bayesian fraction of missing information (BFMI; Betancourt 2016); WAIC and LOO indicators (Vehtari et al. 2016)).

Figure 1 .
Figure 1.Example of computing TOMs of TZ Boo from TESS.(a) Segment of TESS time-series data (sector 50) used to construct the normalized FBLC shown in (b).For details see Section 4.1.(c) Disentangled signal of TZ Boo from the TESS LC.(d) Disentangled signal of 9.5 day period.

Figure 3 .
Figure 3. (a) Marginal posteriors and trace plots of sequential samples generated by the NUTS MCMC sampler for e 3 .(b) Energy transition distribution and marginal energy distribution as a diagnostic indicator of the quality of parameter space exploration in the case of a two LITEs and Applegate composite model.The value BFMI → 1 for all Markov chains indicates a fast random walk with small autocorrelations and therefore an efficient exploration and convergence to the target distribution.The complete figure set for all models (41 images) is available in the online journal.(The complete figure set (41 images) is available.)

Figure 4 .
Figure 4. (a)-(c) Upper panels: best-fitting curve to O − C diagram of TZ Boo for the case of long-term period variation (dashed line) corresponding to conservative mass transfer from the secondary to the primary member and (a) two LITE terms (solid line), (b) Applegate and LITE terms (solid line), or (c) Applegate and two LITE terms (solid line).Lower panels: residuals of the best-fitting curve from observational data.(d) Energy threshold comparison for thin-shell model and two-zone model showing strong feasibility of Applegate mechanism for M 2shell = (0.002-0.021)M ☉ and (0.003-0.021)M ☉ , respectively, in the case of Applegate, LITE, and mass transfer (model 3).The complete figure set (five components) is available in the online journal.(The complete figure set (five images) is available.)

Figure A1 .
Figure A1.Posterior probability and joint probability distributions of the first LITE term resulting from Bayesian regression using the NUTS MCMC sampler for model 1.The dashed vertical lines on the one-dimensional posteriors define the 25th, 50th, 75th, and 97th percentiles of the samples.The complete figure set (seven images) is available in the online journal.(The complete figure set (seven images) is available.)

Table 2
Parameters of the Composite Models