Photo-dynamical Analysis of Circumbinary Multi-planet system TOI-1338: a Fully Coplanar Configuration with a Puffy Planet

TOI-1338 is the first circumbinary planet system discovered by TESS. It has one transiting planet at P$\sim$95 day and an outer non-transiting planet at P$\sim$215 day complemented by RV observation. Here we present a global photo-dynamical modeling of the TOI-1338 system that self-consistently accounts for the mutual gravitational interactions between all known bodies in the system. As a result, the three-dimensional architecture of the system can be established by comparing the model with additional data from TESS Extended Mission and published HARPS/ESPRESSO radial velocity data. We report an inconsistency of binary RV signal between HARPS and ESPRESSO, which could be due to the contamination of the secondary star. According to stability analysis, the RV data via ESPRESSO is preferred. Our results are summarized as follows: (1) the inner transiting planet is extremely coplanar to the binary plane $\Delta I_b \sim 0.12 ^\circ$, making it a permanently transiting circumbinary planet at any nodal precession phases. We updated the future transit ephemerides with improved precisions. (2) The outer planet, despite its non-transiting nature, is also coplanar with the binary plane by $\Delta I_c=9.1^{+6.0 \circ}_{-4.8}$ (22$^\circ$ for 99\% upper limit). (3) The inner planet could have a density of $0.137 \pm 0.026$ g/cm$^{-3}$. With a TESS magnitude of 11.45, TOI-1338 b is an optimal circumbinary planet for ground-based follow-up and transit spectroscopy.


INTRODUCTION
Planets orbiting around both of the stars in binary systems are called the circumbinary planets (CBPs).Over the years the detections of CBPs mainly came from transit surveys like Kepler and TESS (e.g., Doyle et al. 2011a;Welsh et al. 2012Welsh et al. , 2015;;Kostov et al. 2020Kostov et al. , 2021)).Subject to the observational bias of the transit method, most confirmed CBPs are coplanar with their host binary orbital planes within ∼ 4.5 • and close-in (with semi-major axis around 1-2 times of system inner stability limit).Knowing the inclination distribution is essential to understanding the population and occurrence rates of the CBP population (Li et al. 2016;Martin 2019).More faraway and misaligned CBPs are harder to detect with transit photometry, as it will produce irregular transit patterns (Martin & Triaud 2014;Chen & Kipping 2021).However, CBPs elude transiting might dynamically perturb the binary and showcase eclipsing binary variations (ETVs) as potential evidence of their existence (e.g., Qian et al. 2012;Er et al. 2021;Esmer et al. 2022), but such claims are sometimes debated (Pulley et al. 2022).Recent radial velocity surveys may expedite the discovery of more misaligned circumbinary planets (Standing et al. 2022).For example, the BEBPOP survey contributes to the discovery of TOI-1338 c (Standing et al. 2023, , hereafter S23).TOI-1338 system is the second circumbinary multiplanet system around a main-sequence binary, following Kepler-47 (Orosz et al. 2012(Orosz et al. , 2019)).The inner transiting planet TOI-1338 b was discovered by TESS photometry, based on three transits in Sector 3, 6 and 10 (Kostov et al. 2020, hereafter K20).Later an outer planet TOI-1338 c at P ∼ 215.5 day is confirmed by radial velocity (S23).Due to the degeneracy of the RV method, the inclination of TOI-1338 c is not known, but the stability analysis in S23 constrains the mutual inclination of TOI-1338 c within ∼40 • , precluding a highly misaligned configuration.TOI-1338 system is a bright (T = 11.45 mag) CBP system in the southern sky, and TOI-1338 b could have a significant fraction of gas envelope.Therefore, TOI-1338 b is an optimal target for JWST transmission observation, which would provide us with valuable information on the migration history of the planet and, potentially, the atmospheric variability in response to variable incident fluxes.
After the discovery of TOI-1338 b, TESS revisited TOI-1338 system in the Extended Mission during 2021 and 2023, during which seven more transits of TOI-1338 b were captured.However, the observed transit mid-times offset from the predicted transit ephemerides in K20, which assumed a one-planet model, by 0.5-1 transit durations.It could be the result of the perturbation by TOI-1338 c.It's timely to incorporate existing observations to update the future transit ephemerides for the transiting TOI-1338 b.
In this work, we aim to refine the TOI-1338 system parameter.The transit timing/duration variations (TTVs and TDVs) exerted by an external perturbed can help to constrain the three-dimensional architecture of the outer orbit, as has been practiced in some systems around single stars (e.g., Dawson et al. 2014;Mills & Fabrycky 2017;Almenara et al. 2022).However, the circumbinary planet system could be more complicated, as the TTVs and TDVs are mainly sourced from the inner binary's geometric projection at the time of transit (Armstrong et al. 2013(Armstrong et al. , 2014)).Therefore we perform photodynamical modeling of the TESS light curves and published radial velocity data to self-consistently solve the four-body dynamics.The paper is structured as follows.In Section 2 we briefly introduce the new transit light curves of TOI-1338 b in the TESS Extended Mission used in this work.In Section 3 we detail our photodynamical analysis of the TOI-1338 system, and in Section 4 we present the results and future transit forecast, and conclude in Section 5.

NEW TRANSITS OF TOI-1338 b
In the following two subsections, we first describe the new TESS observation in Section 2.1.In Section 2.2 we present the light curve detrending and transit midtimes fitting and compare them with the transit predictions in K20.
The two closest sources to TOI-1338 b have angular separations of 54.1" and 58.9", both of them are fainter than TOI-1338 by 3 mags (Figure 1).None of them fall within the aperture of SPOC in the sectors when CBP transits occur.We also check the depth of the primary eclipses in the sectors of CBP transits from SPOC and GSFC light curves in S3, all of them showing consistent eclipse depth within 1.6%.The SOAR speckle imaging also indicates the absence of nearby sources within 3" (see Tokovinin et al. (2019) for a description of instrumentation).Therefore we conclude that the contamination of nearby sources should be minimal in both SPOC and GSFC light curves.
While most of the transits fall in the window of normal operation, the transit in Sector 30, 62, and 68 suffered from stray light of Earth and Moon (SPOC quality flag= 1E12).The stray light causes the background flux to rise and may cause uncertainties in the estimated background fluxes when performing aperture photometry and, thus introducing systematics in the produced light curves.As shown in Figure 2, the influence of stray light is most pronounced in Sector 30 and 68 light curves: large photometric jitters are present around the mid-transit and the post-transit portion of light curves.The transit in Sector 62 has a fall-off in fluxes near the egress phase.When binned to a 30-minute cadence, the egress/ingress of CBP transit is still well-resolved in Sector 62 and 68, but the egress is hardly resolved in Sector 30.
The variations in transit duration and midtimes are unique features for transiting CBP, and resolving the ingress/egress phase is crucial to precisely determining the magnitude of variations and measuring the movement and coordinates of different bodies.We decided to fit all the transits except for the one in Sector 30 due to its damaged egress phase.Later in Section 3.4 we will show the in-transit systematics will indeed affect the transit midtime determinations if we directly fit the light curves, but the bias is less severe if we fit the light curves globally using photodynamical modeling.

Detrending and Transit Times Fitting
We used the sap flux product to measure the midtime and duration for each transit.The detrending and fitting are devised as follows.The light curves are first fitted with quadratic limb darkening transit model in batman (Kreidberg 2015).The fitted parameters are transit midtime T mid , planet period P , semimajor axis scaled to stellar radius a/R * , impact parameter b, planetary-to-stellar radii ratio R p /R * , and two quadratic limb darkening coefficients q 1 , q 2 (Kipping 2013).We assume circular orbits during the midtime fitting.The χ 2 likelihood function is used to measure the goodness-of-fit of the model compared to the observation.The model is first optimized using Monte Carlo Markov Chain (MCMC) for 30000 iterations.Once a best-fit model is found, the light curves are trimmed to shorter segments centered at the best-fit midtime, with lengths of two times of best-fit transit durations.The trimmed light curves are detrended with a linear slope, with in-transit data given zero weight during the linear fit.The flux errors are also rescaled to make the reduced χ 2 = 1.The scaling factors for transits in each sector are also listed in Table 1.Then the model is re-optimized using MCMC for another 30000 iterations.We also calculated the estimated transit duration from our fitted parameters via Equation 141 in Winn (2010).
The measured transit midtimes and durations are listed in Table 1.Comparing the predicted transit midtimes of K20 and the observed mid-times, we found that while the observation in Sector 37 matches well with the prediction in K20, the rest of the observed transit midtimes deviates from K20 predictions by more than the magnitude of their predicted uncertainties, which is most evident in Sector 65 and 68 where deviations are readily comparable to the transit durations.The deviation between K20 predictions and TESS following observations could be caused by the imprecise system parameters, most likely the planetary orbital parame-ters, derived from the limited three transits observed in the TESS Primary Mission used in K20, and also could source from the gravitational interactions from the TOI-1338 c whose presence is unknown at that time.In Section 4.3 we will show that TOI-1338 c is indeed perturbing TOI-1338 b and contributes to the transit timing deviations we show here.Therefore we carry out photodynamical modeling to update the global parameters of TOI-1338 system in the next section.

PHOTODYNAMICAL MODELING
In this section, we present the photodynamical analysis of TOI-1338.We firstly describe our photodynamical model in Section 3.1 and the system parameter retrieval strategy in Section 3.3; we also discuss some potential bias by transit systematics in Section 3.4 and test the validity of photodynamical solutions in terms of the dynamical stability in Section 3.5.

Model Setup
The system parameters of circumbinary planet systems are often exhaustively explored by photodynamical modeling (e.g., Doyle et al. (2011b), Welsh et al. (2012)).Given initial orbital and mass parameters, the N-body integrator solves the equations of motions of binary and planets and produces in-time observables, i.e., eclipse, transit light curves, and radial velocity time series.The synthetic data are compared to the observations, and system parameter credible intervals are explored by sampling algorithms (e.g., MCMC) with related observational uncertainties.
The dynamical status of a two-planet circumbinary system can be characterized by 22 parameters, equivalent to six binary and twelve planetary osculating orbital elements (P, e, i, w, Ω and mean longitude λ), and four masses (M A , M B , M b , M c ).The ascending node angle of the binary orbit can be further set to zero, which leaves 21 parameters.Synthetic stellar eclipse and planet transit light curves can be produced by the relative positions of each body, which require the stellar and planetary radii (R A , R B , R b ), primary-to-secondary star surface brightness ratio in TESS band f TESS , and two sets of quadratic law limb darkening coefficients, for which we used the triangular sampling in Kipping (2013).In total 6 additional parameters are needed for light curve synthesis.The secondary eclipses show flatbottom in-eclipse profiles and have low SNR, therefore the limb darkening coefficients of the secondary star are not well constrained and thus are set to be the same as the primary star 2 .The dynamical modeling also issues radial velocity curves relative to the system barycenter.Thus, the system barycentric velocity γ is modeled as three instrument-dependent parameters for HARPS andESPRESSO(pre-2019 andpost-2019).In summary, 30 parameters are needed in photo-dynamical modeling in the case of the TOI-1338 system.We outline our photodynamical model as follows.The system is initialized in the Jacobian coordinate, namely, the object is initialized in the order of its semi-major axis to the primary star.We use the IAS15 integrator in the python package rebound (Rein & Liu 2012) to solve the motion of all objects.The integration is done in the center-of-mass coordinate.We use batman (Kreidberg 2015) to generate synthetic light curves from the relative coordinates of different bodies on the sky-projected plane centered on the system barycenter, after linear corrections for the light travel effect.The radial velocity data precision of ESPRESSO achieved ∼1 m/s, calling the need for relativistic correction, therefore the simulated RVs are further corrected with light travel effects, transverse Doppler, and gravitational reddening effects (Zucker & Alexander 2009;Konacki et al. 2010; Sybilski 2 A more sensible treatment for the limb darkening coefficients of the secondary star would be using the theoretical values presented by Claret (2017): u 1 = 0.15, u 2 = 0.47 for T eff = 3300 K, log g=5.0, and [Fe/H]=0.0.Using these theoretical values, the maximum difference between the new secondary eclipse light curves and our adopted model is 10 ppm, only manifesting during the ingress/egress phase.The error is much smaller than the observational error of 200 ppm and thus this barely changes the overall likelihood of the photodynamical fitting and parameter sampling.
et al. 2013).However, later we found that these will not significantly alter the fitted result as the RV variations are absorbed into the barycenter velocity.Tidal effects are not included since it suffered from uncertainty in tidal parameters and is estimated to have an amplitude ∼ 0.2 m/s, which is much lower than the total amplitude of relativistic effects (∼ 5 m/s).
For the case of the TOI-1338 system, the General Relativity (GR)-induced apsidal precession is around 9.01 arcsec per year.The apsidal precession rate induced by tides is > 10 times smaller than the GR term, following the prescription in Correia et al. (2013).We do not include the GR and tidal effects in our dynamical model.
The reference epoch is fixed at BJD=2458300.00,approximately 90 days before the first observed CBP transit in Sector 3. Using the best-fit values in K20, we integrate the TOI-1338 system into our reference epoch and record the result as an initial guess of the model parameters.The orbital parameters of the outer planet are derived from S23, with randomly generated node angles and inclinations but satisfying the system stability limit (mutual inclination smaller than 40 • , see Supplement Figure 8 of S23).
We adopt the differential evolution Monte Carlo Markov Chain (DE-MCMC, Ter Braak 2006) in emcee to sample the parameter space.The range of parameters is displayed in Table 2. Flat priors are assumed for all parameters.The χ 2 function is used to define the goodness-of-fit of the model.This is computed by the square of the differences between observation and synthetic photometric/RV data, divided by the square of observational errors.

Fitted Data
Using the photodynamical model described in the above section, we simultaneously fitted the binary eclipse light curves, transits of TOI-1338 b in TESS light curves, and the radial velocity data published in S23.
We fitted nine TOI-1338 b's primary transit events, excluding one transit event in Sector 30 since it is severely affected by systematics.To speed up the integration speed and prevent the potential bias delivered by the instrumental systematics or out-of-eclipse trend (Orosz et al. 2019), we select nine primary and secondary eclipses, respectively, for the following photodynamical modeling, based on the ranking of their reduced χ 2 .Figure 12 and 13 in Appendix B show the selected nine primary eclipses and secondary eclipses.The transit and eclipse data are detrended according to Section 2.2.Once an initial dynamical result that can reproduce the observed light curves is obtained, the light curves are detrended again with the more precise transit/eclipse midtimes and durations inferred from the dynamical result.
61 HARPS and 123 ESPRESSO RV data points are published by S23.Many of these RV data have been flagged to be potential FWHM/bisector outliers, or due to the data being taken at the time of binary or planetary transits.We exclude RV data with outlier flags for our photodynamical model, leaving 58 HARPS and 103 ESPRESSO RV data to be modeled.

Planetary System Parameters Sampling
We first perform a joint fit to the stellar eclipse, transit photometry, and HARPS/ESPRESSO radial velocity data.At the starting stage, 100 chains evolved from the initial guess with small offsets to explore the parameter space.After the chains remained stable in small ranges, we evolved 100 chains for another 690000 steps.We remove the first 350000 steps as the burn-in phase.The rest of the 340000 steps from 100 chains are concatenated together, and they are sampled every 1000 steps to make the posterior sample with a total size of 34000.The Gelman-Rubin statistics R of all parameters are smaller than 1.05 (Brooks & Gelman 1998).The autocorrelation steps of fitted parameters vary between 1972 and 5213, giving at least 6522 effectively independent samples in total.The median and uncertainties (defined as the 16th and 84th percentiles of the posterior distribution) of each fitted parameter are listed in Table 2.
We examined the residuals of all eclipse, transit, and RV residuals and found the best-fit solution from joint fit provides a good match to the light curve data.However, we found a very significant periodicity correspond-ing to the binary period in HARPS RV residuals (FAP > 0.1%) but none in ESPRESSO residuals.In other words, there is inconsistency in the binary RV motion between HARPS and ESPRESSO RV data.The peakto-peak amplitude of HARPS residual is ∼ 14 m/s, and peaks around when the binary is in its conjunction.This is not likely a problem associated with our photodynamical code as this inconsistency persists when we perform a joint RV-only Keplerian fit analysis (see details in Appendix A).
To circumvent the potential bias on the binary orbital parameters brought by RV inconsistency, we experimented with two additional photodynamical fits performed individually to ESPRESSO and HARPS data (and hereafter we refer to these two experiments as esp-only and harps-only analyses), with the light curve dataset identical to the joint fit.We took 100 samples from the joint fit posterior samples and evolved them for 300000 steps.We trimmed the first 150000 steps as the burn-in phase, and the Gelman-Rubin statistics R for the remaining chains are smaller than 1.04 for all parameters.The effective sample sizes are at least 4263 and 3624 for esp-only and harps-only posteriors, respectively.We kept the last 150000 steps and sampled every 1000 steps to make the posterior distributions.The median and uncertainties of fitted parameters from the esp-only and harps-only analysis are presented in the second and third columns of Table 2. Most parameters from esp-only and harps-only fit are 1-2 σ deviant from joint analysis, while some of them show 3σ discrepancy.For the binary orbital parameters, the joint fit and esp-only fit show close agreement, while the binary eccentricity from harps-only fit is around 3σ larger than those from the other two analyses, which probably stems from the RV inconsistency from HARPS data.The binary mass ratio and orbital parameters of TOI-1338 b inferred from the three analyses are different too, with joint and esp-only analysis having 3σ discrepancy on binary mass ratio and that from harps-only being intermediate between the former two.We provide some intuitions of the potential causes.TOI-1338 is a single-line binary, so only the mass function (related to the mass ratio, absolute mass, and binary period, eccentricity, inclination) can be constrained from RV data.The degeneracy between the mass ratio and the mass function is solved by accounting for the exact timing of planetary transit.For example, the peak-to-peak amplitude of transit timing variations nonlinearity of CBP transit is, in the first order, related to the binary phase and the exact projected stellar location, which is further related to the semi-major axis of the transited star and hence the mass ratio can be constrained (Schwamb et al. 2013;Armstrong et al. 2013;Kostov et al. 2013Kostov et al. , 2014)).The CBP transit durations also depend on the relative velocity difference between the binary star and planet.For the latter one, the velocity of the planet is also related to the total masses of the inner binary, and thus can also contribute to solving the mass ratio from the mass function.This degeneracy is also supported by the clear correlation between the mass ratio and the inner planetary orbit parameters (P b , λ b , eccentricity vectors).Following this thread of thought, we attributed the discrepancy of mass ratio derived from the two analyses to the difference in planetary orbit solutions, since we found the derived mass function from joint and esp-only analysis is consistent within 1σ.
Despite the difference in inner planetary orbit solutions, the best-fit χ 2 of the transit light curve between joint and esp-only analysis has a negligible difference, which is χ 2 transit = 10166.8and χ 2 transit = 10164.6for 10183 cadences, respectively.But given the inconsistency between HARPS and ESPRESSO data, we prefer esp-only solution over the joint-fit one.Between esponly and harps-only solutions, the χ 2 in transit is almost identical, and the harps-only solution is only marginally preferred over the esp-only ones in terms of the χ 2 eclipse .However, we think the outer planet's orbit is better constrained by the ESPRESSO data, for its higher precision and more observations.At the writing of the manuscript, Sebastian et al. (2024) published the absolute dynamical masses of TOI-1338 AB for 1.098 ± 0.017 M ⊙ and 0.307 ± 0.003 M ⊙ , respectively.This measurement is more aligned with our esp-only solution.Additionally, in Section 3.5 we will show that the esp-only solution are more dynamically stable compared to the other two sets of solutions.Therefore, among the three analyses presented in this section, we adopt the solution of esp-only analysis.
In our nominal photodynamical model, we assume constant zero flux dilution for all TESS sectors.For completeness of the photodynamical discussion, we also run a model accounting for different dilution levels in different sectors.Applying a unique dilution factor to each sector is necessary for missions like TESS since the star's location in the CCD plane in each sector is different and so is the amount of flux from the nearby stars expected to contain within the target's aperture.In this simulation, we include the primary and secondary eclipse light curves in the same nine sectors as the transits were observed to help better constrain the dilution factor.The dilution factors are allowed for negative values to reflect the situation where the background fluxes are potentially over-substracted.However, we found the dilution factor of each sector is consistent with zero within 1σ, and the important parameters like stellar mass ratios did not change significantly.We attribute this to nearly consistent or negligible dilution within TOI-1338 in different sectors and proceed with the remaining discussion with the model without dilution consideration.

The Effect of In-Transit Noise
In Figure 3, we compare the posteriors of transit midtimes of all sectors from the direct fitting (Table 1) and photodynamical fitting.Generally, the uncertainties of transit midtime from the global photodynamical model are smaller than those from direct fitting, and the midtime posteriors from two analyses are consistent within 1σ, except for the transits in S62, S65, and S68.The in-transit systematics that cannot be removed from the detrending technique described above could potentially bias the midtime estimates and thus the photodynamical modeling.To investigate whether the in-transit systematics, especially in S62 and S68, have a significant effect on the photodynamical modeling, we took an independent analysis of midtime fitting incorporating the Gaussian process (GP) to see if the systematic-corrected midtimes are consistent with the results from the photodynamical model.The GP could model the noise and transit altogether to yield more unbiased midtime results compared to the direct fitting.This is done by using the python package juliet (Espinoza et al. 2019).
We show the GP-correct midtime posteriors in S62, S65 and S68 in Figure 3.For the S62 and S68, the midtimes derived from the photodynamical modeling are more discrepant with those from direct fitting but are more consistent with the results of GP-corrected midtimes.This indicates that our midtime posteriors from the global photodynamical modeling do not greatly suffer from the systematics in S62 and S68 as the direct fitting does.However, the transit in Sector 65 displays ∼ 2σ difference between photodynamical modeling and other results.We postulate that this is caused by the very sharp flux jump near the egress phase of this transit event (see Figure 2), which could not be sufficiently modeled by the GP.This could potentially make the predicted egress time from the direct fitting end earlier and thus the midtime will be a little earlier too.

Stability Analysis
To investigate the stability of the posterior solutions from our photodynamical modeling, we calculated a dynamical Mean Exponential Growth factor of Nearby Orbits (MEGNO; Goździewski et al. 2001;Cincotta et al. 2003) for all three sets of solutions in Table 2. MEGNO  is a numerical tool to differentiate chaotic and stable motion.Specifically, if the orbits are stable, the final ⟨Y ⟩ will converge to 2, else it will increase as time.
We integrate each solution with the whfast (Rein & Tamayo 2015) integrator for 50 kyr and calculate the MEGNO value ⟨Y ⟩ of the system.The cumulative density function of calculated ⟨Y ⟩ for three sets of solutions are shown in Figure 4.The median ⟨Y ⟩ value of joint-, esp-, and harps-only solutions are 6.13, 1.98, and 7.84, respectively.This result suggests the majority of posterior solutions (∼ 80% for 1 < ⟨Y ⟩ < 3) from esponly analysis is stable within 50 kyr, while the other two sets of posteriors are more likely to become unstable within the integrated time.Examining the unstable solutions with ⟨Y ⟩ > 3 shows that the larger the eccentricity of e c TOI-1338 c is, the more likely the system would get unstable.Note that the joint and harps-only solutions yielded ∼ 0.1 for e c while esp-only solutions have ∼ 0.047.This is also consistent with the conclusion of S23 that the eccentricities of both planets should be lower than 0.1 to keep the system stable.Based on the results from the stability analysis, it is advised to take the esp-only solutions as the final adopted solutions from our photodynamical modeling.

RESULTS
Figure 4. Dynamical stability (in terms of the displacement of MEGNO ⟨Y ⟩ factor from 2) for three sets of photodynamical solutions from Table 2.Each solutions are integrated for 50 kyr.The closer the ⟨Y ⟩ is to 2, the more stable the system is in the integrated time.The stability analysis shows that the majority of esp-only solutions is stable while the other sets of solutions are more prone to be unstable within the integrated time, mainly owing to the different eccentricity of TOI-1338 c from the different solutions (see Section 3.5).
In our photo-dynamical modeling, we analyze the mass constraint of TOI-1338 b, which could potentially have an extremely low density as 0.137 ± 0.026 g/cm −3 , in Section 4.2.We also discuss the mutual inclination of non-transiting planet TOI-1338 c relative to the binary plane in Section 4.3.Then we present the future primary transits forecast in Section 4.4.

Eclipse Timing Variation Analysis
The mass of the circumbinary planet is traditionally constrained by the eclipse timing variations (e.g., Doyle et al. 2011a;Welsh et al. 2012).The planet will gravitationally perturb the binary and cause the apsidal precession of binary orbits, which would appear as a divergence of the binary period derived from primary and secondary eclipses if the binary orbit is eccentric.Futhermore, if the planet is non-coplanar with the binary, the inclination of binary orbit would also go through precession and manifest as the eclipse depth variations, though this would work most effectively for grazing eclipses (see, e.g., Socia et al. 2020;Goldberg et al. 2023).The planet's mass could also be constrained by the shortterm "choppings" in the Observation minus Computed (O-C) diagram, the time-scale of which is comparable to planetary orbit (Rappaport et al. 2013).
In the best-fit model excluding GR, the total apsidal precession of binary orbit is 31.35arcsec per year.The precession rate caused by planet b is 14.56 arcsec per year, slightly larger than the expected GR-induced precession rate (∼ 9.01 arcsec per year) and is around 10 times higher than the tidal-induced precession rate, but is smaller than the precession rate induced by the outer planet (∼ 17.91 arcsec per year).In the 5-year TESS's observation baseline, this precession rate will lead to a divergence of 0.8 minutes between primary and secondary eclipses when they are fitted to common ephemerides.However, this divergence amplitude is not readily observable because it is much smaller than the typical uncertainty of the secondary eclipse midtimes (∼5-6 min), as seen from Figure 5.The short time-scale oscillations of O-C curves have an amplitude of 0.03 min, dominated by the outer-planet perturbation.We also tried out a few cases when TOI-1338 c is misaligned by mutual inclination up to 40 • , and the amplitude of primary ETVs is still smaller than the typical uncertainties of 0.3 min.Therefore the ETVs are not likely to constrain the mass of TOI-1338 b nor the inclination of TOI-1338 c.

The Mass of TOI-1338 b
Our adopted photodynamical solutions in Table 2 show that the mass of TOI-1338 b is 11.3 ± 2.1 M ⊕ .To determine the source of the inner planetary mass constraint, we run two additional groups of modeling with modified fitted data to see whether the mass of TOI-1338 b is constrained or not.The first group of modeling fits the stellar eclipse photometry in the first year of TESS observation (S3-S10) and the full set of RV observations, this choice aligns with K20 and removes the possible mass constraints delivered by long-term ETVs.However, the mass constraints for planet b remain the same as the adopted solutions.The second group of modeling fits all TESS photometry but without RV data, with the mass of the binary and outer planet fixed at the best-fit value of the adopted solution.In this case, the mass distribution of the TOI-1338 b is highly skewed to zero, with a 95% upper limit of 26.1 M ⊕ .Therefore we conclude that the mass reported in Table 2 are constrained from the RV data.
We then look for the RV signal of TOI-1338 b in the current RV dataset.This is done by subtracting the observed RVs from the secondary star and outer planet's RV signal simulated issued from the integrator, the parameters of which are taken from the best-fit solutions.The residuals should correspond to the RV signal of TOI-1338 b.The phase-folded residuals are shown in Figure 6.In the left and middle panel of Figure 6, corresponding to the best-fit TOI-1338 b's RV signal from the joint and esp-only analyses, we observe clear modulation in observed RV residuals folded at the period of TOI-1338 b, with a semi-amplitude of ∼ 1.0 and 1.5 m/s, respectively.To assess the compatibility of this modulation with our model, we have included in Fig- It is important to note that some deviations of the observed RV signals from the model predictions are apparent, particularly in certain phase bins.These deviations may be attributed to the uneven coverage of RV observations over the planet phases, especially in esp-only results where the current RV data are sparser in the latter half of the planetary period (phase ∼ 70 − 94 day).
The harps-only solution provides a mass constraint of TOI-1338 b consistent with zero.As is also seen in the right panel of Figure 6, there is no apparent signal when the RVs from the secondary star and TOI-1338 c are subtracted.After subtracting the best-fit signal of TOI-1338 b, the χ 2 changed from 57 to 58, so there is no improvement made.It is likely that the RV signal of TOI-1338 b is completely submerged in the noise of HARPS RV measurements.Also note that in Section 3.3 we report inconsistency in binary RV signal between HARPS and ESPRESSO with amplitude ∼ 14 m/s.Such large systematics could be detrimental to the search of planets with signals as weak as only ∼1 m/s.
The study of S23 and K20 have both formerly provided mass constraints on TOI-1338 b.S23 uses RV datasets identical to ours and perform Keplerian fit to the RV and find tentative signals at periodicity around 100 days consistent with TOI-1338 b, though the significance of this signal is not high enough and thus they report a 99% upper limit of 21.8 ± 0.9 M ⊕ for TOI-1338 b.Note that their analysis is purely based on RV data, independent of any prior knowledge of orbital parameters TOI-1338 b, including the orbital phases which could be precisely informed by the transit events.K20 derived 33 ± 20 M ⊕ for the planet by photodynamical modeling of 1.5 yr TESS photometry and CORALIE/HARPS RV measurements.The mass constraint mainly comes from the seven HARPS RV measurements obtained at that time.Our method is more similar to the analysis of K20, except that we incorporate more RV measurements from ESPRESSO with higher precision and yield more precise mass constraints.

The Inclination of TOI-1338 c
In Figure 7, we show evidence that the TOI-1338 b is gravitationally perturbed by TOI-1338 c, which manifests as TTVs and TDVs in models where we fixed the mass of TOI-1338 c as zero.The variations in transit timing are more prominent, in several transit events the deviations can be detected with several σ, whereas  O-C) diagram of the midtimes of primary and secondary eclipses, fitted to the common period.The periods fitted from primary and secondary eclipses are 14.6085574± 9×10 −7 day and 14.6085774± 1.58×10 −5 day, respectively.Neither the divergence of primary and secondary ephemerides (0.8 min) nor the dynamical delay of the primary eclipse due to the short-term perturbations (0.03 min) predicted from the best-fit model is observable given the large uncertainties in midtimes of primary and secondary eclipses (0.36 min and 5 min, respectively).the variations in transit durations are close to the measured uncertainties.Moreover, the inclination of the outer planet's orbit also changes the TTVs, which may allow us to put meaningful three-dimensional constraints on the outer planet's orbit, which is inaccessible to the RV method by which TOI-1338 c is discovered.
The mutual inclination between the binary and outer TOI-1338 c can be calculated via The three analyses yielded a mutual inclination of TOI-1338 c around 10 • relative to the binary orbital plane, consistent with a coplanar configuration.However, our adopted model prefers a non-transit nature of TOI-1338 c, even though we did not impose any criterion on it.
To explore how sensitively the data constrain the mutual inclination of TOI-1338 c ∆I c , we run a suite of additional fits, with i c being between 50  , the optimization will need to run longer chains).Optimization ends when we visually confirm that the χ 2 remains constant for more than 50% of the steps.The optimized χ 2 map of the dynamical fitting is shown at the bottom of Figure 8.We find the minimum χ 2 exists in (i c , Ω c )=(95 • , 10 • ), consistent with the solutions provided in Table 2.That is, a lower mutual inclination of TOI-1338 c is favored by the current data, and the fits become worse in larger mutual inclinations, as shown in the top panel of Figure 8.When examining the χ 2 of transit light curves and RV data separately, we found that when TOI-1338 c is fixed at high ∆I c the optimized results could provide a good match to the observed transit midtimes and durations, while the fits of RV data become worse.
However, compared to the nominal photodynamical models in Section 3, only the eclipse photometry is removed in the models in this section, therefore the degree of freedom of these model optimizations is still as high as 10336.At higher ∆I c , the ∆χ 2 is only ∼35-40, which may not be a statistically significant argument against high ∆I c if seen from the perspective of reducedχ 2 .Therefore, based on the current data, lower mutual inclination between binary and outer planet's orbit is marginally favored over high mutual inclination scenarios.

Future Transits Forecast
We present the future primary transit of TOI-1338 b in Table 3.It shows that TOI-1338 b will permanently transit the primary star, which is similar to Kepler-47 b (Orosz et al. 2019).This is partly due to the coplanarity of the orbit of TOI-1338 b relative to the binary planet (∆I b ∼ 0.12 • ) and the sufficiently large primary radius.Integrating the best-fit model with rebound, in Figure 9 we present the impact parameter of TOI-1338 b in each conjunction with the primary and the secondary Top: the evolution of impact parameters of primary (red) and secondary transit (blue) integrated from the best-fit model.It shows that TOI-1338 b will permanently transit the primary star but hardly transit the secondary star.Bottom: the skyprojection view of the orbits of binary and TOI-1338 b.The orbit of TOI-1338 b is integrated 7000 days to show the effect of nodal precession due to interactions with binary stars.The size of the point scales with the line-of-sight distance, thus smaller points are behind the binary orbits.The red and blue regions show the crossing area of stellar disks of primary and secondary stars, accounting for the binary orbital motions.Transits can occur when the planet orbit and the stellar crossing area are intersected.The aspect ratio is 70:1.Due to the small mutual inclination of TOI-1338 b, the variation range of planetary inclination is very small, thus its orbit will always intersect the primary star's orbital crossing region.The secondary star has a smaller radius, thus the stellar and planetary orbits are not mutually intersected, precluding potential transitability.
star.Due to perturbation of the inner binary, the orbital plane of TOI-1338 b will undergo precession and the inclination of TOI-1338 b will follow sinusoidal-like variations.
Our best-fit models suggest that no secondary transit occurs in the current TESS observation.Following Martin & Triaud (2015), the minimum mutual inclination of TOI-1338 b ever transiting secondary star is 0.277 • , whereas our posterior samples yielded a 0.127 +0.025• −0.024 .This partly explains why secondary transits of TOI-1338 b are rare in our posterior models.The orbit of TOI-1338 b will hardly precess into an orientation that the planetary and secondary star's orbits are mutually intersected.
With a TESS mag of 11.45 mag and a period of 95 days, TOI-1338 b is an ideal circumbinary planet target for JWST transmission observation.TESS is scheduled to re-observed TOI-1338 in Year 7 observation, and another three primary transits of TOI-1338 b are expected to occur in the observing window of Sector 89, 93, and 96 between 2025 Match and September.We encourage more follow-up observations to be carried out to TOI-1338 b transits to update and refine future transit ephemerides.

DICUSSION AND CONCLUSION
In this work, we perform a system update of the TOI-1338 with two circumbinary planets.Combined with the latest TESS light curves and published radial velocity data, we carry out photodynamical analysis to selfconsistently account for the mutual interactions between bodies within the systems.We determined a coplanar configuration for TOI-1338 c relative to the binary plane ∆I c = 9.1 +6.0• −4.8 , with a 99% upper limit of 22 • .Therefore the genuine mass of TOI-1338 c is determined to be 75.4+4.0 −3.6 M ⊕ .This constraint mainly comes from the additional number of observed transits of TOI-1338 b and the precise RV measurements, where the gravitational pull from the outer planet starts to manifest in the transit timing variations of the inner planet (see Figure 8).Together with the spin-orbit alignment between the primary star and the orbit of TOI-1338 b (projected spin-orbit angle of 2. • 8 ± 17. • 1) (Hodzic et al. 2020), the whole system is in an alignment between primary's spin axis and binaryplanet orbit configuration.This suggests the whole system forms from a single disk with no breaking or warps.
With our refined mass (11.3 ± 2.1 M ⊕ ) and radius (7.66 ± 0.053 R ⊕ ) for TOI-1338 b, we derived a low bulk density of 0.137 ± 0.026 g/cm −3 .This is comparable to circumbinary planet Kepler-47 c which also has a low bulk density of 0.17 g/cm −3 .In Figure 10 we show TOI-1338 b and other low-mass sub-Saturns in the massradius diagram, with the tide-free thermal evolution models for sub-Saturns from Millholland et al. (2020) at TOI-1338 b's age (∼4 Gyr, K20) and insolation.We estimate the envelope mass fraction, f env , to exceed 50%, implying that TOI-1338 b possesses a substantial envelope, akin to several super-puffs found around single stars (e.g., Kepler-79 d with f env = 49.1 ± 1.9%, Millholland et al. (2020); K2-24 c with f env = 52 +5 −3 %, Petigura et al. ( 2018)).Lee & Chiang (2016) suggests planets with such low density might have migrated from a distant region where the disk is cool and has low opacity, enabling efficient cooling of the accretion and forming a massive envelope even with a low-mass core.This hypothesis aligns with the general picture that most observed CBPs may have undergone migration during their formation process (e.g., Pierens & Nelson 2013;Coleman et al. 2023a,b), due to the local low growth rate of planetary cores near the stability limit, which hinders the in situ formation of CBPs (Pierens et al. 2020(Pierens et al. , 2021)).
With a TESS-band magnitude of 11.45 and a Jband magnitude of 10.95, TOI-1338 is an optimal target for future follow-up and transmission observation.Taking an average equilibrium temperature of 495 K for TOI-1338 b, we calculated the transmission spectroscopy metric of TOI-1338 b to be 86±17, which is slightly below the recommended threshold of 96 for sub-Jovian planets (Kempton et al. 2018).Although previous measurements show featureless transmission spectra for low-density sub-Saturns (Libby-Roberts et al. 2020;Chachan et al. 2020), hypothetically owing to the high altitude hazes or circumstellar rings (Wang & Dai 2019;Gao & Zhang 2020;Ohno & Fortney 2022), it would still be crucial to test this scenario in the infrared band to see if the hazes or rings are genuinely present and inflate the apparent radius of the planet.
We provided an updated future primary transit ephemerides of TOI-1338 b in Table 3.In our posterior models, TOI-1338 b will permanently transit the primary star at all nodal precession phases, owing to the extreme coplanarity with the binary plane and the large radius of the primary star.We encourage more follow-up observations of TOI-1338 b's primary transit to refine the system parameters.This paper includes data collected with the TESS mission, obtained from the MAST data archive at the Space Telescope Science Institute (STScI).Funding for the TESS mission is provided by the NASA Explorer Program.STScI is operated by the Association of Uni-versities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.The TESS data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute.The TESS light curves used in this work can be accessed via MAST (Team 2021).The GSFC TESS light curves can be accessed in Powell, Brian (2022).

Figure 1 .
Figure 1.An example of the Target Pixel File of TOI-1338 (the zeroth target) in Sector 6 and the pipeline aperture (orange squares) from SPOC.The three closest nearby sources to TOI-1338 have angular separations of 54.1", 58.9", and 60.7", respectively, all of them have ∆G > 2.7 mag.This figure is produced by tpfplotter (Aller et al. 2020).

Figure 2 .
Figure 2. The ten primary transits of TOI-1338 b in TESS observations.We binned the data into the 30-minute cadence for clarity of display.Red lines are the synthetic light curves from the best-fit photodynamical model of TOI-1338 b (see Section 3).The residuals between the best-fit model and observation are shown below.Cadences affected by stray light flagged by SPOC are marked as blue.The photodynamical modeling used all primary transits except for the one in Sector 30.

Figure 3 .
Figure 3.The comparisons between the fitted nine transit midtimes posteriors derived from direct fitting (blue, Section 2.2), photodynamical model (red, Section 3.3), and Gaussian process correction (orange, see text in Section 3.4).
ure 6 the predicted RV signals of TOI-1338 b from the best-fit model at each of the RV observation epochs, represented as red points.The figure demonstrates a close match between the predicted signals and the observed modulations.

Figure 5 .
Figure5.The observation minus the calculation (O-C) diagram of the midtimes of primary and secondary eclipses, fitted to the common period.The periods fitted from primary and secondary eclipses are 14.6085574± 9×10 −7 day and 14.6085774± 1.58×10 −5 day, respectively.Neither the divergence of primary and secondary ephemerides (0.8 min) nor the dynamical delay of the primary eclipse due to the short-term perturbations (0.03 min) predicted from the best-fit model is observable given the large uncertainties in midtimes of primary and secondary eclipses (0.36 min and 5 min, respectively).

Figure 6 .
Figure 6.Left to right: the best-fit RV signals of TOI-1338 b from joint, esp-only, and harps-only analyses.The RV data from HARPS, and ESPRESSO prior to and after 2019 are shown in blue, orange, and green, respectively.Two periods of TOI-1338 b are shown in the phase-folded figure.Black symbols with error bars are RVs binned in 11 evenly-spaced intervals within one period.The simulated RV signals of TOI-1338 b in each analysis at the observation epoch are also plotted in red dots.
• and 130 • and Ω c ranging between -40 • and 40 • and a step of 5 • .All other parameters are fixed at the best-fit value of joint analyses in Table 2, except that P b , √ e b cos ω b , √ e b sin ω b , λ b , √ e c cos ω c , √ e c sin ω c , λ c , and P c are further optimized by DE-MCMC.The inclination and node angle of the inner planet are not sensitive to the transit timing and are therefore kept fixed.These parameters are optimized against the observational dataset identical to the adopted solutions, except that binary eclipse light

Figure 7 .
Figure7.The difference between the measured transit midtimes and durations of TOI-1338 and those derived from best-fit solution (black points).The error bars are the measured uncertainties from the observation (listed in Table1).We also show how the transit midtimes/durations of TOI-1338 b would change if we tweaked the TOI-1338 c's planetary mass to be zero (blue) and inclination to be 110 • (green), the deviations shown by the tweaked models indicate TOI-1338 b is gravitationally interacting with TOI-1338 c.

Figure 8 .
Figure 8. Top: optimized ∆χ 2 map of transit light curve of TOI-1338 b and RV data fitting as a function of the mutual inclination between TOI-1338 c and binary orbit ∆Ic.The red vertical dashed line denotes the stability limit (∼ 40 • ) derived from S23, beyond which the system would be unstable.Bottom: The isolines of optimized ∆χ 2 map of transit light curve and RV data fitting assuming different orbital inclination ic and node angles Ωc of TOI-1338 c.The isolines of mutual inclinations of 5 • , 10 • , 20 • , and 40 • are also plotted in gray dashed lines.

Figure 9 .
Figure9.The evolution of orbital plane and the change of impact parameter of TOI-1338 b's primary and secondary transits.Top: the evolution of impact parameters of primary (red) and secondary transit (blue) integrated from the best-fit model.It shows that TOI-1338 b will permanently transit the primary star but hardly transit the secondary star.Bottom: the skyprojection view of the orbits of binary and TOI-1338 b.The orbit of TOI-1338 b is integrated 7000 days to show the effect of nodal precession due to interactions with binary stars.The size of the point scales with the line-of-sight distance, thus smaller points are behind the binary orbits.The red and blue regions show the crossing area of stellar disks of primary and secondary stars, accounting for the binary orbital motions.Transits can occur when the planet orbit and the stellar crossing area are intersected.The aspect ratio is 70:1.Due to the small mutual inclination of TOI-1338 b, the variation range of planetary inclination is very small, thus its orbit will always intersect the primary star's orbital crossing region.The secondary star has a smaller radius, thus the stellar and planetary orbits are not mutually intersected, precluding potential transitability.

Figure 10 .
Figure10.The masses and radii of low-mass sub-Jovian planets (4R⊕ < Rp < 9R⊕, 2M⊕ < Mp < 40M⊕) on the NASA Exoplanet Archive as of 2024 April 7 with masses and radius measured better than 30%, color-coded by the insolation fluxes.TOI-1338 b is highlighted in magenta.We also plot the transiting circumbinary planets in the sub-Saturn regime with larger white-filled circles.Dashed lines show the expected mass and radius for sub-Saturn planets at 4 Gyr old and 10 Earth insolation S⊕ from the tide-free, thermal evolution models ofMillholland et al. (2020), for H/He envelope mass fraction of 30%(blue), 45%(green), and 48%(orange).

Figure 11 .
Figure11.Left: the periodogram of HARPS, ESPRESSO, and all RV observations after subtracting the binary and TOI-1338 c MAP signal from Keplerian orbit fit, the dotted, dot-dashed, and dashed line represent false alarm probability of 0.1%, 1%, and 10% .Right: the phase-folded plot of HARPS ESPRESSO, and all RV residuals at binary period.We show two periods of folded residuals to display the modulation.The binary's Keplerian RV signal is also plotted in red lines to compare with the residual phases.

Figure 12 .
Figure 12.The nine primary eclipse light curves used in the photodynamical modeling.Grey errorbars are TESS 2-min light curves, and red lines are synthetic light curves from the best-fit models.Residuals are plotted below each panel.

Figure 13 .
Figure 13.The nine secondary eclipse light curves used in the photodynamical modeling.Labels are the same as Figure 12

Table 1 .
Comparison Between the Predicted and Observed Transit Midtimes and Durations of TOI-1338 b