Context. Recent period investigations of the postcommon envelop binary PCEB NSVS 14256825 suggest that two circumbinary companions are necessary to explain the observed eclipse timing variations ETVs. Aims. Our objective in this work was to search for the bestfitting curve of two LTTE terms of the ETV diagram by implementing a grid search optimization scheme of Keplerian kinematic and Newtonian Nbody fits alongside a dynamical stability analysis of Nbody simulations. Methods. We compiled two datasets of archival photometric data covering different timelines and updated them with new observations and with three new times of minima calculated from the Transiting Exoplanet Survey Satellite TESS. A grid search optimization process was implemented and the resulting solutions that fell within the 90 confidence interval of the bestfitting curve of the ETV diagram were tested for orbital stability using Nbody simulations and the MEGNO chaos indicator. Results. The Keplerian and Netwonian fits are in close agreement and hundreds of stable configurations were identified for both datasets reaching a lifetime of 1 Myr. Our results suggest that the ETV data can be explained by the presence of a circumbinary planet with mass mSUBbSUB  11 MSUBJupSUB in a nearly circular inner orbit of period PSUBbSUB  7 yr. The outer orbit is unconstrained with a period range PSUBcSUB  20  50 yr from 31 to 71 MMR for a circumbinary body of substellar mass mSUBcSUB  11  70 MSUBJupSUB. The stable solutions of the minimum and maximumreduced chisquare value were integrated for 100 Myr and confirmed a nonchaotic behavior. Their residuals in the ETV data could be explained by a spinorbit coupling model ApplegateLanza. However continuous monitoring of the system is required in order to refine and constrain the proposed solutions.
['binaries: eclipsing', 'planetary systems', 'stars: individual: NSVS 14256825', 'Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - Earth and Planetary Astrophysics']
NSVS 14256825 Period variation and orbital stability analysis of two possible substellar companions
{'NSVS 14256825: Period variation and orbital stability analysis of two possible substellar companions': 'K. Zervas 1 and P.-E. Christopoulou 1  \nDepartment of Physics, University of Patras, 26500 Patra, Greece e-mail: konst.zervas@upatras.gr  \nSeptember 10, 2024', 'ABSTRACT': 'Context. Recent period investigations of the post-common envelop binary (PCEB) NSVS 14256825 suggest that two circumbinary companions are necessary to explain the observed eclipse timing variations (ETVs).  \nAims. Our objective in this work was to search for the best-fitting curve of two LTTE terms of the ETV diagram by implementing a grid search optimization scheme of Keplerian (kinematic) and Newtonian (N-body) fits alongside a dynamical stability analysis of N-body simulations.  \nMethods. Wecompiled two datasets of archival photometric data covering di ff erent timelines and updated them with new observations and with three new times of minima calculated from the Transiting Exoplanet Survey Satellite (TESS). A grid search optimization process was implemented, and the resulting solutions that fell within the 90% confidence interval of the best-fitting curve of the ETV diagram were tested for orbital stability using N-body simulations and the MEGNO chaos indicator.  \nResults. The Keplerian and Netwonian fits are in close agreement, and hundreds of stable configurations were identified for both datasets reaching a lifetime of 1 Myr. Our results suggest that the ETV data can be explained by the presence of a circumbinary planet with mass mb = 11 MJup . in a nearly circular inner orbit of period Pb = 7 yr. The outer orbit is unconstrained with a period range Pc = 20 -50 yr (from 3:1 to 7:1 MMR) for a circumbinary body of substellar mass ( mc = 11 -70 MJup . ). The stable solutions of the minimum- and maximum-reduced chi-square value were integrated for 100 Myr and confirmed a non-chaotic behavior. Their residuals in the ETV data could be explained by a spin-orbit coupling model (Applegate-Lanza). However, continuous monitoring of the system is required in order to refine and constrain the proposed solutions.  \nKey words. binaries: close - binaries: eclipsing - stars: individual: NSVS 14256825 - planetary systems.', '1. Introduction': "A well-known post-common envelop binary (PCEB), NSVS 14256825 (V1828 Aql, TIC 404635917) was discovered as a 13.2 mag variable star in the Northern Sky Variability Survey (NSVS) (Wo'zniak et al. 2004), and it has been recognized as an eclipsing binary by Wils et al. (2007) ( ∆ V ∼ 0.8 mag) with a very short period of P = 0.110374230(2) d and an sdOB + dMtype light curve (LC; OB subdwarf and a M dwarf companion). From high precision BVRcIc LCs, Zhu et al. (2011) derived a photometric mass ratio of 0.22(0.02) and a primary mass around the value for the He core flash at M 1 = 0 . 46 M ⊙ and M 2 = 0 . 1 M ⊙ = 105 MJup for the fully convective secondary star. Almeida et al. (2012) classified NSVS 14256825 as sdOB + dM and obtained the first detailed physical and geometrical parameters from UBVRcIcJH photometric and spectroscopic observations (inclination of the system i = 82 . 5 · ± 0 . 3 · , M 2 / M 1 = 0 . 260 ± 0 . 0012) as M 1 = 0 . 419 ± 0 . 07 M ⊙ , M 2 = 0 . 109 ± 0 . 023 M ⊙ , α = 0 . 80 ± 0 . 04 R ⊙ (separation between the components). Nehir & Bulut (2022), using new BVRI photometry (five nights, 2018-2019) and the old spectroscopy, presented the following updated parameters: i = 82.2 · ± 0 . 1 · , M 1 = 0 . 351 ± 0 . 04 M ⊙ , M 2 = 0 . 095 ± 0 . 01 M ⊙ , and α = 0 . 74 ± 0 . 03 R ⊙ .  \nA large number of PCEBs show apparent period variations, which are frequently attributed to giant circumbinary planets (Zorotovic & Schreiber 2013; Heber 2016) through eclipse timing variation (ETV) caused by the light-travel time e ff ect (LTTE)  \n(Irwin 1952; Borkovits et al. 2015). Marsh (2018) reports at least ten PCEBs with proposed circumbinary brown dwarfs or planets detected with ETV modeling.  \nZorotovic & Schreiber (2013) investigated the origin of these planets with binary population synthesis simulations and concluded that it is very unlikely that they formed in the primordial circumbinary disk (first-generation scenario) and survived the energetic common envelope (CE) evolution. The same authors proposed the formation of the planets from the CE material (second-generation scenario) as being the most likely for PCEBs, although their planetary formation is still an open question, with Bear & Soker (2014) arguing in favor of the firstgeneration hypothesis.  \nPulley et al. (2022) evaluated seven PCEBs with claimed circumbinary planets detected using the ETV method. Remarkably, they revealed that none of the more than 30 circumbinary models that have been proposed through the years have accurately predicted eclipse times within a year. Additionally, most of these models have failed the tests of dynamic stability, making the existence of these circumbinary planets questionable. We give, as an example, HW Vir, which is the prototype of PCEBs, with a subdwarf (sdB) primary star and an M dwarf (dM) secondary. Beuermann et al. (2012b) presented a two-planetary companion model that proved to be stable for at least 10 Myr. However, recent studies (Esmer et al. 2021; Brown-Sevilla et al. 2021; Mai & Mutel 2022) have found dynamically unstable configurations on short timescales.  \nNotably, NSVS 14256825 has a long and conflicting history of reports of possibly having circumbinary companions, with five circumbinary models proposed in the past. Beuermann et al. (2012a) proposed an eccentric ( e = 0 . 5) circumbinary planet of mass 12 MJ and a period of 20 yr, whereas Almeida et al. (2013) suggested a two-planetary model with masses 2.9 MJ and 8.1 MJ ; periods of 3.5 yr and 6.9 yr; and eccentricities of 0.0 and 0.52 for datasets covering the time span from 1999 to 2018 and 2007 to 2012, respectively. Wittenmyer et al. (2013) ran N-body simulations for the two-planetary model and proved that it is unstable on a short timescale. Also, Hinse et al. (2014b) could not confirm the two-planetary model of Almeida et al. (2013), while their single LTTE solution was unconstrained due to a short time baseline. The ETV analysis of Nasiroglu et al. (2017) resulted in a single circumbinary brown dwarf with a of mass 15 MJ , a period of 9.9 yr, and an eccentricity of 0.175 for the time span 2007-2016. This solution is in close agreement with the results of Zhu et al. (2019) ( m 3 = 14.2 MJ , P 3 = 8 . 83 yr, e 3 = 0 . 12) and Nehir & Bulut (2022) ( m 3 = 13.2 MJ , P 3 = 8 . 83 yr, e 3 = 0 . 13) for datasets covering the time span 2008-2018 and 2007-2019, respectively. New data for the following years have indicated a divergence from these single LTTE term solutions. Wolf et al. (2021) concluded that at least two circumbinary bodies were necessary to explain the ETV data, as a single LTTE term with a period of 14 years was insu ffi cient to fit the data. Similarly, Pulley et al. (2022) reached the same conclusion, finding that an updated ETV diagram (2007-2021) could not be adequately explained by a quadratic ephemeris and a circumbinary companion in a near circular orbit ( e = 0 . 02) with a period of 7.65 years.  \nIn this study, we collect all available ETV data updated with new times of minima (TOMs) calculated from the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) LC data and new observations. Our objective was to search for the best-fitting curve of two LTTE terms in the ETV diagram by implementing a grid search optimization scheme alongside a dynamical stability analysis of N-body simulations, which we believe has not been attempted yet for this PCEB in this extent.  \nIn the following section, we review alternative explanations for the ETVs in PCEBs. Then, we describe the new TOMs and the compilation of all available ETV data in Section 3. In Section 4, we analyze our search for the best-fitting curve of two LTTE terms in the ETV diagram and provide the results for different models, alongside a dynamical stability analysis and alternative mechanisms. Our discussion of the results can be found in Section 5.", '2. Possible explanations for the ETVs in PCEBs': "The LTTE concept faces challenges not only in terms of dynamical models for hypothetical planetary systems but also with respect to optimization techniques and phenomena inherent in the binary. Specifically, the main datasets used for ETV analyses (i.e., eclipse timings) originate from diverse sources and methods, involve various software packages and algorithms, and are often accompanied by underestimated errors (Pribulla et al. 2012).  \nAlternative mechanisms have been proposed as explanations for the apparent ETV variations, such as mass transfer, apsidal motion (Lacy 1992), magnetic e ff ects (Applegate 1992; Völschow et al. 2016; Lanza 2020), and angular momentum loss through gravitational radiation (Paczy'nski 1967) or magnetic breaking (Rappaport et al. 1983) in one or both stars. Since sdB binaries are detached systems, mass transfer cannot occur. Apsidal motion would not be expected in PCEBs with α < 1 R ⊙ ;  \nthus the orbit can be considered circular (e = 0), and the rotation of the components synchronized with the orbit. Parsons et al. (2014) suggest that orbital eccentricities much less than 0.001 can introduce a maximum ETV of 5.5 s (for P ∼ 0 . 1 d), and thus it is considered a minor contributing factor to the observed ETVs. However, apsidal motion is the combined result of perturbations of a Keplerian orbit that include tidal and rotational contributions, a general relativity contribution, and possibly a thirdbody contribution. Analytical formulae can be found in Parsons et al. (2014); Bulut et al. (2017); Almeida et al. (2020); Dimo ff & Orosz (2023). Due to their proximity, PCEBs frequently undergo complex gravitational e ff ects that depart from idealized spherical models. Due to dynamic interactions between the stars as they evolve and lose their spherical symmetry, caused by factors such as rotation, magnetic fields, and tidal forces, precession of the orbits may result. In situations, where the gravitational term is not taken into account, the rate of apsidal motion can be used as a probe of the internal structure of the stars. At the same time, although it was not ruled out, Almeida et al. (2020) discovered that apsidal motion due to tidal and rotation e ff ects was less plausible than substellar component models in the case of the PCEB GKVir. Therefore, we do not take this e ff ect into further consideration in the current investigation.  \nBecause of the small masses, barely any recorded ETV signals may be attributed to gravitational radiation despite the short orbital periods. Considering the high surface gravities of the components, star winds should also be extremely weak and make minimal contributions. Given that the secular variations can occur at a timescale significantly longer than the baseline of data, it is very easy to confuse periodic changes with the former. The overall cyclical behavior of ETV cannot be explained by angular momentum loss (AML) from magnetic braking or gravitational radiation, which drives the binary to a constant reduction in the period without a quasi-cyclical component.  \nMagnetic activity cycles of the secondary components can imitate ETV behavior provided that the star has a su ffi cient energy budget (Applegate mechanism). According to the Applegate model (Applegate 1992) a magnetic active star in a binary can transfer angular momentum toward the surface from its interior and away from it during a complete magnetic cycle, leading to oscillations of the orbital period. The magnetic star's energy budget must be su ffi cient to adjust its quadrupole moment for this mechanism to work, resulting in a variable luminosity. Völschow et al. (2016) proposed a more accurate and realistic two-zone model where they treated a stellar density profile of the inner and outer regions for two di ff erent densities, including the change of the quadrupole moment both in the shell and the core, adding to previous improvements of the mechanism (Lanza & Rodonò 2004; Brinkworth et al. 2006; Lanza 2006; Tian et al. 2009). A novel model has been developed recently by Lanza (2020). It is based on the exchange of angular momentum between the binary orbital motion and the spin of the magnetically active secondary. This mechanism could help explain the observed ETVs and potentially get over the original Applegate mechanism's energy constraints. Although these alternatives often fail to explain the magnitude of the observed ETVs or to satisfy the required observational evidence of each mechanism (Pulley et al. 2022) in most cases, their contribution can become significant (RR Cae, DE CVn; Pulley et al. 2022).", '3. Data acquisition': "We compiled 305 eclipse minimum times from the literature as follows: 152 from Table 4 of Nasiroglu et al. (2017) (these in-  \nclude 18 by Wils et al. (2007), 32 by Beuermann et al. (2012a), nine by Kilkenny & Koen (2012), and ten byAlmeida et al. (2012)); 84 by Zhu et al. (2019), 41 by Wolf et al. (2021), 20 by Pulley et al. (2022), and eight by Nehir & Bulut (2022) (only the primaries). All TOMs are in BJD-TDB. Details of the observations (error, eclipse type, observatory) can be found in the papers mentioned. In addition, we updated the TOMs with our own observations carried out in 2017 using the 14-inch SchmidtCassegrain telescope at the Patras University 'Mythodea' Observatory (Papageorgiou & Christopoulou 2015a). Three new TOMs were calculated from TESS sector 51 with a 2-min cadence for the year 2022, available from the Mikulski Archive for Space Telescopes (MAST), 1 and two new TOMs were collected from the online B.R.N.O. 2 database during the years 2023 and 2024 (Table 1).  \nOur own observations were taken with an RC filter. The standard processes of bias subtraction, flat fielding, and dark subtraction were used to reduce the CCD frames. Di ff erential magnitudes were obtained by choosing a non-variable star in the field of view comparable to the target in brightness (GSC 00504 00639) and a comparison star (2MASS 1084555518). Image reduction and di ff erential photometry were performed using a fully automated pipeline (Papageorgiou & Christopoulou 2015b) that incorporates Pyraf (Science Software Branch at STScI 2012) and the Astrometry.net packages (Lang et al. 2010).  \nThe TESS LC in Fig. 1a was plotted with normalized SAP fluxes of the best quality (quality flag zero) data points using the Lightkurve Python package for Kepler and TESS data analysis Lightkurve Collaboration et al. (2018), from which we generated a phased-folded and binned LC (FBLC) for each time series segment (Fig. 1b). Employing a modified Kwee-Van Woerden (Deeg 2020) method for eclipse minimum timing with reliable error estimates, one primary TOM was calculated from each FBLC following the process of TOM determination described in Zervas et al. (2024). All the HJD times were transformed into TDB-based BJD ones by using the publicly available code of Eastman et al. (2010).  \nA total of 312 TOMs were used in the ETV analysis (Table B.1), excluding low-precision timings and using only the timings with errors smaller than ∼ 17sec (0.0002 d). Consistent with previous data categorizations (Beuermann et al. 2012a; Nasiroglu et al. 2017; Pulley et al. 2022), we established two datasets. The first one includes the less precise NSVS and All Sky Automated Survey (ASAS) data, o ff ering a broader time frame from 1999 to 2024 (dataset A), while the second one with 307 TOMs excludes these data, narrowing down the analysis to the years 2007 to 2024 (dataset B). Using the linear part of ephemeris from Pulley et al. (2022),  \nBJD = 2454274 . 20875(5) + 0 . 110374157(4) E , (1)  \nwe calculated and plotted the ETV diagram (Fig. 2) of each dataset.", '4.1. Analysis outline': "We adopted the formulation of Go'zdziewski et al. (2012) for the LTTE term:  \nLTTEp = Kp GLYPH<20> sin ω p cos( Ep ( t ) -ep ) + cos ω p q 1 -e 2 p sin Ep ( t ) GLYPH<21> ,  \nTable 1: New times of minimum light for NSVS 14256825.  \n(2)  \nwhere Kp , ep , ω p are the semi-amplitude of the LTTE signal, eccentricity, and argument of the pericenter, respectively, for body p ( p = b , c for inner and outer orbit, in accordance with the common convention). The orbital period Pp and the pericenter passage Tp depend on the eccentric anomaly Ep ( t ), and hence they do not explicitly appear in equation 2. This is a revised kinematic (Keplerian) formulation of the LTTE e ff ect (Irwin 1952) for multiple circumbinary companions, where the eclipse ephemerides is expressed with regard to Jacobi coordinates, with the origin at the center of mass (CM) of the binary.  \nAccording to Go'zdziewski et al. (2012), the osculating orbital elements and masses derived with this revised formulation best match the true N-body initial condition of the system with mutually interacting planets. However, these Keplerian orbits are only an approximation of the true (Newtonian) orbits in the case of more than two bodies, and as a result, one needs to determine whether the degree of approximation is significant. This discrepancy has been highlighted first in radial velocity studies of extrasolar planets (Laughlin & Chambers 2001; Go'zdziewski et al. 2003) and becomes substantial in cases of massive, strongly interacting planets close to low-order mean motion resonances (MMRs).  \nIn total, our least-squares fit involves twelve parameters: two for the linear ephemeris of the binary (the epoch T 0 and period Pbin ) and five orbital elements for each body ( ep , Pp , ω p , Tp , α p ), where α p is the semi-major axis and p = b , c for each LTTE orbit. We followed a grid search approach similar to that of Beuermann et al. (2013). Our optimization process consisted of two sets of runs in the eb , ec plane utilizing a combination of the Nelder-Mead Simplex and Levenberg-Marquardt algorithms using the lmfit Python package (Newville et al. 2023). The goal was to identify the best-fitting curve in the ETV diagram of models featuring two LTTE terms.  \nAdditionally, each solution that achieved a reduced chisquare ( χ 2 ν ) value below a preset limit (listed below) was integrated with the rebound N-body code (Rein & Liu 2012) using the ias 15 integrator (Rein & Spiegel 2015), requiring a lifetime τ = 1Myr and a MEGNO chaos indicator (Cincotta et al. 2003) value ⟨ Y ⟩ ∼ 2 as an additional acceptance criterion. Therefore, an acceptable model can be considered one that provides a good fit to the data and is secularly stable with non-chaotic orbital behavior.  \nThis optimization process was independently applied to two distinct ETV datasets (dataset A, dataset B), each covering a different time span. The initialization of each integration was set up with the transformation of the initial condition from the Jacobian kinematic frame to the N-body Cartesian osculating frame, with the origin at the CM, at the osculating epoch T 0. Furthermore, we used these initial conditions and refined them in terms of the exact N-body model. An iterative Levenberg-Marquardt process was implemented, and we examined how the χ 2 value of the fit  \nν  \nFigure 3 illustrates the di ff erence between the synthetic LTTE signals of a Jacobian (Keplerian) and an N-body (Newtonian) fit for dataset A. The maximum time di ff erence of 1.2 sec is comparable to the discrepancy of Keplerian and Newtonian two-planet models for the ETV diagram of the post-common envelope binary NN Ser, as has been shown in Marsh et al. (2014), while these di ff erences are about ten times larger for the cataclysmic binary HU Aqr (Go'zdziewski et al. 2015).  \n<!-- image -->  \nFig. 1: (a) TESS LC time series composed of three segments (sector 51) (b) phase-folded binned LCs for each segment.  \n<!-- image -->  \n<!-- image -->  \nFig. 2: ETV diagrams of NSVS 14256825. (a) ETV diagram of dataset A spanning 26 years (1999-2024). Archival data ∼ black open circles; NSVS + ASAS data ∼ blue open circles; Mythodea data ∼ orange open circles; TESS data ∼ green open circles; B.R.N.O. data ∼ red open circles. (b) ETV diagram of dataset B spanning 18 years (2007-2024). Color indexing is the same as in (a).  \n<!-- image -->  \ndepends on variations of all orbital elements (Laughlin & Chambers 2001). In that way, a fully self-consistent Newtonian fit to the ETV data could be made.  \nThe steps of the fitting process are outlined as follows:  \n- 1. Run 1: Grid search optimization in the eb , ec plane with fixed eccentricity values in the range [0, 0.98] and a step size of 0.01 (9801 models). The rest of the parameters were subsequently optimized, initiating from values chosen randomly within physically accepted limits.\n- 2. Run 2: The resulting solutions of Run 1 served as the initial points for a new optimization run. However, this time all the parameters, including eccentricities, were adjusted during the fitting process.\n- 3. N-body simulations for 1 Myr of the obtained Keplerian solutions of both runs while selectively excluding models that deviate beyond the 90% confidence level of the reduced chisquare ( χ 2 ν ) for two degrees of freedom in order to acquire well-fitted curves. Specifically, we considered solutions with χ 2 ν ≤ χ 2 ν, best + ∆ χ , where χ 2 ν, best corresponds to the chi-square value of the best-fitting curve and ∆ χ = 4.6.\n- 4. Self-consistent N-body (Newtonian) orbital fits initialized from each Keplerian fit of step 3 and subsequent integrations for 1Myr.  \nFig. 3: Di ff erences of the LTTE signals derived from a typical Jacobian fit and from respective osculating N-body models integrated numerically with the inferred initial condition at the osculating epoch T 0 = 54274 . 20875 BJD of the zeroth cycle, where the two models seem to agree.  \n<!-- image -->  \nIn the first optimization run (Run 1), the eccentricities are fixed, and the rest of the parameters are free to adjust, while in the second run (Run 2), all parameters are allowed to vary. Throughout all orbital simulations, we treated the central binary as a single object and considered only co-planar, edge-on cases. This is the simplest scenario, and it is commonly assumed in sta-  \nbility analyses of exoplanets, as it minimizes the masses of the planets relative to the binary, thereby promoting stability (Beuermann et al. 2012b; Horner et al. 2012; Marsh et al. 2014; Mai & Mutel 2022; Esmer et al. 2023). However, there are cases, such as the putative HU Aqr planetary system, that appear to be more complex (highly non-coplanar), resulting in stable orbits in contrast to the co-planar scenario (Go'zdziewski et al. 2015). Fundamental parameters of the binary are from Nehir & Bulut (2022) and were used for the rest of the analysis when needed.", '4.2. Period variation and orbital stability': 'Optimization results of Keplerian models are shown in Table 2 and Table 3 for each dataset and Runs 1 and 2, respectively. Each optimization run is characterized by three critical values: the best reduced chi-square value ( χ 2 ν, best ) of the best-fitting curve to ETV data and the minimum- ( χ 2 ν, min ) and the maximum-reduced chisquare ( χ 2 ν, max ) values among all stable solutions. The same description is followed in Tables 4, 5 for the Newtonian models. The reported parameter uncertainties correspond to the standard 1σ error, and their calculation resulted from the covariance matrix.  \nIt is notable that none of the best-fitting curves ( χ 2 ν, best ) are stable for more than a few hundred years. However, there are hundred of stable orbits ( τ = 1Myr, ⟨ Y ⟩ ∼ 2) within the 90% confidence level of each of the best reduced chi-square values. Illustrated in Fig. 4-Fig. A.4 are the lifetime distributions of eccentricities as well as the periods and the fitting curves corresponding to the three critical cases of Keplerian and Newtonian models for each dataset.  \nFor dataset A, 99 stable Keplerian orbits and 74 stable Newtonian orbits were found among 5838 models, each with a chisquare value falling within the 90% confidence level of the bestfitting curve ( χ 2 ν, best = 5 . 26) during optimization Run 1. Subsequently, following optimization Run 2, 153 stable Keplerian orbits and 93 stable Newtonian orbits were identified out of 8123 models, all meeting the same criteria.  \nSimilarly, in dataset B, after optimization Run 1, we found 198 stable Keplerian orbits and 171 stable Newtonian orbits among 6055 models, each with a chi-square value lying within the 90% confidence level of the best-fitting curve ( χ 2 ν, best = 5 . 76). The greatest number of stable solutions was discovered after optimization Run 2 of dataset B, revealing 454 stable Keplerian orbits and 423 stable Newtonian orbits among 8216 models, all having reduced chi-square values within the 90% confidence level of χ 2 ν, best = 5 . 26.  \nThe lifetime distributions of Pb , Pc (subplots b and d of Figs. 4-A.3) reveal the existence of a constrained inner orbit with period Pb = 7 yr and an unconstrained outer orbit with periods in the range of 3:1 to 7:1 MMR. The same conclusion about the constrained (unconstrained) inner (outer) orbits can be deduced by the histograms of the Keplerian (Figure A.5) and Newtonian (Figure A.6) models for both datasets.  \nBoth datasets seem to constrain the inner planet to a nearly circular orbit with a mean period of Pb ∼ 7 yr, while the outer companion has a circular orbit for either dataset and optimization run. When comparing the mean outer periods across datasets and optimization runs, a noticeable deviation of almost 27 years was observed between the two datasets of optimization Run 1 (Figure A.5d, Figure A.6d), while this discrepancy decreases to a deviation of at least seven years for the results of optimization Run 2 (Figure A.5h, Figure A.6h). Nevertheless, the well-defined peak in the outer period histograms of dataset A suggests a more con-  \nrained outer period with a mean value of Pc ∼ 23 yr, in contrast to the broader distributions of dataset B with a mean value of Pc ∼ 40 yr.  \nThe N-body osculating initial conditions derived from the kinematic (Keplerian) fits provide χ 2 ν values that are comparable with the Jacobian solutions. Figure 6 illustrates a kinematic model and its formal transformation to the N-body Cartesian osculating frame, from which Figure 3 was produced. Similarly, the refined self-consistent N-body solutions result in χ 2 ν values comparable to those of the Keplerian fits. The only exception is the χ 2 ν, max solutions of optimization Run 1 (Tables 2, 4), which is justified because the eccentricities of self-consistent N-body models oscillate, in contrast to their fixed values in Run 1. Therefore, the kinematic solutions and the transformed kinematic initial conditions are as consistent with the observations as the selfconsistent Newtonian fits.  \nWe investigated further the orbital stability of the eight critical solutions ( χ 2 ν, min , χ 2 ν, max of both datasets and all optimization runs) by integrating each one for 100 Myr. All orbits were found to be stable in this timescale, which is considered the lifetime of the PCEB phase of sdB binaries (King & Schenker 2002; Schreiber & Gänsicke 2003). In Fig. 7, we plot indicatively the orbital evolution of eccentricities and semi-major axes, initialized by a Keplerian solution (Table 2, dataset B, χ 2 ν, min ), over 100 Myr.  \nAdditionally, in Fig. 7c, we display the LTTE e ff ect of this indicative solution for the first 152 years, phase folded in the ETV diagram. Each synthetic curve corresponds to an integration time equal to the outer period (38 years). The black synthetic curve is the Keplerian (kinematic) fit, which appears to be in absolute agreement with the dynamical fit for the first 38 years (green curve). The rest of the synthetic curves are also in close agreement with the initial Keplerian fit. This allowed us to test the validity of the Keplerian fits to the data.  \nFinally, in Fig. 8 we plot the temporal evolution of the meanmotion resonant angles θ 1 = 3 λ o -λ i -2 ϖ i , θ 2 = 3 λ o -λ i -2 ϖ o , and θ 3 = 3 λ o -λ i -ϖ i -ϖ o (Beaugé et al. 2003; Michtchenko et al. 2008) for the first 10000 years of a stable solution close to 3:1 MMR (Table 2, χ 2 ν, min , dataset A). It is apparent that these critical angles librate around zero degrees, indicating a dynamical resonant system.', '4.3. Alternative explanations for ETV': 'We have investigated if the cyclic period modulation can also be driven by magnetic activity cycles (Applegate 1992) of the secondary component of NSVS 14256825, as it is a convective red dwarf star with an e ff ective temperature of 3077 K and a mass of 0 . 095 M ⊙ (Nehir & Bulut 2022). Using the online tool of Völschow et al. (2016) 3 , we calculated the energy demand ( ∆ E / E 2) for two di ff erent models: (i) the finite-shell two-zone model (Völschow et al. 2016) and (ii) the thin-shell model (Tian et al. 2009). We also calculated the requirements for the spinorbit coupling model (Lanza 2020). To explain the observed ETV modulation, the magnetic activity of the secondary must meet the criterion ( ∆ E / E 2 << 1) for the energy demand. All the calculated energy ratios for a period modulation of 7 to 48 yrs were found to be larger than the threshold. Thus, we investigated the requirements for the spin-orbit coupling model (ApplegateLanza) to account for residual changes in the O-C stable model fits. For the two datasets, A and B, these are 22s and 43s over  \nTable 2: Grid search optimization results of Keplerian two-planet LTTE fit after Run 1 of the fitting process for both datasets (A and B).  \nNotes. Columns χ 2 ν, best display the best-fitting curve parameters, whereas columns χ 2 ν, min and χ 2 ν, max contain parameters corresponding to the minimum- and maximum-reduced chi-square values of stable configurations, respectively. Digits in parentheses are for the uncertainty at the last significant place.  \nTable 3: Grid search optimization results of Keplerian two-planet LTTE fit after Run 2 of the fitting process for both datasets (A and B).  \nNotes. Columns χ 2 ν, best display the best-fitting curve parameters, whereas columns χ 2 ν, min and χ 2 ν, max contain parameters corresponding to the minimum- and maximum-reduced chi-square values of stable configurations, respectively. Digits in parentheses are for the uncertainty at the last significant place.  \nT = 18 yr and T = 26 yr, respectively. Following the rationale of the equations of Mai & Mutel (2022) and Lanza (2020), the change in rotational energy ∆ Erot ∼ -(2 -3) × 10 32 J is smaller than the maximum energy available (4 . 4 -6 . 3) × 10 32 J from luminosity changes over the timescale for O-C variability in the secondary star Emax = L 2 × T . Therefore, ∆ Erot / Emax is ∼ 0 . 5, and it could account for the residual changes in the ETV model fits. Given the spectral type M5-M8 for the secondary of NSVS 14256825, this is in agreement with the conclusion of Bours et al. (2016) that as magnetic activity decreases toward a later spectral type,  \nwe expect to drive much smaller O-C variations. Nevertheless it is still a high number, meaning that half of the energy generated by the secondary star should account for the non-axisymmetric gravitational quadrupole moment. We also calculated the tidal synchronization timescale following the assumptions and equations of Lanza (2020) and found ts = 3 . 4 -34 yr for the dimensionless tidal quality factor Q = 10 5 -10 6 that respectively characterizes the e ffi ciency of tidal energy dissipation in the active (secondary) component. As these times are similar to the observed O-C changes, the timescale for spin-orbit coupling may  \nTable 4: Self-consistent Newtonian two-planet orbital fits initialized from the results of optimization Run 1 for both datasets (A and B).  \nNotes. Columns χ 2 ν, best display the best-fitting curve parameters, whereas columns χ 2 ν, min and χ 2 ν, max contain parameters corresponding to the minimum- and maximum-reduced chi-square values of stable configurations, respectively. The term Mb , Mc is the resulting mean anomaly of the inner and outer orbits, with each configuration being set up in the N-body osculating frame at the osculating epoch T 0. Digits in parentheses are for the uncertainty at the last significant place.  \nTable 5: Self-consistent Newtonian two-planet orbital fits initialized from the results of optimization Run 2 for both datasets (A and B).  \nNotes. Columns χ 2 ν, best display the best-fitting curve parameters, whereas columns χ 2 ν, min and χ 2 ν, max contain parameters corresponding to the minimum- and maximum-reduced chi-square values of stable configurations, respectively. The term Mb , Mc is the resulting mean anomaly of the inner and outer orbits, with each configuration being set up in the N-body osculating frame at the osculating epoch T 0. Digits in parentheses are for the uncertainty at the last significant place.  \nbe determined by tidal synchronization instead of by magnetic activity timescales.', '5. Discussion and conclusions': "Weimplemented a grid search approach consisting of two sets of optimization runs (Run 1: for fixed eccentricities, Run 2: adjusting all resulting parameters from Run 1) in the eb , ec plane utiliz-  \ning a combination of the Nelder-Mead Simplex and LevenbergMarquardt algorithms. The goal was to find the best-fitting curve ( χ 2 ν, best ) in the ETV diagram of two datasets (dataset A, dataset B), implementing a revised kinematic (Keplerian) formulation of the LTTE e ff ect and self-consistent N-body (Newtonian) models. Additionally, we ran three-body simulations of the resulting solutions that lie within the 90% confidence interval of the best-  \nFig. 4: Top: Lifetime colormap plots in the (a) eb , ec and (b) Pb , Pc planes for 5838 Keplerian models (99 stable) within the 90% confidence level of χ 2 ν, best as they resulted from optimization Run 1 for dataset A. Bottom: Same as top but for 5838 Newtonian models (74 stable). The color coding displays the range of lifetimes from white to dark blue, while the orange and green crosses correspond to the positions of the minimum- and maximum-reduced chi-square value of the stable configurations. The best-fitting curve to the data is identified by the red cross, which is unstable in all cases. The MMR solutions have been labeled in the Pb , Pc plane for the nominal Pb = 7 yr.  \n<!-- image -->  \ntting curve, searching for stable ( τ = 1Myr) and non-chaotic ( ⟨ Y ⟩ ∼ 2) orbits.  \nAs a result, hundreds of stable orbits, with a lifetime of at least 1 Myr, were found for two circumbinary bodies of substellar mass within the 90% confidence interval of the best-fitting curve. However, none of the Keplerian and Newtonian χ 2 ν, best solutions were dynamically stable for more than a few hundred years. The histograms and lifetime distributions of inner and outer periods revealed the existence of a constrained inner orbit with period Pb = 7 yr, mass mb = 11 MJup , and an unconstrained outer orbit with a period in the range of 3:1 to 7:1 MMR ( Pc = 20 -50 yr) and a mass range of mc = 11 -70 MJup . Recent ETV studies of multiple-planet systems have reported that orbits close to low-order MMRs are strongly unstable (Horner et al. 2012; Go'zdziewski et al. 2015; Hinse et al. 2014a; Esmer et al. 2021) for timescales as short as thousands of years. The only exception based on ETV analysis is NY Vir (Esmer et al. 2023), with a period variation attributed to two putative Jovian planets close to 2:1 MMR, although this is uncertain due to the short time interval of the observations.  \nFurthermore, the observational timeline of ETV diagrams are usually narrow relative to the putative orbital periods, and as a consequence, the inferred masses of hypothetical planets reach the brown dwarf and the red dwarf limits (Hinse et al. 2012; Go'zdziewski et al. 2015). Such cases may introduce significant deviations between synthetic signals derived from the Keplerian and osculating Newtonian elements, leading to qual-  \ntatively di ff erent best-fitting configurations constrained by the available data. In the case of NN Ser, the Newtonian fits of Marsh et al. (2014) are in qualitative agreement with the Keplerian fits of Beuermann et al. (2013). However, the following analyses (Bours et al. 2016; Pulley et al. 2022; Özdönmez et al. 2023) of the system revealed that the initially stable system in 2:1 MMR does not fit the updated time measurements. Here, both datasets confirm the presence of a circumbinary planet with Pb ∼ 7 yr, but deviations in the outer period were observed: Pc ∼ 23 yr for dataset A and Pc ∼ 40 yr for dataset B. This discrepancy suggests that the extended timeline of dataset A, spanning 26 years, may provide better constraints on the orbital parameters compared to the narrower 18-year time frame of dataset B.  \nIn any case, the fitting curves of stable orbits are not perfectly fitted to the complete datasets in comparison to the bestfitting curves. The low precision NSVS and ASAS data seem to follow the general trend of strong peak-to-peak amplitude reduction over time, although the fitting curves of the minimum- and maximum-reduced chi-square of stable solutions do not fit well within their respective time frames (1999-2007). Moreover, both datasets exhibit visually similar fits for the minimum- ( χ 2 ν, min ) and maximum-reduced ( χ 2 ν, max ) chi-square of stable solutions in the 2007-2024 interval, but they do not fit perfectly through the latest data of 2021-2024. In general, the self-consistent N-body models account for mutual gravitational interactions between the circumbinary bodies, which could a ff ect the ETV trends. We showed that the temporal variation of the mid-eclipse times of  \nbest  \n(a) χ 2 ν, best (dynamically unstable)  \n<!-- image -->  \n(b) χ 2 ν, min (dynamically stable)(d) χ 2 ν, best (dynamically unstable)  \n<!-- image -->  \nFig. 5: Keplerian fits for both datasets [A: (a)-(f), B: (g)-(l)] and optimization runs [1: (a)-(c), (g)-(i) ,2: (d)-(f), (j)-(l)]. The bestfitting curves ( χ 2 ν, best ) result in dynamically unstable configurations, while the fitting curves of the minimum- and maximum-reduced chi-square value ( χ 2 ν, min , χ 2 ν, max ) define the χ 2 ν limits among all stable configurations. The parameters of each case can be found in Table 2 and Table 3.  \n<!-- image -->  \n(g) χ 2 ν, best (dynamically unstable)  \n<!-- image -->  \n<!-- image -->  \n(j)  \n2  \nχ  \nν,  \n(h) χ 2 ν, min (dynamically stable)  \n(i) χ 2 ν, max (dynamically stable)  \n<!-- image -->  \nν,  \nthe binary, as predicted by a best-fitting dynamical two-planet model, are in close agreement with the initial Keplerian fit. In addition, our Keplerian fits produce qualitatively the same results as the Newtonian ones. As a result, we attribute these ETV residuals to the sparse data of the corresponding time frames and / or additional mechanisms (Applegate-Lanza, spin-orbit coupling model). We applied the model of Lanza (2020), and although it depends on poorly constrained binary system parameters, such as the tidal energy dissipation and the internal mass distribution of the stars, we found that it may contribute to the observed ETV, apart from the potential circumbinary bodies, revealing the complex nature of the orbital structure of NSVS 14256825, as in  \nthe case of HS 0705 + 6700 (Mai & Mutel 2022) and V471 Tau (Kundra et al. 2022).  \nThe plentiful stable solutions that were found, coupled with their relatively well-fitted ETV curves, could serve as evidence supporting the existence of two circumbinary companions of substellar mass orbiting NSVS 14256825. We also tested the stability of the eight critical stable solutions ( χ 2 ν, min , χ 2 ν, max of both datasets and of both optimization runs) for an integration interval of 100 Myr. All of these solutions were found to be stable in this timescale, which is considered to be the lifetime of the PCEB phase of sdB binaries.  \n(l) χ 2 ν, max (dynamically stable)  \n<!-- image -->  \n(dynamically unstable)  \n(e) χ 2 ν, min (dynamically stable)  \n<!-- image -->  \n<!-- image -->  \n(c) χ 2 ν, max (dynamically stable)  \n<!-- image -->  \n(f) χ 2 ν, max (dynamically stable)  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFig. 6: ETV diagrams with synthetic curves. Left: Keplerian fit of two-planet Jacobian solution (Table 2, χ 2 ν, min , dataset A). Right: N-body synthetic curve for the osculating initial condition derived through the formal transformation of the Jacobian elements (Keplerian fit) to the N-body Cartesian osculating frame, centered at the CM of the binary.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFig. 7: Orbital evolution of (a) eccentricities and (b) semi-major axes for the inner and outer planets over 100 Myr (dataset B, Run 1, χ 2 ν, min ). (c) The temporal variation of the χ 2 ν, min fitting curve (dataset B, Run 1) over the first 152 years is represented by five sets of parameters, each corresponding to a complete outer orbit revolution (38 years). The black synthetic curve is the Keplerian (kinematic) fit from which this dynamical model was initialized, and it appears to be in absolute agreement with the first 38 years of orbital evolution (green curve) and in close agreement with the subsequent fits (blue and red curves).  \n<!-- image -->  \nWecould not suggest one conclusive model of two circumbinary planets since one cannot rely solely on the statistics of a reduced chi-square. Regarding the reduced chi-square for each dataset and optimization run, the eight critical stable solutions outline the extrema of a group of solutions. However, within these limits exist hundreds of solutions with various parameter values that are not centered around a best-fit solution. Dataset A seems to constrain the model parameters for both companions, but the fitting curves of the critical stable solutions are far from satisfactory for the earliest NSVS and ASAS data. On the other hand, dataset B provides more satisfactory fits, but there are cases where the mean outer period exceeds 30 years in both optimization runs and the inferred masses reach values greater than 30 MJup .  \nAmong the ten PCEBs exhibiting LTTE signals (for the complete list, refer to Pulley et al. 2022), half of them (HW Vir, NY Vir, NN Ser, HS 0705 + 6700, Kepler-451) exhibit two or more cyclic variations attributed to circumbinary planets. Their ETV and dynamical stability analyses share some common features, among which are the continuously updated LTTE solutions, peak-to-peak amplitude reduction (Pulley et al. 2022), and  \nstable regions of higher χ 2 than the best-fit solution (Potter et al. 2011; Dreizler et al. 2012; Esmer et al. 2022; Özdönmez et al. 2023). From our analysis, it seems that NSVS 14256825 is not an exception.  \nWe examined the potential for detecting the planets of ETV using Gaia (DR3; Gaia Collaboration et al. 2023) astrometry. The detection limit of an outer planet can be calculated by following the formalism of Sahlmann et al. (2015) and the NSVS 14256825 distance, 752.6 pc according to the Gaia DR3 parallax measurements. Among the stable models, using the masses and the periods of, for example, the χ 2 ν, max dataset B (Run2) solution ( mb = 11 . 12 MJ , Pb = 7 . 04 yr, mc = 8 . 57 MJ , Pc = 24 . 91 yr), the gravitational pull of the planet will displace the system's barycenter, with a semi-major axis of ab ∼ 19 mas and ac ∼ 34 mas, respectively. This represents the planet's astrometric signature, which to be detectable must be greater than the scan length accuracy per optical field pass σ f o v corresponding to the luminosity G (Perryman et al. 2014). For the NSVS 14256825 (G = 13.22) σ f o v ∼ 42 µ as, and thus Gaia can detect the planets of NSVS 14256825.  \nFig. 8: Temporal evolution of resonant angles (a) θ 1, (b) θ 2, and (c) θ 3 of a stable solution close to the 3:1 MMR (Table 2, χ 2 ν, min , dataset A) over the first 10000 years of integration.  \n<!-- image -->  \nAs recently proposed for the case of HW Vir by Baycroft et al. (2023), Gaia full-epoch astrometry (circa 2025) will provide an alternative method to ETVs for confirming circumbinary planets of PCEBs. In the meantime, it is essential to continue monitoring the system in order to expand the observational timeline with new data, providing further validation and refinement of the proposed solutions.  \nAcknowledgements. We thank the anonymous referee for numerous constructive comments and suggestions that greatly improved this paper.", 'References': "Almeida, L. A., Jablonski, F., & Rodrigues, C. V. 2013, ApJ, 766, 11  \n- Almeida, L. A., Jablonski, F., Tello, J., & Rodrigues, C. V. 2012, MNRAS, 423, 478\n- Almeida, L. A., Pereira, E. S., Borges, G. M., et al. 2020, MNRAS, 497, 4022 Applegate, J. H. 1992, ApJ, 385, 621\n- Baycroft, T. A., Triaud, A. H. M. J., & Kervella, P. 2023, MNRAS, 526, 2241 Bear, E. & Soker, N. 2014, MNRAS, 444, 1698\n- Beaugé, C., Ferraz-Mello, S., & Michtchenko, T. A. 2003, ApJ, 593, 1124\n- Beuermann, K., Breitenstein, P., Debski, B., et al. 2012a, A&A, 540, A8\n- Beuermann, K., Dreizler, S., & Hessman, F. V. 2013, A&A, 555, A133\n- Beuermann, K., Dreizler, S., Hessman, F. V., & Deller, J. 2012b, A&A, 543, A138\n- Borkovits, T., Rappaport, S., Hajdu, T., & Sztakovics, J. 2015, MNRAS, 448, 946\n- Bours, M. C. P., Marsh, T. R., Parsons, S. G., et al. 2016, MNRAS, 460, 3873\n- Brinkworth, C. S., Marsh, T. R., Dhillon, V. S., & Knigge, C. 2006, MNRAS, 365, 287\n- Brown-Sevilla, S. B., Nascimbeni, V., Borsato, L., et al. 2021, MNRAS, 506, 2122\n- Bulut, ˙ I., Bulut, A., & Demircan, O. 2017, MNRAS, 468, 3342\n- Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Physica D Nonlinear Phenomena, 182, 151\n- Deeg, H. J. 2020, Galaxies, 9, 1\n- Dimo ff , A. J. & Orosz, J. A. 2023, AJ, 166, 114\n- Dreizler, S., Beuermann, K., & Hesman, F. V. 2012, MEMSAI, 83, 498\n- Eastman, J., Siverd, R., & Gaudi, B. S. 2010, PASP, 122, 935\n- Esmer, E. M., Ba¸stürk, Ö., Hinse, T. C., Selam, S. O., & Correia, A. C. M. 2021, A&A, 648, A85\n- Esmer, E. M., Ba¸stürk, Ö., & Selam, S. O. 2023, MNRAS, 525, 6050\n- Esmer, E. M., Ba¸stürk, Ö., Selam, S. O., & Ali¸s, S. 2022, MNRAS, 511, 5207  \nGaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1  \nGo'zdziewski, K., Konacki, M., & Maciejewski, A. J. 2003, ApJ, 594, 1019  \nGo'zdziewski, K., Nasiroglu, I., Słowikowska, A., et al. 2012, MNRAS, 425, 930  \nGo'zdziewski, K., Słowikowska, A., Dimitrov, D., et al. 2015, MNRAS, 448,  \n1118  \n- Heber, U. 2016, PASP, 128, 082001\n- Hinse, T. C., Go'zdziewski, K., Lee, J. W., Haghighipour, N., & Lee, C.-U. 2012, AJ, 144, 34\n- Hinse, T. C., Horner, J., & Wittenmyer, R. A. 2014a, Journal of Astronomy and Space Sciences, 31, 187  \nHinse, T. C., Lee, J. W., Go'zdziewski, K., Horner, J., & Wittenmyer, R. A. 2014b, MNRAS, 438, 307  \n- Horner, J., Hinse, T. C., Wittenmyer, R. A., Marshall, J. P., & Tinney, C. G. 2012, MNRAS, 427, 2812  \nIrwin, J. B. 1952, ApJ, 116, 211  \nKilkenny, D. & Koen, C. 2012, MNRAS, 421, 3238  \n- King, A. R. & Schenker, K. 2002, in Astronomical Society of the Pacific Conference Series, Vol. 261, The Physics of Cataclysmic Variables and Related Objects, ed. B. T. Gänsicke, K. Beuermann, & K. Reinsch, 233  \nKundra, E., Hambálek, L'., Vanaverbeke, S., et al. 2022, MNRAS, 517, 5358 Lacy, C. H. S. 1992, AJ, 104, 2213  \n- Lang, D., Hogg, D. W., Mierle, K., Blanton, M., & Roweis, S. 2010, AJ, 139, 1782\n- Lanza, A. F. 2006, MNRAS, 369, 1773\n- Lanza, A. F. 2020, MNRAS, 491, 1820\n- Lanza, A. F. & Rodonò, M. 2004, Astronomische Nachrichten, 325, 393\n- Laughlin, G. & Chambers, J. E. 2001, ApJ, 551, L109\n- Lightkurve Collaboration, Cardoso, J. V. d. M., Hedges, C., et al. 2018, Lightkurve: Kepler and TESS time series analysis in Python, Astrophysics Source Code Library\n- Mai, X. & Mutel, R. L. 2022, MNRAS, 513, 2478\n- Marsh, T. R. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte, 96\n- Marsh, T. R., Parsons, S. G., Bours, M. C. P., et al. 2014, MNRAS, 437, 475  \nMichtchenko, T. A., Beaugé, C., & Ferraz-Mello, S. 2008, MNRAS, 387, 747  \n- Michtchenko, T. A., Ferraz-Mello, S., & Beaugé, C. 2011, Workshop Series of\n- the Asociacion Argentina de Astronomia, 3, 247\n- Nasiroglu, I., Go'zdziewski, K., Słowikowska, A., et al. 2017, AJ, 153, 137\n- Nehir, Ç. & Bulut, ˙ I. 2022, New A, 92, 101703\n- Newville, M., Otten, R., Nelson, A., et al. 2023, lmfit / lmfit-py: 1.2.2\n- Özdönmez, A., Er, H., & Nasiroglu, I. 2023, MNRAS, 526, 4725\n- Paczy'nski, B. 1967, Acta Astron., 17, 287\n- Papageorgiou, A. & Christopoulou, P. E. 2015a, AJ, 149, 168\n- Papageorgiou, A. & Christopoulou, P. E. 2015b, in Astronomical Society of the Pacific Conference Series, Vol. 496, Living Together: Planets, Host Stars and Binaries, ed. S. M. Rucinski, G. Torres, & M. Zejda, 181\n- Parsons, S. G., Marsh, T. R., Bours, M. C. P., et al. 2014, MNRAS, 438, L91  \nPerryman, M., Hartman, J., Bakos, G. Á., & Lindegren, L. 2014, ApJ, 797, 14  \n- Potter, S. B., Romero-Colmenero, E., Ramsay, G., et al. 2011, MNRAS, 416, 2202\n- Pribulla , T., Vaˇnko, M., Ammler-von Ei ff , M., et al. 2012, Astronomische Nachrichten, 333, 754\n- Pulley, D., Sharp, I. D., Mallett, J., & von Harrach, S. 2022, MNRAS, 514, 5725\n- Rappaport, S., Verbunt, F., & Joss, P. C. 1983, ApJ, 275, 713\n- Rein, H. & Liu, S. F. 2012, A&A, 537, A128\n- Rein, H. & Spiegel, D. S. 2015, MNRAS, 446, 1424\n- Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003\n- Sahlmann, J., Triaud, A. H. M. J., & Martin, D. V. 2015, MNRAS, 447, 287\n- Schreiber, M. R. & Gänsicke, B. T. 2003, A&A, 406, 305\n- Science Software Branch at STScI. 2012, PyRAF: Python alternative for IRAF Tian, Y. P., Xiang, F. Y., & Tao, X. 2009, Ap&SS, 319, 119\n- Völschow, M., Schleicher, D. R. G., Perdelwitz, V., & Banerjee, R. 2016, A&A, 587, A34\n- Wils, P., di Scala, G., & Otero, S. A. 2007, Information Bulletin on Variable Stars, 5800, 1  \nWittenmyer, R. A., Horner, J., & Marshall, J. P. 2013, MNRAS, 431, 2150 Wolf, M., Kuˇcáková, H., Zasche, P., et al. 2021, A&A, 647, A65 Wo'zniak, P. R., Williams, S. J., Vestrand, W. T., & Gupta, V. 2004, AJ, 128, 2965 Zervas, K., Christopoulou, P.-E., & Papageorgiou, A. 2024, ApJ, 961, 97 Zhu, L., Qian, S., Liu, L., et al. 2011, in Astronomical Society of the Pacific Conference Series, Vol. 451, 9th Pacific Rim Conference on Stellar Astrophysics, ed. S. Qain, K. Leung, L. Zhu, & S. Kwok, 155 Zhu, L.-Y., Qian, S.-B., Fernández Lajús, E., Wang, Z.-H., & Li, L.-J. 2019,  \nResearch in Astronomy and Astrophysics, 19, 134 Zorotovic, M. & Schreiber, M. R. 2013, A&A, 549, A95", 'Appendix A: Additional figures': 'Fig. A.1: Same as Figure 4 but for 8123 Keplerian models (153 stable) within the 90% confidence level of χ 2 ν, best as they resulted from optimization Run 2 for dataset A (top) and for 8123 Newtonian models (93 stable, bottom).  \n<!-- image -->  \nFig. A.2: Same as Figure 4 but for 6055 Keplerian models (198 stable) within the 90% confidence level of χ 2 ν, best as they resulted from optimization Run 1 for dataset B (top) and for 6055 Newtonian models (171 stable, bottom).  \n<!-- image -->  \nFig. A.3: Same as Figure 4 but for 8216 Keplerian models (454 stable) within the 90% confidence level of χ 2 ν, best as they resulted from optimization Run 2 for dataset B (top) and for 8216 Newtonian models (423 stable, bottom).  \n<!-- image -->  \n(a) χ 2 ν, best (dynamically unstable)  \n<!-- image -->  \n(d) χ 2 ν, best (dynamically unstable)  \n<!-- image -->  \n<!-- image -->  \n(g) χ 2 ν, best (dynamically unstable)  \n<!-- image -->  \n(j) χ 2 ν, best (dynamically unstable)  \n<!-- image -->  \n(l) χ 2 ν, max (dynamically stable)(b) χ 2 ν, min (dynamically stable)  \n<!-- image -->  \n(e) χ 2 ν, min (dynamically stable)  \n<!-- image -->  \n<!-- image -->  \nEpoch (cycles)  \n(h) χ 2 ν, min (dynamically stable)  \n(i) χ 2 ν, max (dynamically stable)Fig. A.4: Newtonian fits for both datasets [A: (a)-(f), B: (g)-(l)] and optimization runs [1: (a)-(c), (g)-(i) ,2: (d)-(f), (j)-(l)]. The best-fitting curves ( χ 2 ν, best ) result to dynamically unstable configurations, while the fitting curves of the minimum and maximum reduced chi-square value ( χ 2 ν, min , χ 2 ν, max ) define the χ 2 ν limits among all stable configurations. The parameters of each case can be found in Table 4 and Table 5.  \n<!-- image -->  \n<!-- image -->  \n(c) χ 2 ν, max (dynamically stable)  \n<!-- image -->  \n(f) χ 2 ν, max (dynamically stable)  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFig. A.5: Histograms resulted from optimization runs. Top: Normalized histograms of eb , ec , Pb , Pc for the Keplerian models as they resulted from the optimization Run 1 for both datasets A, B. The dashed lines indicate the center value of the bin with the maximum frequency. Bottom: Same as top panel but for optimization Run 2.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFig. A.6: Same as Figure A.5 but for the Newtonian models.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->', 'Appendix B: Times of minima of NSVS 14256825': 'Times of minima (1999-2024) of NSVS 14256825 (literature and new).  \nTable B.1: Times of minima (1999-2024) of NSVS 14256825 (literature and new).  \nTable B.1 - continued from previous page  \nNotes. Minima with ( ∗ ) are not included in dataset B.  \nReferences. (1) Wils et al. (2007); (2) Beuermann et al. (2012b) ; (3) Kilkenny & Koen (2012); (4) Almeida et al. (2013);(5) Nasiroglu et al. (2017); (6) Zhu et al. (2019); (7) Wolf et al. (2021); (8) Nehir & Bulut (2022); (9) Pulley et al. (2022); (10) This work (Table 1).'}