Investigation on the Orbital Period Variations of NN Ser: Implications for the Hypothetical Planets, the Applegate Mechanism and the Orbital Stability

We present 36 new mid-eclipse times obtained between 2017 and 2023 using the T100 telescope in Turkey, extending the time span of the 𝑂 − 𝐶 diagram to 25 years. Once again, these new observations show significant deviations from previous published models that were able to explain the observed variations of the binary period. We investigate two plausible explanations for this variability: the LTT effect due to the presence of one or two invisible low-mass (planetary) companion(s) in distant circumbinary orbits; other mechanisms, like e.g. the Applegate mechanism, associated with the magnetic cycles of the M-dwarf component of the WD+dM binary. Through MCMC analyzes, we demonstrate that the observed 𝑂 − 𝐶 variability can be explained by the presence of a planet with a minimum mass of ∼ 9 . 5 𝑀 𝐽 . This circumbinary planet orbits around the binary system with a period of about 19.5 years, maintaining a stable orbit for a timeline of 10 Myr. By adding a weak LTT signal from a secondary hypothetical planet we achieve statistically better results. However, the orbits of the bodies in a two-planet system remain stable only for a small range of the parameter space. The energy required to power the Applegate and other Applegate-like mechanisms is too high to explain the period variations observed. Thus, on the one hand there is substantial evidence supporting the existence of a planet in the NN Ser system, but on the other hand there are also compelling indications that cast doubt on the existence of a second hypothetical planet.


INTRODUCTION
NN Ser is a well-studied short-period eclipsing pre-cataclysmic binary system containing a very hot white dwarf (0.535  ⊙ ) and an M4 spectral type main-sequence dwarf (0.111  ⊙ ) (Haefner 1989;Wood & Marsh 1991;Catalan et al. 1994;Qian et al. 2009;Parsons et al. 2010a;Bours et al. 2016, and references therein).Its brightness decreases by ∼ 5.8 mag from  ∼ 17 mag (GAIA  = 16.53 mag) during the primary eclipse with an orbital period of ∼ 3.12 hr.In addition, the system shows an orbital period variation.Haefner (1989) investigated the orbital period variation of the eclipsing binary NN Ser, and updated the orbital parameters in 2004 (Haefner et al. 2004).Brinkworth et al. (2006) reported that the orbital period variation of the system may arise from either angular momentum loss through magnetic braking or the existence of a third object.The variation of the orbital period has been mostly explained by the light travel time (LTT) effect only, which results from the presence of one or two planets (Chen 2009; Qian et al. 2009;Beuermann et al. 2010;Horner et al. 2012a;Beuermann, Dreizler, & Hessman 2013;Mustill et al. 2013;Bours et al. 2016).However, other mechanisms like magnetic cycle, magnetic breaking, gravitational radiation have been dismissed as a plausible explanation.According to the report by ★ E-mail: aykut.ozdonmez@atauni.edu.trBeuermann, Dreizler, & Hessman (2013), the two-planet model was dynamically stable for a lifetime greater than 100 Myr utilising the system parameters, i.e. the orbital periods of 7.65 and 15.47 years, and the minimum masses of 6.97   and 1.73   for the inner and outer planet, respectively.In contrast to the abovementioned study, Mustill et al. (2013) discussed the orbital evolution of NN Ser, and demonstrated that the orbital configuration including two planets in the binary system evolves dynamically unstable on short timescales.Recent studies on the orbital period variation and orbital stability suggested that the stable system configuration including two planets consistent with the  −  diagram can be hypothetically constructed even for timescales larger than 1 Myr (Marsh et al. 2014;Bours et al. 2016).Navarrete et al. (2018) reported that the Applegate mechanism could cause additional fluctuations in eclipsing time variations.The recent observations since 2016 are inconsistent with all existing models, contrary to the predictions of previous models in the literature (see Fig 5 .in Pulley et al. 2022).Hardy et al. (2016) found convincing evidence for the existence of dust around NN Ser.Although this finding alone does not provide conclusive proof for the presence of planet(s) around NN Ser, it significantly enhances the possibility of planet formation from common-envelope material within a 'second-generation' scenario.In addition, a resonant interaction between the binary and the circumbinary disc may cause the loss of angular momentum leading to a slow change in the binary's orbital period (Artymowicz et al. 1991;Kashi & Soker 2011;Völschow, Banerjee, & Hessman 2014;Schleicher et al. 2015;Chen & Podsiadlowski 2017).
This study makes an observational contribution to the literature by presenting 36 new mid-eclipse times obtained between 2017 and 2023.The  −  diagram covers a period of 25 years by adding these eclipse times.The details of the observations, modelling primary eclipse light curves, and the compiled mid-eclipse times are given in Section 2. This study presents a comprehensive investigation of the variation in orbital period through analysis of the  −  diagram, orbital stability and the Applegate mechanism.The analysis methods and results are presented in Section 3, and further discussed in Section 4.

TIMING DATA
We performed optical photometric observations of NN Ser between Feb 25, 2017 and Feb 18, 2023 using a SI 1100 CCD camera (field of view of 21.5 ′ × 21.5 ′ and pixel scale of 0.31 arcsec pixel −1 ) with a 4 K × 4 K chip attached to the 1 m telescope of TUG T100 (TÜBİTAK National Observatory) in Antalya, Turkey.During the observations, the readout time was decreased to 2-3 seconds by using the sub-array mode (300 × 300 pixel size) together with the binning 2 × 2 mode of the camera.Thus, the time resolution of the light curves (LC) was increased.We did not use any filter during the observations in order to obtain high counts within an exposure time of 5-10 seconds.The CCD images were reduced using standard procedures, i.e. bias subtraction, flat fielding, and cosmic ray correction.Aperture photometry was performed on the reduced CCD images with the same methodology as in Er, Özdönmez, & Nasiroglu (2021).
The system becomes fainter about 5.8 mag, reaching ∼ 23 mag during the primary eclipse (Haefner et al. 2004).The brightness at mid-eclipse is lower than the 1.0m telescope limit.Thus, the base magnitude of the eclipse was set to be 5.8 mag lower than the median magnitude outside the eclipse.The errors of the base magnitudes were determined by calculating the median values of the errors observed outside of the eclipse, which were then adjusted by adding three times the standard deviation.To determine the mid-eclipse time from the primary eclipse LC, we first divided the light curve into two parts; ingress and egress.Then, we modelled these light curve parts with two Boltzmann functions.As shown in Fig. 1, the mid-eclipse time was calculated as the mean value between the mid-ingress and mid-egress times.The error in the mid-eclipse time is the total errors of the uncertainties of the mid-ingress and mid-egress times acquired from the models.Six out of 36 primary eclipse light curves do not include either ingress or egress parts.For these incomplete LCs, the mid-eclipse times were determined by adding half of the eclipse duration to the mid-ingress times or subtracting from the mid-egress times.The eclipse duration was calculated to be 4.2 ± 0.35 minutes, which represents the average time difference between the mid-times of ingress and egress for 30 complete light curves.
The mid-eclipse times presented in this study and those compiled from the literature are listed in Table 1.It should be noted that the mid-eclipse times obtained from the incomplete light curves are labelled as "T100-h".In Table 1, we present a total of 223 mid-eclipse times including 36 new primary mid-eclipse times obtained from our observations between 2017 and 2023.Beuermann et al. (2010) reanalyzed the mid-eclipse time calculations and made corrections to some of the values from the previous results.We used these revised mid-eclipse times and their statistical errors from Beuermann et al. (2010).On the other hand, all mid-eclipse times presented by Parsons  2014) are obtained from the secondary eclipse light-curves, but it was concluded that these times are not much affected by the apsidal precession, and hence they won't significantly differ from the general trend of the  −  diagram.These mid-eclipse times weren't used by Bours et al. (2016) in their analysis of the orbital period variation.Besides, Bours et al. (2016) didn't list the mid-eclipse times they used, so we couldn't compile the mid-eclipse times from their observations.For the the mid-eclipse times without any published error, we adopted a 30 s uncertainty, which is the largest data error as listed.The largest and average errors of mid-eclipse times obtained from our observations are about 10 s and 3 s, respectively.We note that the mid-eclipse times obtained from the LCs using telescopes larger or smaller than 2m diameter have approximately the same average errors.

ANALYSIS AND RESULTS
The results are presented in three sub-sections outlining the methods of analysis used; i.e. modelling the  −  diagram using the Markov-Chain Monte Carlo (MCMC) technique, the orbital stability analysis by N-body integration, and searching for magnetic activity by testing the Applegate and Applegate-like mechanisms, respectively.We have investigated these topics in depth to get a better understanding of the underlying mechanisms behind the orbital period variation and the structure of the NN Ser system.

Modelling O-C times
The mid-eclipse time variations caused by a planet similar in size to Jovian are relatively small, a few tens of seconds.It is expected that this variation is significantly larger than both observational and systematic timing uncertainties.One way to minimize bias, which may be caused by a eventual scattering, is to restrict the use of timing data considering their uncertainty and the measurement method.For this reason, two types of data sets (A and B) were used in our study similar to those used in the previous studies of NN Ser.Data Set A (hereafter DS-A) includes all 223 mid-eclipse times available to date in the literature and our new measurements listed in Table 1.
In Data Set B (DS-B), we excluded both the measurements having uncertainties larger than 2 seconds and those that were obtained from secondary eclipse.DS-B consists of 97 mid-eclipse times.The average errors of mid-eclipse times for DS-A/B are 2.58 and 0.52 s, respectively.
The presence of the LTT effect can result in significant sinusoidal deviations in the orbital period of a binary system.Consequently, the linear ephemeris should be continually revised to maintain precision over time.Parameters of the linear ephemeris for DS-A are obtained as follows.
Here,   ℎ is the linear ephemeris,  0 is the initial ephemeris at the zero cycle (i.e. = 0), and   is the orbital period of the binary system.The obtained parameters of linear ephemeris from DS-B have very similar values as the parameters in Eq. 1.It also gives the same cycle (i.e.no shift) relative to the commonly used linear ephemeris given by Beuermann et al. (2010).
The  −  diagrams of DS-A/B are constructed by this linear ephemeris (see Figure 2).
To test the presence of hypothetical bodies causing the orbital period variation, we used two models: Model I consisting of a quadratic term () and a LTT term ( 1 ) of a 3rd body for a one-planet system is described through Model II additionally includes an another LTT term ( 2 ) arising from 4th body in a two-planet system; i.e. (Model I +  2 ).The parameter  refers to a quadratic term which is computed as  /2, where  represents the derivative of the period over time (i.e., dP/dt).  represents the LTT term of the ith body as defined in Irwin (1952).
For the   , we used a modified formulation described by Goździewski et al. (2012) of the form.
Here, the semi-amplitude of the LTT signal of the To calculate the mass (  ) of the ith body, we solved the mass function.
Here,  12 sin   is the projected semi-major axis of the binary system around the barycentre of the planetary system, and   is the inclination of the th body's orbit.For the total mass of the binary, we adopted   = 0.646  ⊙ determined by Parsons et al. (2010a).
In order to express the eclipse time variability, we used the identical fitting process, i.e. using a Genetic algorithm (GA) and MCMC methodology based on a likelihood function (L), as described in our previous study (Er, Özdönmez, & Nasiroglu 2021).Uniform prior samples of all free parameters have been randomly assigned within the specified ranges; i.e., ,   ,   ,  0, ,   > 0 days,   ,    [−0.75, +0.75],    [0.08, 0.15] days and Δ 0  [−10, +10].The free parameter   in units of days is included in the L function to account for systematic uncertainties (  ), and it scales the raw uncertainties of eclipsing times (  ) in quadrature; such that We first obtained most plausible initial parameters with GA (Charbonneau 1995), and then we run MCMC to sample the posterior distribution by using the affine-invariant ensemble sampler implementation of the emcee package (Goodman & Weare 2010) provided by Foreman-Mackey et al. (2013).For MCMC, we used 512 initial conditions (walkers) and assessed how these walkers behaved over 50,000 and 120,000 steps in chains of each independent variable both in Model I and II, respectively.The best-fitting parameters and their associated uncertainties are determined by evaluating the 16th, 50th, and 84th percentiles of the marginalized distributions of maximized L. The most-plausible parameters for Model I -II are listed in Table 2 and 3, respectively.These results indicate that the bodies in the models are all Joviantype planets due to their masses (see Table 2-3).The 1D distributions of all models utilized on both DS-A/B exhibit a distinct, prominent peak.The samples also demonstrate a relatively uniform 2D distribution centered around a particular solution.In order to conserve space, the corner plots representing the 1-D and 2-D posterior probability distributions of the system parameters sampled by MCMC are shown only for Model I using DS-A and for Model II using DS-B in Fig. A1 and Fig. A2, respectively.

Testing Stability of Orbits
To examine the orbital stability of NN Ser for our multi-component models, we used the N-body orbital integration package of RE-BOUND1 (Rein & Liu 2012) including a Mean Exponential Growth factor of Nearby Orbits (MEGNO, Cincotta & Simó 2000) indicator and a Wisdom-Holman symplectic integrator (WHFast, Rein & Tamayo 2015).REBOUND utilizes N-body integration to simulate the motion of celestial objects such as planets and stars under the influence of gravity.It provides two significant results in this study: (1) The MEGNO chaotic parameter surface maps: It provides a MEGNO indicator <  > by testing the chaotic behaviour of the system components for a range of both semi-major axis and eccentricity values over a period of given time.The MEGNO indicator <  > indicates whether the orbits are chaotic or not.For the given initial conditions and time, if <  > is ≤ 2, the system will remain stable.Values of <  > greater than 2 indicate a chaotic (unstable) orbital configuration.In such a system, the resonance term between the components is high enough to deform the orbits, but it still takes some time for the orbital parameters to diffuse to orbit-crossing values.If at least one particle is ejected or collides, a value of 10 is assigned to <  >.
(2) Orbital stability timeline: It integrates the orbits for a given time, and It demonstrates the variations of the orbital parameters like the semi-major axis and eccentricity as a function of time.This can be used to determine how the planets interact, when a planet will escape the system or collide, and whether the orbit will remain stable for a given period of time.
For both simulations, the central binary star was treated as a single mass (  = 0.646  ⊙ ), and all orbits were restricted to co-planar.We set the timestep to be 0.13% of the shortest orbital period of the inner planet, i.e. 0.01 yr (see Rein & Tamayo 2015, for some discussion on the timestep choices).It is assumed that the limit distance for escaping from the system is 20 AU.For these conditions and using the component parameters listed in Table 2 -3, we performed the dynamical stability simulations to obtain both the MEGNO values for 10 Myr and the orbital stability timeline at most 10 Myr.
In the case of the MEGNO stability map of Model I, we integrated a one-planet system, and varied the semi-major axis  3 and the eccentricity  3 of the planet.For Model II, we individually evaluated  3,4 and  3,4 of both the outer and inner planet in a binary system with two planets.Our simulated dynamical result shows that the one-planet system constructed for the best-fit parameters of Models I remains stable at least 10 Myr.In addition, both MEGNO chaos parameters using the component parameters with their errors for DS-A/B remain in the stable region, as expected for a one-planet system.The existence of a second planet in the system (i.e.Model II) makes the interaction of the components more chaotic.Since the uncertainties of the system parameters determined by Model II using DS-A are very high, the result of the stability analysis predicting unstable orbital configuration is doubtful.The chaotic variation of the  4 and  4 parameters of the inner planet with the destabilisation of the orbit corrupts the configuration of the system within only 2000 yr.Therefore, the plots are not shown for the abovementioned models.
For Model II DS-B, we have derived the MEGNO stability map over a time scale of 10 Myr using the orbital solution of the best-fitting parameters (Fig. 3).In particular, the boundaries of the stable region look rather fractal.The MEGNO chaos parameter is ≤ 2 for the bestfitting parameters of both planets, implying non-chaotic behaviour (i.e. a stable orbit).However, especially for some parameter values of the inner planet within the errors, the orbital stability is disrupted and the inner planet has escaped or collided within 10 Myr.In order to determine the time scale on which orbits can persist undisrupted for up to 10 Myr, we obtained the orbital stability timeline of the systems constructed by randomly selecting 200 models from the MCMC chain that also explain the  −  diagram.The distribution of both the semi-major axis and the eccentricities of the planets in these 200 systems is shown in Fig. 4, considering their stability time scale.It reveals that 21, 36 and 5 system configurations are corrupted  Wood & Marsh (1991).In this study, the mid-times were obtained using the T100 telescope (referred to as T100), while the times obtained from the incomplete light curves are denoted as T100-h.The RMS values are also given for each  −  plot.
within 1, 1-5 and 5-10 Myr due to unstable orbits, respectively.However, in 138 out of 200 system configurations, the orbits remain stable over a longer time scale of more than 10 million years.In the system configurations constructed from the MCMC chain models that are stable for more than 10 Myr, the distances between the inner and outer planets are larger, while the eccentricities of the massive outer planets, which have more influence on the resonance factor, are smaller.

Magnetic mechanism
The presence of substellar components in binaries may not be certain since cyclical variations of the binary orbital period may also be caused by magnetic mechanisms first proposed by Applegate (1992).It was suggested that solar-like magnetic cycles in the secondary star can lead to changes in the angular momentum of the system, resulting in a cyclic modulation of the orbital period.The Applegate mechanism has been revised and improved by analyses of the magnetic mechanisms and their effects on orbital period variations by other researchers.Tian, Xiang, & Tao ( 2009) developed an approximate equation by integrating the concepts of the Applegate mechanism and an improved model that includes changes in the quadrupole moment due to rotational and magnetic energy (Lanza, Rodono, & Rosner 1998).This formula establishes the relationship between the energy required to drive the magnetic mechanism and the observed eclipse time variations.By adopting the formulation of Brinkworth et al. (2006) in different approximations, Völschow et al. (2016) produced the "two-zone" model, assuming different densities in the shell and core, as well as a detailed numerical model where the framework is applied to realistic stellar density profiles.This model requires far higher energy to generate the Applegate mechanism.Lanza (2020) took a different approach to the magnetic mechanism by introducing a new model.This model is based on the coupling of the spin of the active star with the orbital motion of the binary system.The coupling is directly influenced by the non-axisymmetric stellar quadrupole moment, rather than by tides on much longer timescales.This cyclic exchange of angular momentum between the stellar spin and the orbit leads to the modulation of the orbital period.This proposed mechanism requires significantly less energy than previously proposed models.We calculated the required energy Δ as a fraction of the available energy in the magnetically active secondary star,   , to generate the corresponding magnetic mechanism from the formulations of these three studies, which have different perspectives on the magnetic mechanism.For the calculations, we used the system parameters listed in Table 2-3 and the astrophysical parameters of NN Ser in Völschow et al. (2016).It should be noted that the energy ratios have been calculated only for the substellar components that have a LTT signal with a low amplitude (see Table 4).In all Applegate-like models, the magnetic mechanism is sufficient to explain the orbital period variation when the ratio Δ/  < 1.

DISCUSSION AND CONCLUSION
We compiled 223 mid-eclipse times of NN Ser including 187 previously published in the literature, covering a time span of 25 years.The timescale of the periodic variation in the  −  does indeed provide a clue to the mechanisms that might exist.For example, magnetic braking and gravitational radiation occur on large timescales, such as 10 8 − 10 9 years, because they evolve slowly and steadily.Since the timescale of the periodic variation in the  −  diagram of NN Ser is less than ∼ 25 years, it is highly unlikely that magnetic braking and gravitational radiation from the  −  diagram can be detected.However, the presence of the LTT effect or the Applegate mechanism can be observed on shorter timescales of a few years to decades.We modelled the timing variation in the most-recent − diagram by considering quadratic ephemeris and possible LTT signals using MCMC.We found that the LTT signal in the one-planet model originates from a planet more massive (  ∼ 9.4) than those reported in the literature, so most of the model parameters differ from the previous studies (Pulley et al. 2022, references therein).The system parameters obtained from the solution of a one-planet system (i.e.Model I) using both DS-A/B are consistent with each other.However, the use of DS-B instead of DS-A for Model I reduces the systematic uncertainty from 2.80 to 1.52 seconds.While the root mean square (RMS) value is 4.21 seconds for the Model I obtained from DS-A, the RMS value of 3.83 seconds is calculated for the model using DS-B.As data with large errors are removed from DS-B, both the average error of the timing data and the number of scattered data from the model are reduced, resulting in a slight decrease in the systematic uncertainty and the RMS value (see Table 2-3).It can be concluded that the dataset used for the solution of the hypothetical single-planet system has no significant alteration on the results.This conclusion is not valid for the parameters of the two-planet solution (i.e.Model II).The reasons are as follows: 1) System parameters for Model II using DS-A have very large errors.Moreover, the errors of some parameters are close to their parameter values.2) The semiamplitude of LTT term of the 4th body when using DS-A is close to the RMS value of the  −  and the average error of the timing data.3) Both the RMS value and the systematic uncertainty are not significantly reduced when the additional LTT term of the 4th body is added, especially for DS-A.The aforementioned problems are less apparent when DS-B is used to derive the parameter solution for the two-planet system.Therefore, the presence of the LTT signal of the 4th body should be searched only with DS-B.
Using DS-B in Model II produces the best statistical results to explain the orbital period variation of NN Ser.Reaching a conclusion based on statistical results alone may lead to astrophysical fault because the orbits may be highly unstable on short timescales and/or several other physical mechanisms can contribute to the timing variations.When we analyse the orbital stability of the system, a one-planet system is stable throughout almost the entire evolutionary timeline.However, when a second planet is included into the system, the interaction between the components becomes more chaotic, amplifying the resonant terms.Based on our analyses, it has been determined that even minor alterations in the parameters of the two-planet system can disrupt its orbital stability.For 20.5% of the orbital configurations formed from the MCMC chain models, orbital stability is maintained only for a duration ranging from 1 to 10 Myr.In addition, a two-planet system constructed using 69% of the parameter solutions remains stable throughout a significantly longer time scale of 10 Myr.These results are also consistent with the MEGNO stability map obtained for a two-planet system.Changing the parameters within their respective error ranges in the MEGNO map derived for 10 Myr causes the transition of the planetary orbits from stable to disrupted.It is difficult to make strong statement on the orbital stability of NN Ser when only the best-fitting model is stable and most around it in the uncertainty range is unstable.
To examine the possibility of cyclical magnetic activity playing a role in generating a hypothetical signal instead of the LTT signal, we calculated the energy ratios Δ/  for the smaller LTT signal in our models by using three different formulations of the magnetic mechanism.Based on the formulation proposed by Tian, Xiang, & Tao (2009), the ratio Δ/  was calculated to be approximately 1 for the secondary (inner) planet in both our models DS-A/B.This value suggests that the magnetic mechanism could potentially contribute to fluctuations in the eclipsing time variations.In all the other cases, the ratio Δ/  was found to be significantly greater than 1 (see Table 4).Consequently, the Applegate mechanism alone is insufficient to explain the observed variations of the binary period.The existence of planet(s) is also supported by the detection of dust around NN Ser (Hardy et al. 2016).The investigated signal of the 4th planet may also potentially be attributed to a resonant interaction between the binary system and the circumbinary disc.
Based on our models, we predict a decline phase of cyclical variation in the − diagram expected to occur in 2024-2025, particularly for Model II.However, the upward trend resulting from the quadratic ephemeris is expected to continue steadily for all models.To validate our predictions, further observations in the upcoming years will be essential.

Figure 1 .
Figure1.The light curve of NN Ser obtained from observations of TUG T100 at Feb 18, 2023 given as an example.The ingress and egress parts of the light curves were fitted individually with a Boltzmann function, as described in Sec. 2. In the plot, the fits are represented by red and blue dashed lines for ingress and egress parts, respectively, and the obtained mid-times are denoted by colored "X".The vertical blue dotted lines are the limits of the ingress and egress parts of the LC.

Figure 2 .
Figure 2.  −  diagrams and their residuals for Model I (left panels) and Model II (right panels), using DS-A (upper panels) and DS-B (lower panels).Linear ephemeris times are calculated using equation 1.The solid black lines indicate the best-fit model, which includes the LTT term of the hypothetical circumbinary Jovian planet(s) that is represented by colored dotted lines as well as a quadratic trend denoted by different colored dotted lines.The shaded gray area represents the ±3 posterior spread which is calculated from 1000 randomly selected parameter samples from the MCMC posterior.All data is represented with distinct colors and labeled with corresponding abbreviations, which serve as references; B10 Beuermann et al. (2010); B13 Beuermann, Dreizler, & Hessman (2013); F20 Faillance et al. (2020); H04 Haefner et al. (2004); H89 Haefner (1989); M14 Marsh et al. (2014); P10 Parsons et al. (2010b); P22 Pulley et al. (2022); PM02 Pigulski & Michalska (2002); Q09 Qian et al. (2009); WM91Wood & Marsh (1991).In this study, the mid-times were obtained using the T100 telescope (referred to as T100), while the times obtained from the incomplete light curves are denoted as T100-h.The RMS values are also given for each  −  plot.

Figure 3 .
Figure3.MEGNO chaos parameter surface maps were generated for a duration of 10 Myr using the two-planet solutions from Model II DS-B, covering a range of eccentricities and semi-major axis.Specifically, the analysis was conducted for (a) the 3rd body, i.e., the outer planet, and (b) the 4th body, i.e., the inner planet.The white circles on the maps represent the best-fit model parameters, accompanied by their corresponding errors.

Figure 4 .
Figure 4.The semi-major axis and eccentricities of the planets in the system configurations constructed by randomly selecting 200 models from the MCMC chains of Model II DS-B.The colour bar indicates the stability time scale for each pair of values.

Figure
Figure A1.1-D and 2-D projections of the posterior probability distributions of the free parameters inferred from the Model I for Data set A. The samples are thinned by selecting one on every 100 samples after removing burnout sample.Δ represent differences between posterior and calculated values of the parameter.Contours are for the 16th, 50th, and 84th percentile of samples in the posterior distribution.This figure is derived using corner.py(Foreman-Mackey 2016).

Figure
Figure A2.1-D and 2-D projections of the posterior probability distributions of the free parameters inferred from the Model II for Data set B. The other descriptions are the same as in Fig. A1

Table 1 .
The Mid-eclipse Times of NN Ser, its error and references.
Nasiroglu et al. (2017)2)   .The other orbital parameters of the th body are eccentricity (  ), longitude of pericenter (  ), orbital period (  ), and time of pericenter passage ( 0, ).To avoid weakly constrained   and   for quasi-circular and moderately eccentric orbits, we used Poincare elements; i.e.  ≡     and  ≡     .The eccentric anomaly of   includes the orbital period of the th body,   , and the time of pericenter passage,  0, .Further information can be found inGoździewski et al. (2012)andNasiroglu et al. (2017).

Table 2 .
The most plausible Model I and system parameters obtained by using DS-A/B, along with the corresponding RMS values of  −  times.

Table 3 .
The most plausible Model II and system parameters obtained by using DS-A/B, along with the corresponding RMS values of  −  times.