Photodynamical Modeling of the Compact, Multiply Eclipsing Systems KIC 5255552, KIC 7668648, KIC 10319590, and EPIC 220204960

: We present photodynamical models of four eclipsing binary systems that show strong evidence of being members of higher-order multiple systems via their strong eclipse timing variations and/or via the presence of extra eclipse events. Three of these systems are from the main Kepler mission, and the other is from the K2 mission. We provide some ground-based radial velocities measurements for the three Kepler systems and make use of recent light curves from the TESS mission. Our sample consists of two 2 + 1 systems and two 2 + 2 systems. The first 2 + 1 system, KIC 7668648, consists of an eclipsing binary ( P bin = 27.8 days) with late-type stars ( M 1 = 0.8403 ± 0.0090 M ⊙ , R 1 = 1.0066 ± 0.0036 R ⊙ and M 2 = 0.8000 ± 0.0085 M ⊙ , R 2 = 0.8779 ± 0.0032 R ⊙ ) with a low-mass star ( M 3 = 0.2750 ± 0.0029 M ⊙ , R 3 = 0.2874 ± 0.0010 R ⊙ ) on a roughly coplanar outer


Introduction
Roughly half of nearby Sun-like stars are in binary systems (e.g., [1]).Furthermore, perhaps 10% of those binary systems are themselves members of triple or quadruple systems (e.g., see [2]).Large space-based photometric surveys such as Kepler [3][4][5] and TESS [6] have discovered thousands of eclipsing binaries (e.g., [7][8][9]).Owing to the often long, continuous photometric coverage with high signal-to-noise, many of these eclipsing binaries have been shown to have nearby companions found from eclipse timing variations (e.g., [10][11][12]) or from having multiple eclipse events with distinct periodicities (e.g., [13]).It is well known that binary stars are the foundation upon which stellar astrophysics is built because it is through eclipsing binaries, we can measure stellar masses and radii with high precision.Triple stars offer even more, especially the low-mass, short-period systems with extra eclipses, none of which were known prior to the Kepler mission.Very succinctly, triple systems can often provide more information, and this allows a better determination of the system properties.This additional precision allows us to probe deeper physics than could be done with simple binary stars.Furthermore, the sample of triple-star architectures can be used to statistically test the predictions of tidal theory and theories of the formation of multiple-star systems.
The addition of a third body companion with a non-zero mass certainly complicates the orbital motions of the binary star components: generally, the orbits are no longer purely Keplerian.But this complexity provides additional information that we can exploit.In particular, the orbits are quite sensitive to the exact positions and momenta of the three bodies and the observed light curve is very unforgiving: if the orbits, masses, and radii are not known very precisely (typically with relatives of less than a few percent) the model light curve will not match the observations.This sensitivity is a curse in that it can be exceedingly difficult to match the data (e.g., the title of a 2014 paper by Marsh et al. [14] says this nicely: "KIC 2856960: the impossible triple star").However, when the correct model is found, it is a blessing since the stringent observational constraints impose very tight limits on system parameters and would perhaps allow us to measure the effects to which we were previously insensitive.But from where exactly does the extra information in a triple system come?
In an isolated eclipsing binary system, there are two locations per orbit that the relative positions and velocities of the stars are directly measured with photometry, namely at the conjunction phases when the eclipses occur (the velocity constraint comes from the duration of the eclipses).But in an eclipsing triple system, when additional eclipses/transits occur they tell us the locations and motions of the bodies at those tertiary eclipse times.Because the orbital periods of the binary and the outer orbit are generally not related, tertiary eclipses can occur at any binary phase.Consequently, an eclipsing 3-body system offers qualitatively new information, stroboscopically probing more of the binary orbit.It is this additional information that allows the extremely precise mass and radii determinations.
In this work we present comprehensive photodynamical models of KIC 5255552 and KIC 7668648, both of which show extra eclipse events.KIC 5255552 turns out to be a 2 + 2 system (e.g., a pair of close binaries that orbit the system barycenter), and our model here is the first one that can account for all the extra eclipse events and the large eclipse timing variations.KIC 7668648 is a compact 2 + 1 system (e.g., a binary with a circumbinary companion that orbits both components of the binary).Owing to the numerous extra eclipse events, along with the inclusion of ground-based radial velocity measurements and additional data from TESS, we can determine the masses and radii of all three stars with ≈1% relative uncertainties.KIC 10319590 is a 2 + 1 system that shows dramatic changes in the eclipse depths over time (the eclipses disappeared completely after the first ≈400 days of the Kepler mission), but without any observed extra eclipse events above the noise level of the data.EPIC 220204960 is a 2 + 2 system that shows strong evidence of tidal interactions between the two binaries, again without any observed extra eclipse events above the noise level of the data.
The layout of this article is as follows.In Section 2, we summarize the available Kepler and TESS photometry and discuss ground-based spectroscopic observations of the three Kepler systems.In Section 3 we discuss our photodynamical model and code that is used to analyze the observational data.In Section 4 we discuss the actual modeling of the four systems of interest.In Section 5 we briefly discuss the stellar parameters and age constraints from model isochrones, some basic long-term dynamical properties, and comparisons with previous work.Finally, Section 6 contains a brief summary.

Observational Data
Table 1 summarizes the observational data available for our four targets.We give below general summaries of the sources for our observational data.Additional details for each source are given in Section 4.

Kepler
NASA's Kepler Mission was launched 7 March 2009 with the goal of discovering and characterizing Earth-sized planets in the habitable zone of their host stars [3][4][5].Science operations started 13 May 2009 (Quarter 1, BJD 2,454,964.51)and continued through the end of nominal spacecraft operation on May 11, 2013 (Quarter 17, BJD 2,456,424.00).Nearly 200,000 targets were monitored with high precision at a ≈30 minute cadence (hereafter long cadence) with a duty cycle close to 92% in some cases.A small subset of the targets was observed with a cadence of 1 min (hereafter short cadence).For all the targets discussed here we downloaded the latest version of the processed Kepler data in binary FITS format (Data Release Version 25) from the Mikulski Archive for Space Telescopes (MAST 1 ).
The Kepler light curves sometimes have individual cadences with bad or corrupted data.Each cadence has a data quality flag, where a flag value of zero indicates an observation with no known problem, and a flag value of >0 indicates potential problems owing to certain spacecraft events (such as a reaction wheel desaturation event) or other events (such as cosmic rays near an aperture).The MAST Kepler Archive Manual 2 gives the definition of the qualified flags.In our experience, individual cadences with quality flags with values > 16 are generally unreliable and it has been our practice to remove those cadences.After the removal of individual cadences with quality flags > 16, the Kepler light curves generally need detrending to remove instrumental drifts and occasionally modulations due to star spots.We used the iterative process described in Orosz et al. [15] arrive at the final normalized light curves.Briefly, the light curves are initially detrended "by hand" using an interactive code to mask out eclipses and fit cubic splines to the out-of-eclipse parts of the light curve, where the number of nodes for the splines is typically between 10 and 50.Once a good solution is found, the best-fitting model is used to precisely determine the various times of ingress and egress and, hence, the duration of each event.We keep only portions of the light curve centered on eclipse and transit events, where, as a rule of thumb, the length of a "chunk" centered on an eclipse event is twice the duration of the eclipse when possible.A fifth-order polynomial is fit to each chunk where data that occur between ingress and egress are given zero weight (other orders for the polynomials were tried, but generally, the fifth-order polynomials gave the best results).The data in each chunk are normalized via division by the respective best-fitting polynomial.Finally, the normalized chunks are reassembled to obtain the final light curve.

K2
The K2 mission was a follow-on to the Kepler mission [16] and was active between 4 February 2014 and 26 September 2018.Over that time span, there were 20 different pointings along the ecliptic where the two remaining reaction wheels and pressure from the solar wind were used to maintain the spacecraft pointing.The mission produced the same pixel-level and light-curve data products as the original Kepler mission.In addition, many versions of the light curves with various levels of additional processing were provided by various community-driven efforts and made available as High-Level Science Products on the MAST 3 .

TESS
NASA's Transiting Exoplanet Survey Satellite (TESS, see Ricker et al. [6]) has been conducting a photometric survey of most of the sky since the start of nominal science observations on 25 July 2018.TESS observes a 24 • × 96 • field (referred to as a "Sector") for ≈27 days at a time.During the two-year primary mission, the Southern Hemisphere and the Northern Hemisphere were each divided into 13 Sectors.Most of the Kepler field was observed in Sector 14 and again in Sector 26.A small number of targets were observed with a 2-minute cadence.Those that were not targeted for short cadence were observed at a cadence of 30 min in the Full Frame Images (FFIs).As part of the various mission extensions, the Kepler field was observed in Sectors 40, 41, 53, 54, and 55 with a 10 min cadence.For the objects that we discuss here, the TESS light curves were computed using the eleanor package of Feinstein et al. [17], which makes use of the TESScut package devised by Brasseur et al. [18].The light curves were detrended in a similar manner as the Kepler light curves described above.

Ground-Based Follow-Up Observations
Spectroscopic observations of KIC 5255552, 7668648, and 10319590 were obtained as part of the NOAO Survey program 2011A-0022 (P.I. A. Prša) as described in Kirk et al. [7].The echelle spectrograph attached to the Mayall 4 m telescope at the Kitt Peak National Observatory was used, which provided a spectral resolving power of R ≈ 20, 000 and a wavelength coverage between 4600 and 9050 Å, although only the wavelength region between 4875 and 5858 Å was used to extract the radial velocities.The typical observing sequence consisted of back-to-back exposures of around 900 s, with calibration exposures of a ThAr lamp at each telescope position.The IRAF (Image Reduction and Analysis Facility, version V2.17) software package [19,20] was used to perform the basic CCD calibration tasks and to extract, normalize, and merge the final spectra.The typical signal-to-noise ratios were usually between about 15 and 50 per resolution element.
The radial velocities were computed using the "broadening function" (BF) technique devised by Rucinski [21,22], as this technique generally gives better results than a simple cross-correlation analysis for double-lined spectra.Prior to the analysis, the merged spectra were normalized to their continua using the "continuum" task in the IRAF package.The wavelength range of 4875 to 5850 Å was used to compute the radial velocities, as this region avoids strong Balmer lines and strong interstellar and telluric lines, and contains numerous metal lines.High signal-to-noise spectra of various bright radial velocity standard stars were used as templates.Gaussian functions were fitted to the BFs to determine their centroids and widths.

Times of Eclipse, Transit, and Occultation Events
After performing the final detrending step on the light curves as described above, we are left with normalized chunks of data centered on individual eclipse, transit, and occulta-tion events.To find the mid-times of these events, we ideally want to have a smooth curve that has the same shape as the observed eclipse profiles and "slide" it across the data in time until the best match is found.Gaussian functions or simple polynomials will not exactly match the profile shape, especially when the signal-to-noise ratio is large.Instead, we used the ELC (Eclipsing Light Curve) code of Orosz & Hauschildt [23] to produce realistic model curves and its Markov sampler based on the Differential Evolution Monte Carlo Markov Chain (DE-MCMC) algorithm of Ter Braak [24] to sample parameter space and to find the uncertainties in the fitted times.To model an eclipse event, one needs to specify the time of conjunction, which obviously determines the time of mid-eclipse.This parameter can be given a narrow prior.The shape of the eclipse event is determined by the orbital period (for fixed masses, this sets the scale for the duration), the inclination, the radius of the eclipsed star, the ratio of radii, and limb-darkening parameters, and if necessary, the eccentricity parameters.These parameters can be given wide priors since for this purpose, we want a model with an eclipse profile that matches the observed eclipse profile, and we do not care if the properties of the stars themselves are realistic.Thirty-two chains are typically evolved for 100,000 generations, and posterior samples are taken every 5000 generations.The final adopted time of eclipse is taken to be the median of the posterior sample and its uncertainty to be the larger of |T 84% − T median | or |T 16% − T median |, where T 84% and T 16% is the times at which 84% and 16% of the samples are smaller, respectively.

The Photodynamical Model
For this discussion, a "photodynamical model" is a model that combines a numerical solution to the appropriate gravitational equations of motion of a system with two or more bodies, along with a detailed scheme to synthesize light curves, velocity curves, and other observable properties of the system.Our photodynamical model is implemented in our ELC code [23].Many details of the numerical scheme used to solve the equations of motion are given in Orosz et al. [15].We summarize some of the key features of the model and describe some recent modifications.
The numerical integrator uses Cartesian coordinates relative to the system's barycenter.Given the masses of each body and the Keplerian elements of each orbit (e.g., the period P, the time of barycenter conjunction T conj , the eccentricity e, the argument of periastron ω, the inclination i, and the nodal angle Ω) at a specified reference time T ref , the Cartesian coordinates (where x and y are in the sky plane and z is the radial coordinate with the positive z-axis is pointed towards the observer) are computed using algorithms given in Murray & Dermott [25].Three basic configurations are currently available.The first configuration consists of primary central body with up to 9 companions, each in successively larger orbits using Jacobi coordinates.There are no restrictions on the masses of the bodies (e.g., the central body does not have to be the most massive).This configuration includes the 2 + 1 mode where a close binary has a more distant companion.The second configuration is the "double binary" or 2 + 2 mode with 4 bodies consisting of two close binaries in a much larger orbit about the system barycenter.The third configuration is the "triple plus binary" or "3 + 2" mode with 5 bodies that has a close binary with a tertiary companion in a wider orbit paired with another close binary.The standard Newtonian equations of motion are supplemented with extra force terms to account for precession due to General Relativity (GR) and/or precession due to tidal bulges on the stars [26].There is also a Cartesian-to-Keplerian converter so that one may track changes to the orbital elements over time.
The numerical solution of the equations of motion essentially yields a function that gives either the sky coordinates of every body as a function of time or the radial coordinates of every body as a function of time (as discussed in Orosz et al. [15] these coordinates are corrected for light travel time).The time derivatives of the coordinates as a function of time are also returned.The times of eclipses, transits, and occultations are taken to be the times when the projected separation between the centers of two given bodies on the sky plane is minimized.In a similar manner, the times of ingress and egress of the various eclipse, transit, or occultation events are taken to be the times when the projected separation between the centers of two given bodies on the sky plane is equal to the sum of the respective radii of the two bodies.
The stars are assumed to be spherical for the photodynamical mode of ELC.The algorithm described by Short et al. [27] is used to compute the light curve.This algorithm allows for several of the common limb-darkening approximations to be used [e.g., the "quadratic" law, the "square-root" law, the "Power-2" law (the one used in this work), and the "logarithmic" law], and can account for more than two overlapping bodies (e.g., syzygy events).
ELC has several optimizers and Markov samplers to explore the parameter space.The ones used in this work include (i) the genetic algorithm based on the code given in Charbonneau [28], (ii) a simple Monte Carlo Markov Chain (MCMC) algorithm outlined by Tegmark et al. [29], (iii) the Differential Evolution MCMC (DE-MCMC) algorithm of Ter Braak [24], (iv) an implementation of the Nested Sampling algorithm described by Skilling [30], and (v) a recently developed "hybrid" code dubbed "gennestELC" that combines elements of the genetic and the Nested Sampling algorithms.Given an Ndimensional parameter space, scaled to the unit hypercube, one chooses N live "live points" and computes their likelihood values.The genetic algorithm is used to evolve the live points for a user-adjustable number of generations.After this evolution a bounding hyper-ellipsoid is computed, and a new set of live points is picked by uniform sampling from the hyper-ellipsoid.In the Nested Sampling algorithm, the volume of the hyperellipsoid is scaled at a rate proportional to 1/N live (see Skilling [30]) so that the actual volume (its "statistical volume") is much larger than the actual volume.In the gennestELC algorithm, the volume of the hyper-ellipsoid is allowed to shrink faster than the statistical rate but slower than the actual rate without any scaling.This algorithm can find good solutions in complex parameter spaces much faster than the standard Nested Sampling algorithm can when the number of free parameters is more than about 25.
All the algorithms discussed above are bounded in the sense that the user must select a lower bound and an upper bound for each fitting parameter.Obviously, care must be taken so that the "true" value of the parameter is between the lower and upper bound.It is often more efficient to fit various combinations of two parameters rather than the parameters themselves.For example, one can fit for √ e cos ω and √ e sin ω rather than for e and ω, or the sum and difference of the radii of the stars in the binary, and so on.
The likelihood is usually the standard χ 2 (one can also use the average deviation as a likelihood).The time-dependent observables that one usually fits for are the light curves (where up to 8 different bandpasses can be fit simultaneously), the radial velocity curves for up to 5 bodies and the times of eclipses, transits, and occultations.Observable properties that are not usually time-dependent, such as the observed rotational velocity of the stars, their spectroscopically measured surface gravities, etc. can also be used as extra constraints.
Finally, the ability to use the MIST [MESA (Modules for Experiments in Stellar Astrophysics) Isochrones and Stellar Tracks] models [31,32] in the likelihood has recently been implemented.For a large grid of various metallicities and ages, the isochrones supply the stellar radius and temperature as a function of the stellar mass.For each isochrone, the distance in a mass-radius-temperature cube from each star to the model curve is computed and summed.The distance is computed for each model in the grid, and the overall smallest summed distance d iso,min is found.The age and metallicity of the best-fitting isochrone are also saved.Finally, d iso,min is scaled by a weight W iso and the product is added to the overall χ 2 .When W iso is small (say less than 0.05), the isochrones do not really influence the fit since only a small quantity is added to the overall χ 2 .In this case, the isochrones can be used after the fitting is carried out as an independent check on the overall consistency and plausibility of the fit.In cases where some stellar parameters may not be particularly well constrained (for example, a three-body system where the third star does not eclipse or transit), the weighting factor W iso can be much larger (say 20) to ensure the stellar masses, radii, and temperatures are all simultaneously consistent with evolutionary models.At this point, no attempt is made to interpolate between models in the mass-radius-temperature cube.Finally, the user can choose to use only the mass-radius plane and ignore the temperatures.

The Systems of Interest
4.1.KIC 10319590 4.1.1.Overview KIC 10319590 was first identified as an eclipsing binary with a period of 21.34 days by Prša et al. [8] using the first 44 days of Kepler data.Slawson et al. [9] found significant eclipse time variations (ETVs) that suggested a dynamical origin (as opposed to the light travel time effect or LTTE).KIC 10319590 was among the more remarkable systems of the 39 triple-star candidates found by Rappaport et al. [12] in that the binary stopped eclipsing entirely after the first ≈400 days of Kepler observations (see Figure 1).Borkovits et al. [10] performed an analytic study of the ETVs of KIC 10319590 and found an outer period of 451 days, orbital eccentricities of 0.026(1) and 0.17(1) for the inner and outer orbits, respectively, and a retrograde outer orbit with a mutual inclination of 135.4 ± 0.3 • .In subsequent work, Borkovits et al. [11] found a mutual inclination of 40.2(4) degrees for the outer orbit.KIC 10319590 was observed nine times using the 4 m telescope at Kitt Peak and the echelle spectrograph between 28 June 2011 and 18 June 2013.A spectrum of the F8V star HD 174912 was used as the template for the BF analysis.The spectra are double-lined (e.g., two strong peaks are detected), with the strongest set of lines attributable to the primary star in the inner binary.The second set of lines does not belong to the secondary in the inner binary but rather to the tertiary star.The two peaks were blended in the spectra taken on 5 June 2012, 15 June 2013, and 18 June 2013.A double-Gaussian model was fitted to the remaining six BFs to determine the radial velocities (see Figure 2), after suitable heliocentric corrections were applied (see Table A6).Using the areas under the BF peaks as a proxy for the flux ratio [22,33], we find a flux ratio (in this case, tertiary flux divided by primary flux) in the spectral region between 4875Å and 5850Å of 0.18 ± 0.05.Selected Broadening Functions (BFs) of KIC 10319590 (filled circles) and their respective best-fitting double-Gaussian models (solid lines).The stronger peak is due to the primary star in the inner binary, and the other peak is due to the tertiary.
The photodynamical model has 30 free parameters.There are orbital parameters for the inner binary orbit, namely the period P bin , a time of primary conjunction T conj,bin , the eccentricity parameters e bin cos ω bin , e bin sin ω bin , and the inclination i bin .The nodal angle of the binary is fixed at Ω bin = 0.There are similar parameters for the outer orbit: P 3 , T conj,3 , e 3 cos ω 3 , e 3 sin ω 3 , i 3 , and Ω 3 .The masses are parameterized by M 1 , Q bin ≡ M 2 /M 1 , and Q 3 ≡ (M 1 + M 2 )/M 3 .The radii are parameterized by R 1 , ρ bin ≡ R 1 /R 2 and ρ 3 ≡ R 1 /R 3 .The stellar effective temperatures are T 1 , τ bin ≡ T 2 /T 1 , and T 3 .The power-2 limb-darkening law with the Maxted transformation for the coefficients was used (see [34]), and the limb-darkening coefficients h 1 and h 2 for the primary and for the secondary were allowed to vary.The apsidal motion constants k 2,1 and k 2,2 for the primary and secondary star, respectively, were free parameters.Finally, there are four additional free parameters to account for "seasonal contamination": s 0 for Quarters 1, 5, 9, 13, and 17; s 1 for Quarters 2, 6, 10, and 14; s 2 for Quarters 3, 7, 11, and 15; and s 3 for Quarters 4, 8, 12, and 16.The reference time for the osculating parameters is BJD 2,454,800.We fit the Kepler light curve and the radial velocity curve.We also include the measured times of all primary and secondary eclipses in the likelihood function.
We used the genetic algorithm optimizer for the initial exploration of parameter space.Once a good solution was found, the DE-MCMC code was run 10 separate times for 60,000 generations each, where each run had a different seed for the random number generator.After a burn-in period of 2500 generations, posterior samples were generated by sampling models every 3750 generations.The ten individual posterior samples were combined, yielding sample sizes of 13,600 for each fitting parameter and derived parameter.Table 2 gives the median values of the posterior distributions for the fitting parameters and the corresponding 1σ uncertainty.Except for the apsidal constants k 1,2 and k 2,2 the fitted and derived parameters are reasonably well constrained.Table 3 gives the median values of the posterior distributions and the corresponding uncertainties of several derived parameters.The initial Keplerian and Cartesian initial conditions for the best-fitting model are given in Table A9 in the Appendix C.  The best-fitting model light curve and the residuals of the fit are shown in Figure 1 and the best-fitting radial velocity curve and the residuals are shown in Figure 3. Figure 4 shows the Common-Period O-C (CPOC) diagram with the measured eclipse times and the best-fitting model times.Both types of eclipses were not measurable in the Kepler data past roughly day BJD 2,455,360.In the best-fitting model, the last primary eclipse is on day BJD 2,455,520.1,and the secondary eclipse is on day BJD 2,455,381.6.

KIC 7668648 4.2.1. Overview
KIC 7668648 was first noted in the catalog of Prša et al. [8].Using ≈125 days of Kepler data, Slawson et al. [9] showed that there were significant ETVs (≈±20 min).This object was one of 39 triple-candidate systems identified by Rappaport et al. [12], who performed a basic analytic analysis of the eclipse times that determined periods of P 1 = 27.8184days and P 2 = 204.7 days and masses of 0.02 ≤ M 3 ≤ 0.41 M ⊙ for the third star and 0.22 ≤ M bin ≤ 3.65 M ⊙ for the binary, both at 90% confidence.Borkovits et al. 2015 andBorkovits et al. 2016 [10,11] noted the presence of extra eclipse events (although they did not display them) and did a much more comprehensive analytical analysis of the eclipse times found from the full Kepler light curve.The latter reported masses of M 3 = 0.27(8) M ⊙ for the third star and M bin = 1.6(5)M ⊙ for the binary, where the 1σ uncertainty in the last digit is given in the parentheses.They also found a mutual inclination (e.g., the angle between the orbital plane of the binary and the plane of the outer orbit) of 42(1) degrees.Zhang et al. 2017 [35] showed a total of 11 transit events (their Figure 3) that are occultations and transits of the third star.They also show six additional possible transits (their Figure 4), but these appear to be spurious.Finally, Zhang et al. 2017 [35] attempted to model the binary eclipses without accounting for dynamical effects and, as a result, were unable to come up with a comprehensive and self-consistent model.Indeed, Borkovits et al. 2016 [11] noted that for systems like KIC 7668648 with large ETVs and extra eclipse events ". . . the accurate interpretation of such a system can be carried out only by simultaneous modeling of its photometric and dynamical properties, . . .".We provide such a model below.

Photodynamical Model
KIC 7668648 was observed in Long Cadence during Quarters 1 to 16 of the nominal Kepler mission.This target was placed on the list of Short Cadence targets starting in Quarter 17 and a little over one month of data was taken before the nominal mission operations ceased owing to a failed gyroscope.The full Kepler light curve is shown in Figure 5.There is a remarkable change in the eclipse depths where the eclipses went from being ≈15% deep at the start of the Kepler mission to being ≈50% deep at the end.TESS observed KIC 7668648 in Sectors 14 and 26 at 30 min cadence, and Sectors 40, 41, 53, 54, and 55 at 10 min cadence.The source, also known as TIC 158215322, is relatively faint with a TESS magnitude of 14.8, and consequently, the light curves extracted with eleanor have relatively low signal-to-noise.Three primary eclipses and three secondary eclipses were identified.Table A2 in the Appendix A gives the eclipse times derived from the Kepler and TESS data.KIC 7668648 was observed seven times using the 4m telescope at Kitt Peak and the echelle spectrograph between June 6, 2012 and June 16, 2013.A spectrum of the G8V star HD 181655 was used as the template for the BF analysis.The spectra are double-lined, with the strongest set of lines attributable to the primary star in the inner binary, while the other set of lines is attributable to the secondary star in the inner binary.A double-Gaussian model was fitted to the seven BFs to determine the radial velocities (see Figure 6), after suitable heliocentric corrections were applied (see Table A5).Using the areas under the BF peaks as a proxy for the flux ratio (in this case, secondary flux divided by primary flux), we find a flux ratio in the spectral region between 4875Å and 5850Å of 0.75 ± 0.05.Selected Broadening Functions (BFs) of KIC 7668648 (filled circles) and their respective best-fitting double-Gaussian models (solid lines).The stronger peak is due to the primary star in the inner binary, and the other peak is due to the secondary star in the inner binary.
The photodynamical model has 30 free parameters and is very similar to the one used for KIC 10319590, with three changes: (i) tides have a minimal effect on the binary, so only the extra force terms related to GR were used; (ii) the radii of the primary and secondary stars were parameterized by their fractional radii, where f 1 ≡ R 1 /a and f 2 ≡ R 2 /a; and (iii) the Power-2 limb-darkening coefficients for the third star h 1 and h 2 were free parameters.The reference time for the osculating parameters is BJD 2,454,750.We fit the Kepler light curve and the radial velocity curve.We also include the measured times of all eclipse events in the likelihood function.There are some pretty significant instrumental changes in the eclipse depths seen in the TESS light curve, so consequently, only the measured eclipse times were used in the fitting.
The fitting procedure for KIC 7668648 was very similar to the fitting procedure for KIC 10319590 described above.Once a good initial solution was found, the DE-MCMC code was run 10 separate times for 30,500 generations, where each run of the DE-MCMC code had a different seed for the random number generator.After a burn-in period of 1250 generations, posterior samples were generated by sampling models every 750 generations.The ten individual posterior samples were combined to yield overall sample sizes of 32,000 for each fitting parameter and derived parameter.The median values of the posterior distributions and the corresponding 1σ uncertainties for the fitting parameters are given in Table 4. Table 3 gives the median values of the posterior distributions and the corresponding uncertainties of some key derived parameters.The parameters are generally well constrained, and many have very small formal errors.For example, the stellar masses have uncertainties of about 1% and the stellar radii have uncertainties of about 0.4%.The initial Keplerian and Cartesian initial conditions for the best-fitting model are given in Table A8 in the Appendix C. The best-fitting model radial velocity curve models are shown in Figure 7.The bestfitting light-curve model is displayed in Figure 5.As noted before, there is a striking change in the depths of the eclipses from the start of the Kepler mission to the end of its nominal operation.In addition, there is another feature that we have no easy explanation for, namely, the roles of the primary and secondary changed between the beginning of the light curve and the end.This is shown in Figure 8.The first eclipse event seen in the light curve (near day −23, labeled P 1 in the figure) is deeper than the next eclipse event (near day −9, labeled S 1 in the figure), so the former was called a "primary" eclipse and the latter was called a "secondary" eclipse.The last pair of eclipse events seen in the Kepler light curve occur near days 1396.5 (labeled P 52 in the figure) and 1410 (labeled S 52 in the figure).There are 51 orbital cycles between events P 1 and P 52 and between events S 1 and S 52 .Thus, one would have to label event P 52 as a "primary" eclipse and event S 52 as a secondary eclipse, even though the latter is clearly deeper than the former.  of KIC 7668648 (these were observed on Long Cadence).Bottom: Close-up views of a pair of eclipses 51 cycles later (these were observed in Short Cadence).Please note that both events are overall much deeper than the first pair of events.In addition, the "secondary" eclipse at cycle 52 is deeper than the corresponding "primary" eclipse.In each case, the observations are shown with the black points and the best-fitting model is the red line.The residuals of the fit are shown in the thin panels.
As noted earlier, there are numerous tertiary eclipse events seen in the light curve.Figure 9 shows the transits of the primary (top of the figure) and occultations of the tertiary star by the primary (bottom of the figure).Figure 10 shows the corresponding figure for the secondary+tertiary events.In general, the model fits these observed events well.Owing to the complex geometry brought about by the strong dynamical interactions between the tertiary star and the binary, the tertiary events can have a wide variety of impact parameters and, hence, different depths.Finally, KIC 7668648 has one of the more extreme ETV signals of any known binary, where the instantaneous period (e.g., the interval between consecutive primary eclipses or the interval between consecutive secondary eclipses) can vary by amounts as large as ≈5 h.For perspective, the average orbital period is around 27.8 days, so 5 h is ≈0.75% of the orbital period.Figure 11 shows the model and observed CPOC signals.As discussed below, the apsidal period is around 27 years, so about half of the apsidal cycle is covered between the beginning of the Kepler observations and the end of the TESS observations.The amplitudes of the primary and secondary CPOC signals are ≈13 h.In spite of the large ETVs, the residuals of the fits are less than 30 s for the Kepler observations.This object was first identified as a binary by Prša et al. [8] using the first 44 days of Kepler data.Using ≈125 days of Kepler data, Slawson et al. [9] were able to determine an orbital period of P = 32.46days.Using the entire Kepler data set spanning ≈1450 days, Borkovits et al. (2013) [10] performed an analytic study of the eclipse time variations (ETVs).They found an outer orbital period of P 2 ≈ 863 days and orbital eccentricies of e 1 ≈ 0.3 and e 2 ≈ 0.4.Finally, they noted the presence of extra eclipses but did not display them or discuss them in detail.Zhang et al. [36] noted that the three sets of extra eclipses were broadly consistent with being caused by a third body on the ≈863 day period identified by Borkovits, but were otherwise unable to find a comprehensive and self-consistent solution.Likewise, Getley et al. [37] performed a stability analysis of KIC 5255552 and attempted to model the light curves.There is a set of four extra eclipses that span ≈25 days, and these events proved to be especially challenging to model since a single star on an ≈863 day orbit would not linger near conjunction for that long.Getley et al. [37] speculated that KIC 5255552 might contain two binaries with only the first binary eclipsing.As we show below, that is indeed the case.

Photodynamical Model
KIC 5255552 was observed by the Kepler mission from the start of Kepler Quarter 1 through the end of Kepler Quarter 16.Module 3 of the Kepler detector failed on 9 January 2010, and unfortunately KIC 5255552 was not observed during Quarters 5, 9, 13, and 17 (see Figure 12).The source, also known as TIC 120316928, was observed by TESS in Sectors 14 and 26 at 30 min cadence and Sectors 40, 41, 53, and 54 at 10 min cadence.A total of four primary eclipses and five secondary eclipses appear in the TESS light curve.Table A1 in the Appendix A gives the eclipse times derived from the Kepler and TESS data.KIC 5255552 was observed 11 times using the 4 m telescope at Kitt Peak and the echelle spectrograph between 24 June 2011 and 14 June 2013.Two spectra from 15 June 2013 were unusable owing to extremely poor signal-to-noise ratios.A spectrum of the G8V star HD 181655 was used as the template for the BF analysis of the remaining nine spectra.Except for the spectrum from 13 June 2013, the BFs have one strong peak and a second peak that is considerably weaker.These two peaks can be attributed to the primary and secondary of the eclipsing binary, respectively.A double-Gaussian model was fitted to the eight BFs with double peaks to determine the radial velocities (see Figure 13), after suitable heliocentric corrections were applied (see Table A4).Using the areas under the BF peaks as a proxy for the flux ratio (in this case, secondary flux divided by primary flux), we find a flux ratio in the spectral region between 4875Å and 5850Å of 0.18 ± 0.04.
Figure 13.Selected Broadening Functions (BFs) of KIC 5255552 (filled circles) and their respective best-fitting double-Gaussian models (solid lines).The stronger peak is due to the primary star in the first binary and the weaker peak is due to the secondary star in the first binary.
After a considerable effort, we were unable to find a suitable fit to the light curve using the 2 + 1 model that was successfully used for KIC 10319590 and KIC 7668648.Thus, we instead used the "double binary" model, which has 49 free parameters.We have three orbits: (i) the orbit of the eclipsing binary (hereafter "binary #1"), with a period P bin,1 , a time of primary conjunction T bin,1 , eccentricity parameters e bin,1 cos ω bin,1 and e bin,1 sin ω bin,1 , inclination i bin,1 , and nodal angle Ω bin,1 ≡ 0; (ii) the orbit of the non-eclipsing binary (hereafter "binary #2") with similar Keplerian elements (in this case, the nodal angle is a free parameter); and (iii) the outer orbit with similar Keplerian elements (where again the nodal angle is a free parameter).The stellar masses for binary #1 are parameterized by the primary mass M 1 and the binary mass ratio Q bin,1 ≡ M 2 /M 1 .The stellar masses for binary #2 are parameterized by 4 .The stellar radii for binary #1 are parameterized by R 1 and ρ bin,1 ≡ R 1 /R 2 , and the stellar radii for binary #2 are parameterized by ρ 3 ≡ R 1 /R 3 and ρ 4 ≡ R 1 /R 4 .The stellar effective temperatures are T 1 and τ bin,1 ≡ T 2 /T 1 for binary #1 and T 3 and T 4 for binary #2.We have light-curve data in two bandpasses.Since all four stars are involved in eclipse events in the Kepler light curve, each star has a pair of limb-darkening coefficients h 1 and h 2 for the Power-2 limb-darkening law with the Maxted transformations.The TESS light curve only has eclipse events involving the stars in binary #1, and this gives a pair of limb-darkening coefficients h 1 and h 2 for each star in that binary.We also have the seasonal contamination parameters for the Kepler as described above and similar parameters for the TESS light curve: c 0 for the time interval between days 3680 and 3705 (in units of BJD −2,455,000), c 1 for the time interval between days 4010 and 4030, c 2 for the time interval between days 4415 and 4440, and c 3 for the time interval between days 4760 and 4795.As was the case for TIC 7668648, it was found that tides have a minimal effect, so only the extra force terms related to GR were used.We fit the Kepler light curve, the TESS light curve, and the radial velocity curves for the stars in binary #1, the measured eclipse times for the primary and secondary eclipses of binary #1, and for the measured times for the extra eclipse events (whose origins were not known at first).
Finding a solution for KIC 5255552 proved to be extremely challenging.The basic properties of binary #1 are roughly known since it is eclipsing.Likewise, some of the basic properties of the outer orbit are roughly known from a basic analysis of the eclipse times of binary #1.On the other hand, the geometrical properties of binary #2 are almost completely unknown, except for two observed constraints: (i) the binary apparently does not eclipse, which eliminates some small and complex regions of parameter space, and (ii) the two stars in binary #2 are less luminous than the secondary in binary #1 since only two peaks are seen in the BFs.We discovered that one cannot really treat the second binary as a point mass for detailed work.That is, two nearly identical models that have everything the same apart from the time of primary conjunction for the second binary would have differences in the eclipse times for the first binary that was much larger than the uncertainties in the measured times.This explains in part why fits with three-body models, in an effort to get "close" to the solution as a starting point for a four-body model, do not really work as well as one might have expected.Finally, as others have noted [36,37], the timing of the extra eclipse events depends sensitively on precise positions of the bodies involved, and small deviations in various orbital elements can result in a large change in the timing of the extra events.
We initially tried two strategies.First, if one ignored the extra eclipse events and tried to get a very good model for the eclipses in binary #1, then at some point, the model would have the extra eclipses at the correct times.How one weights the measured times of the extra events can also be altered, so by setting the uncertainties to 1 day instead of ≈0.004 days, the extra events are almost "ignored".Several runs with the genetic algorithm and later the hybrid genetic algorithm/nested sampling code were done, but without success.The other approach was to make the uncertainties on a few of the extra events very small, so that those events were always matched after gennestELC converged.This approach also did not really work since subsequent attempts to refine these solutions with the DE-MCMC code almost always resulted in little overall improvement.
In the end, the successful strategy for finding a solution was relatively straightforward but required an incredible brute-force computing effort.Several different searches were conducted with gennestELC where the prior ranges for the period of binary #2 and the time of conjunction of binary #2 were severely restricted.The uncertainties on the extra eclipse events were set to 1 day.Then, after many runs were completed over a wide range of possible orbital periods for binary #2, the solutions were ranked by the overall χ 2 , and the DE-MCMC code was deployed using the top-ranked models as starting places.A few of these starting models were "close" enough so that a solution that matched the light curve with all the extra eclipse events was found.
After a good solution was found, the DE-MCMC code was used in a manner similar to what was used for KIC 7668648 and KIC 10319590 as described above.The code was run 10 times with different initial random number seeds for a total of 30,000 generations.After a burn-in period of 4000 generations, posterior samples were taken every 2000 generations.The samples were combined, yielding sample sizes of 22,400 for each fitting parameter and each derived parameter.Table 5 gives the median values of the posterior distributions for the fitting parameters and their corresponding 1σ uncertainty.Table 6 gives the same quantities for selected derived parameters.The initial Keplerian and Cartesian initial conditions for the best-fitting model are given in Table A7 in the Appendix C. Figures 14-16 show close-up views of the Kepler light curve near the times when the extra eclipse events occur.Figure 14 shows the sequence of four extra events that span ≈ 25 days that Getley et al. [37] discussed.The first two events in this sequence are caused by the primary star of binary #2 transiting the secondary of binary #1.The last two events in the sequence are the secondary star in binary #2 transiting the secondary star in binary #1.The three extra eclipse events near day 780 shown in Figure 15 are occultations of the stars in binary #2 by stars in binary #1.Finally, the two extra eclipse events near day 1375 shown in Figure 16 are transits of the primary star in binary #1 by the secondary star in binary #2. Figure 17 shows the TESS data (which features four primary eclipses and five secondary eclipses of binary #1) and the best-fitting model.Figure 18 shows the radial velocities and best-fitting radial velocity curve for binary #1, and Figure 19 shows the CPOC diagram for binary #1.The individual signals have modulations of up to ≈5 h.In addition, the diverging signals indicate significant apsidal motion, where the rate is ≈0.218 degrees per cycle.EPIC 220204960 (hereafter EPIC 2202) was observed for ≈80 days in Campaign 8 of the K2 mission [16].Rappaport et al. [38] found that EPIC 2202 consists of two EBs with periods of 13.27 days and 14.41 days.They obtained follow-up spectroscopic observations and measured radial velocities for all four components.The two EBs individually show strong ETVs, and given all of this Rappaport et al. [38] concluded that EPIC 2202 is a compact quadruple system in a binary-binary configuration (a 2 + 2 similar to that of KIC 5255552) with all 4 stars having masses near 0.4 M ⊙ .They performed a series of numerical simulations and an analytic analysis of the outer orbit (which assumed roughly coplanar orbits) and concluded the outer orbital period is very likely between 300 and 500 days with a possibility of being as long as two to four years (e.g., ≈730 to 1460 days).Using the photometric and radial velocity data provided by Rappaport et al. [38], we perform our own independent photodynamical analysis with the goal of seeing what constraints, if any, we could place on the properties of the outer orbit.

Photodynamical Model
Saul Rappaport kindly sent us the processed K2 light curve they used in their paper.We normalized the flux to 1.0 using a simple spline fitting routine where data in the various eclipses were given zero weight.The radial velocity data were taken from Rappaport et al. [38] (their Table 3).In our notation, the primary of "binary #1" has the radial velocity measurements denoted by "Star A-1" in their Table 3 and the secondary of binary #1 has the radial velocity measurements denoted by "Star A-2" in their Table 3.The primary and secondary of "binary #2", respectively, have measurements from the columns labeled "Star B-1" and "Star B-2" from that same table.The uncertainties on all the radial velocity measurements were taken to be 3 km s −1 , with two exceptions.First, the radial velocities of all four stars are relatively close together from the observation on day 2713.9(in units of BJD −2,455,000), which indicates the cross-correlation peaks were even more blended than usual.Initial fits showed evidence of "peak pulling" for both binaries (see [21]).That is, the star with the larger redshift at that time had a negative residual, and the star with the small redshift had a positive residual.Thus, for that day, the uncertainties on the measurements for all four stars were set to 6 km s −1 .Likewise, the radial velocities for binary #2 from day 7698.0showed similar evidence of peak pulling, so the uncertainties on those measurements were set to 6 km s −1 .Binary #1 was near the quadrature phase on that day, so the radial velocity measurements for that binary appear to be unaffected.
The photodynamical model has 42 free parameters.The Keplerian elements for binary #1 and binary #2 have the same parameterization as they did for KIC 5255552.For the outer orbit, parameterization of the period, time of conjunction, and the eccentricity elements were the the same as they were for KIC 5255552.However, instead of fitting for the outer orbit's inclination i out and nodal angle Ω out , we fit for the sum and difference of those quantities: Σ out,i,Ω ≡ i out + Ω out , ∆ out,i,Ω ≡ i out − Ω out .For each binary, the masses were parameterized by the sum of the masses of each component and the mass ratio (e.g., Σ bin 1 ,mass ≡ M 1 + M 2 , Q bin,1 ≡ M 2 /M 1 , with similar quantities for binary #2).Likewise, the stellar radii for each binary were parameterized by their sum and ratio, as in for binary #2.The stellar temperatures are parameterized by T 1 , T 2 , T 3 , and T 4 .Each star has a pair of limb-darkening coefficients h 1 and h 2 for the Power-2 limb-darkening law with the Maxted transformations.The source is heavily blended with a nearby brighter star and a single contamination parameter c 0 was used for the entire light curve.Finally, extra force terms for GR and tidal effects were used, and consequently, each star had an apsidal constant as a free parameter: k 2,1 , k 2,2 , k 2,3 , and k 2,4 .
We again used a brute-force method to explore parameter space.Both binaries eclipse, so consequently, we have reasonably good constraints on the orbital periods, the times of conjunction, the eccentricity parameters, and the inclinations.The nodal angle of binary #1 was fixed at zero, and the nodal angle of binary #2 was given a prior range of −90 • ≤ Ω out ≤ 90 • .On the other hand, there are very few constraints on the outer orbit that are immediately obvious.The outer period probably cannot be less than ≈100 days or so or otherwise the system would be unstable.Since there are strong ETVs seen in both binaries, the outer period cannot be excessively long.Given all of this we adopt fairly large prior ranges for the parameters of the outer orbit: 1700 ≤ T out ≤ 4400 in units of BJD −2,455,000, 100 ≤ P out ≤ 3000 in units of days, −0.61 ≤ e out cos ω out ≤ 0.61, −0.61 ≤ e out sin ω out ≤ 0.61, 30 ≤ Σ out,i,Ω ≤ 290 in units of degrees, and −180 ≤ ∆ out,i,Ω ≤ 180 in units of degrees.The gennestELC code was run 22 separate times, each with different initial seeds for a random number generator.After convergence, the DE-MCMC code was run for 100,000 generations, using the best model from gennestELC as the basis for generating initial models.After a burn-in period of 25,000 generations, posterior samples were taken every 12,500th generation.The separate posterior samples were combined, yielding sample sizes of 20,608 for the fitting and derived parameters.
The best-fitting model light curve is shown in Figure 20 and the best-fitting model radial velocity curves are shown in Figure 21.The initial Keplerian and Cartesian initial conditions for the best-fitting model are given in Table A10 in the Appendix C.Many of the posterior distributions are complex, often with two or more modes.Consequently, when reporting values of the fitted and derived parameters (Tables 6 and 7, respectively) we quote the sample median P med (where P is a generic parameter) and a range P 16% to P 84% , where 16% of the sample is smaller than P 16% and 16% of the sample is larger than P 84% .The posterior distributions for some of the key orbital elements are shown in Figures 22-24.The posterior distributions for some derived quantities are shown in Figure 25 (the eccentricity e and the argument of periastron ω), Figure 26 (the stellar masses and radii), and Figure 27 (the mutual inclination between the orbital planes of the two binaries and the mutual inclination between binary #1 and the outer orbit).

Stellar Parameters and Age Constraints
We can place the stars in the four systems studied herein in the mass-radius plane and see what, if any, constraints can be placed on their ages (hereafter denoted as τ).For this analysis, we use the MIST models [31,32].For a given age (where log(τ/yr) ranges from 5.0 to 10.3 in 0.05 dex steps) and metallicity (parameterized by [Fe/H] ≡ log 10 (N Fe /N H ) * − log 10 (N Fe /N H ) ⊙ , which ranges from −4.00 to −2.00 in 0.50 dex steps and from −2.00 to +0.50 in 0.25 dex steps), various stellar parameters are computed as a function of the stellar mass.In this way, we can plot various model curves in the massradius plane, and assuming the stars in a given system have the same age and that no mass exchange has happened, the measured masses and radii of all the stars should be consistent with a single isochrone model in the ideal case.In KIC 5255552, all four stars have mass measurements with ≈2% uncertainties and radius measurements with uncertainties < 1%.All three stars in KIC 7668648 have measured masses with uncertainties of ≈1% and measured radii with formal uncertainties of ≈0.3%.Thus, taking these measurements and their uncertainties at face value, the constraints on the ages could potentially be relatively tight.Figure 28 shows the mass-radius plane for these two systems (top panels).A solar metallicity isochrone (e.g., [Fe/H] = 0) with log(τ/yr) = 9.6 (≈4 Gyr) provides a decent match to all four stars in KIC 5255552.For KIC 7668648, a solar metallicity isochrone with log(τ/yr) = 10.2 (≈16 Gyr) provides a good fit.However, for obvious reasons, this isochrone is not physically realistic.Likewise, an isochrone with a metallicity of [Fe/H] − 0.25 and an age of log(τ/yr) = 10.1 (≈12.6 Gyr) provides a good match in the mass-radius plane, but as before, the age is unrealistic.Isochrones with τ = 10 Gyr and a metallicity of [Fe/H] − 0.5 and with log(τ/yr) = 9.95 (≈8.9 Gyr) and a metallicity of [Fe/H] = −0.75give good matches in the mass-radius plane.For each panel, the observed masses and radii and their uncertainties for each star are indicated with the red crosses, and the isochrone ages are indicated by the various colored lines, with black being the youngest age, then blue, then green, and then finally red being the oldest age.The solar metallicity isochrone ([Fe/H] = 0) with an age of log(τ/yr) = 9.6 (≈4 Gyr) provides a close match for all four stars in KIC 5255552.KIC 7668648 appears to be old (τ ≈ 10 Gyr) and somewhat metal poor, where the metallicities greater than [Fe/H] − 0.5 give ages older than the galaxy.Solar metallicity isochrones with ages of log(τ/yr) = 9.75 (≈5.6 Gyr) and log(τ/yr) = 9.85 (≈7.1 Gyr) provide plausible matches in the mass-radius plane for KIC 10319590.The uncertainties on the stellar masses and radii for EPIC 2202 are too large to give meaningful constraints.
The masses of the three stars in KIC 10319590 have uncertainties of ≈4%.The two stars in the eclipsing binary have radii determined to ≈1% accuracy.Since there are no extra eclipse events involving the third star, its radius is poorly determined.The masses and radii of the four stars in EPIC 2202 are determined with accuracies of ≈10% and ≈5%, respectively.The bottom panels of Figure 28 show these two systems in the mass-radius plane.Solar metallicity isochrones with ages of log(τ/yr) = 9.75) (≈5.6 Gyr) and log(τ/yr) = 9.85 (≈7.1 Gyr) provide plausible matches in the mass-radius plane for KIC 10319590.The uncertainties on the stellar masses and radii for EPIC 2202 are too large to give any meaningful constraints on the age of the system.

Dynamical Properties
Although a thorough investigation of the dynamical properties of the four systems is beyond the scope of this paper, we can compute some basic quantities that illustrate how the orbits can change over time.Owing to gravitational interactions between the various stars in these systems, the orbits are not perfect Keplerian orbits, and the orbital "elements" (e.g., the eccentricity e, the argument of periastron ω, the inclination i, etc.) will change.We can use the dynamical integrator in ELC and its Cartesian-Keplerian converter to integrate a model over a long time scale and record the values of e cos ω, i, Ω, and the mutual inclinations the various orbital planes over time.A simple Fast Fourier Transform routine can be used to find periodicities in these curves.Hereafter, we will call the time it takes the argument of periastron ω to change by 360 • as the "apsidal period" (in practice, this is the main periodicity of the e cos ω or e sin ω curves).Likewise, we will call the most significant periodicity of the i vs. time curve or the Ω vs. time curve the "nodal period".After the equations of motion are solved over the requested time interval, the eclipse-finding routine can find the number of eclipse events (e.g., primary eclipses, secondary eclipses, transits of the tertiary star in front of the primary star, etc.) in that interval.The "eclipse fraction" is then defined as the number of actual eclipse events divided by the "expected" number, which is simply the time interval divided by the average orbital period over that interval.For the three Kepler systems that have well-constrained parameters, we used the initial conditions from all the posterior models and integrated them over long intervals to determine the periodicities of the various curves.Table 8 gives the results for the 2 + 1 systems KIC 7668648 and KIC 10319590, and Figure 29 shows these curves for the respective best-fitting models.Table 9 gives the results for the 2 + 2 systems KIC 5255552 and EPIC 2202, where for EPIC 2202, we only give rough results for the best-fitting model since the code that computes the various periodicities of models in the posterior sample does not reliably work when the parameters of the outer orbit are poorly constrained.Figure 30 shows the long-term behavior of the key dynamical elements for the best-fitting models for the 2 + 2 systems.For KIC 7668648, we see that the inclination of the two orbits and the nodal angles of the two orbits do not change by more than ≈3 • .Consequently, the mutual inclination between the two orbital planes is never more than 2.5 • , and the primary and second eclipse fractions are each 1.0.The various eclipse events involving the third star occur between ≈26% and ≈29 of the time.The nodal precession period of the inner binary (and of the outer orbit) is only 19.61793 ± 0.00045 years.The apsidal period of the inner binary is somewhat longer at 26.0619 ± 0.0021 years, while the apsidal period of the outer orbit is 55.3400 ± 0.0011 years.
For KIC 10319590, the change in the inclination of the inner binary over time is much more dramatic where the inclination is in the range of 70 • to 149 • .The nodal angle of the binary has a similarly large range of between 0 • and 70 • .On the other hand, those two quantities for the outer orbit have a much smaller range of around ±10 • .Since the nodal angle of the outer orbit hovers around 30 • , the mutual inclination between the two orbitals is never less than about 25 • .The nodal period is about 108 years, and the eclipse fraction is about 0.76.Given these two quantities, the current non-eclipsing state of the binary will last ≈0.24× 108 ≈ 26 years.Owing to the relatively large mutual inclination, there were no eclipse events involving the third star in the posterior sample.Finally, the way the argument of periastron of KIC 10319590 changes over time is more complicated than it is for KIC 7668648.For the inner binary of KIC 10319590, the apsidal period (e.g., the time for ω to change by 360 • ) is 349 ± 21 years.However, a smaller amplitude modulation with a period of around 55 years is also evident in the e cos ω curve.The apsidal period is ≈1600 years for the outer orbits.
For KIC 5255552, the various inclinations i and the nodal angles Ω vary by less than ≈6 • over the time frame investigated.As a result, the mutual inclination between the two binary orbital planes and the mutual inclinations between the outer orbit and the binary planes are always less than ≈4 • .As seen in Table 9, binary #1 in KIC 5255552 always eclipses, whereas the eclipse fractions of binary #2 are ≈0.16.Figures 14-16 show the various extra eclipse events seen in KIC 5255552.The secondary star in binary #2 transited both the primary star in binary #1 and the secondary star in binary #1.Surprisingly, these types of events are rare with a transit fraction of ≈0.2% (see Table 9).The secondary star in binary #1 transited the primary star in binary #2, and vice versa.Table 9 shows that these events are somewhat more common with transit fractions of ≈8% and ≈10%.Finally, the secondary star in binary #2 transited the secondary star in binary #1, and this type of event has a transit fraction of only ≈0.2%.Given the low transit fractions of the various observed extra eclipse and occultation events, it seems we were exceedingly lucky that these events occurred during the relatively short Kepler mission, and that KIC 5255552 was not in the field of view of the failed CCD module when the events occurred.
As noted above, the code that determines the periodicities of the various dynamical quantities using models in a posterior sample does not reliably work for cases like EPIC 2202, where the properties of the outer model are poorly constrained.Consequently, we give rough periodicities of the model whose initial conditions are given in Table A10 in the Appendix C. KIC 5255552 and KIC 7668648 have nearly coplanar orbital planes.On the other hand, the mutual inclination between the two orbital planes in KIC 10319590 is much larger and is always within ≈12 • of 40 • .There is a mechanism known as Kozai cycles with tidal friction, whereby the orbital separation of a binary can shrink due to the combined effects of secular perturbations from a distant companion and tidal friction that become more and more important as the fractional radii (e.g., R/a) of the stars increases with time [39][40][41].An expected statistical outcome of this process is that short-period binaries will have distant companions with a distributions of mutual inclinations between the two orbital planes strongly peaked near 40 • and 140 • .The present-day configuration of KIC 10319590 seems to be consistent with this expected outcome.Although the uncertainties are much larger for EPIC 2202, the distribution of the mutual inclinations between the orbital plane of the first binary and the orbital plane of the outer orbit has peaks near 40 • and 140 • .Even though EPIC 2202 is a 2 + 2 system and not a 2 + 1 system, could a similar Kozai cycles with tidal friction process be shrinking the binary orbits while widening the outer orbit?A thorough analytical and numerical study of this process in 2 + 2 systems would be interesting, and is certainly beyond the scope of the present work.

Comparison to Previous Studies
When comparing the results of a comprehensive photodynamical model to previous work, one must use care since many of the parameters change with time.Thus, a quantity like the orbital period at a specific reference time that can have a very small formal error could be many σ off when compared to a value found for a different reference time or a value with a different meaning (for example, a period that is some kind of an average over a long interval).Also, the previous studies of the three Kepler systems did not have access to radial velocity measurements.Nevertheless, it is still interesting and useful to make a few cursory comparisons.
KIC 10319590 is a relatively uncomplicated 2 + 1 system with a fairly large period ratio (P out /P bin = 21.4 at the reference time of BJD 2,454,800) and no extra eclipse events.On the other hand, the number of useful Kepler data is limited since the binary stopped eclipsing early in the Kepler mission.Borkovits et al. [11] found masses of 1.5 ± 0.2 M ⊙ for the binary and 0.8 ± 0.1 M ⊙ for the third star, respectively.We find a total mass for the binary of 1.86 ± 0.04 M ⊙ and a tertiary mass of 1.05 ± 0.04 M |odot .Both quantities are ≈2σ larger than the corresponding quantities found by Borkovits et al. [11].They found a mutual inclination between the two orbital planes of 40.2 ± 0.4 • .Again, taking into account that this quantity changes with time, their value compares favorably with our value of 40.84 ± 0.12 • .Regarding the long-term orbital dynamics, Borkovits et al. [11] found an apsidal period of 68 years and a nodal period of 110 years.Their apsidal period is somewhat different than our value of 46.7 ± 0.3 years.On the other hand, we find a nodal period of 107.7 ± 0.3 years, which is in reasonably good agreement with their value.
KIC 7668648 is a 2 + 1 system like KIC 10319590, but is much more compact (the period ratio is P out /P bin = 7.5 at the reference time of BJD 2,454,750, compared to 21.4 for KIC 10319590).Also, KIC 7668648 has extra eclipse events, therefore making the modeling more challenging.Borkovits et al. [11] found a total mass for the binary of 1.6 ± 0.5 M ⊙ , which agrees with our value of 1.64 ± 0.01 M ⊙ .They found a tertiary mass of 0.27 ± 0.08 M ⊙ , which agrees favorably with our value of 0.275 ± 0.003 M ⊙ .Regarding the long-term orbital dynamics, Borkovits et al. [11] found a mutual inclination between the two orbital planes of 42 ± 1 • .This value seems to be too large, given the presence of numerous tertiary eclipse and occultation events.We find I = 2.4440 ± 0.0007 • at the reference time of BJD 2,454,750.They found an apsidal period of 54 ± 6 years, compared to our value of 25.826 ± 0.005 years.They quote a nodal period for the binary of 25 years with no uncertainty.We find a nodal period for the binary of 19.7601 years with a small formal uncertainty of 0.0005 years (e.g., around 4 h).
Comparisons for KIC 5255552 are more fraught since it is actually a 2 + 2 system and not a 2 + 1 system.Using a model that assumed a 2 + 1 system, Borkovits et al. [11] found a binary mass of 1.7 ± 0.2 M ⊙ , which compares favorably with our value of 1.70 ± 0.02 M ⊙ .They found a tertiary mass of 0.7 ± 0.1 M ⊙ .In our case, we find a total mass for binary #2 of 0.99 ± 0.01 M ⊙ .
For EPIC 2202, Rappaport et al. [38] attempted a few different techniques to measure the masses and radii of the component stars, for example, by modeling each binary separately and attempting to model them jointly.Their mass measurements had relative uncertainties on the order of 10%, similar to what we find here.Their mass measurements and our mass measurements agree at the ≈1σ level.Their radius measurements also had relative uncertainties of around 10%, compared to ≈5% for the present work.Based on dynamical simulations, they concluded the outer orbital period is likely between about 300 and 500 days.The posterior sample for this quantity generated by the extensive MCMC analysis has peaks near 500 days and 1000 days, with a tail that extends out to around 3000 days.As Rappaport et al. pointed out, a few more radial velocity measurements and/or a few more eclipse observations would greatly constrain the properties of the outer orbit.

Summary
We presented full photodynamical models of the 2 + 1 systems KIC 7668648 and KIC 10319590.The stellar masses and radii in KIC 7668648 are tightly constrained.Based on an isochrone analysis, the system is very old and must be somewhat metal-poor.The system is compact (P out /P bin = 7.5), and the two orbital planes are within a few degrees of being coplanar.KIC 10319590 does not have extra eclipse events, so the stellar masses and radii are less well constrained.The mutual inclination of the two orbital planes oscillates around ≈40 • , consistent with a triple system that underwent Kozai cycles with tidal friction.KIC 5255552 turns out to be a 2 + 2 system where one of the binaries currently does not eclipse.We presented here, for the first time, a comprehensive model that explains the eclipse timing variations of the eclipsing binary and all the extra eclipse events.All three orbital planes are within ≈4 • of being coplanar.Finally, we performed an extensive MCMC analysis of EPIC 2202 to constrain the geometrical properties of the system.Using the presently available data, the mutual inclination between the orbital planes of the two eclipsing binaries is only weakly constrained, where there is a hint of a mutual inclination of around 10 • .On the other hand, the posterior distribution for the mutual inclination between the orbital plane of the shorter period eclipsing binary and the outer orbit has two distinct peaks near 30 • and 150 • , with values between 50 • and 100 • ruled out.Further investigation will be needed to see if this particular bimodal distribution is the result of a type of Kozai cycles with a tidal friction process that operates in 2 + 1 systems.Finally, additional radial velocity measurements and additional eclipse observations would improve the constraints on the various fitting and derived parameters for all four systems, but perhaps none more than those of EPIC 2202.

Appendix A. Eclipse Time Measurements
We present tables of the measured eclipse and transit times for the Kepler systems.

Figure 1 .
Figure 1.Top: The normalized Kepler light curve of KIC 10319590 (points) with the best-fitting model (red line).Bottom: The residuals of the fit.

Figure 2 .
Figure 2.Selected Broadening Functions (BFs) of KIC 10319590 (filled circles) and their respective best-fitting double-Gaussian models (solid lines).The stronger peak is due to the primary star in the inner binary, and the other peak is due to the tertiary.

Figure 3 .
Figure 3. Top: The radial velocity curve of the KIC 10319590 primary star (black points) with the best-fitting model (black line) and the radial velocity curve of the KIC 10319590 tertiary star (red points) with the best-fitting model (red line) Bottom: The residuals of the fit.

Figure 4 .
Figure 4. Top: The Common-Period Observed Minus Computed (CPOC) diagram for the KIC 10319590 primary eclipses (black points) and secondary eclipses (red triangles), with the respective best-fitting models (black and red lines).Bottom: The residuals of the fit.

Figure 5 .
Figure 5. Top: The normalized Kepler light curve of KIC 7668648 (points) with the best-fitting model (red line).Bottom: The residuals of the fit.Note that the last portion of the light curve near day 1400 is in short cadence (1 min "exposures").

Figure 6 .
Figure 6.Selected Broadening Functions (BFs) of KIC 7668648 (filled circles) and their respective best-fitting double-Gaussian models (solid lines).The stronger peak is due to the primary star in the inner binary, and the other peak is due to the secondary star in the inner binary.

Figure 7 .
Figure 7. Top: The radial velocity curve of the KIC 7668648 primary star (black points) with the best-fitting model (black line) and the radial velocity curve of the secondary star star (red points) with the best-fitting model (red line) Bottom: The residuals of the fit.

Figure 8 .
Figure 8. Top:Close-up views of the first pair of primary (top left) and secondary (top right) eclipses of KIC 7668648 (these were observed on Long Cadence).Bottom: Close-up views of a pair of eclipses 51 cycles later (these were observed in Short Cadence).Please note that both events are overall much deeper than the first pair of events.In addition, the "secondary" eclipse at cycle 52 is deeper than the corresponding "primary" eclipse.In each case, the observations are shown with the black points and the best-fitting model is the red line.The residuals of the fit are shown in the thin panels.

Figure 9 .
Figure 9. Top: Transits of the primary star by the tertiary star in KIC 7668648.The transit near day 296.3 is grazing.Bottom: Occultations of the third star by the primary star in KIC 7668648.In each case, the observations are shown with the black points and the best-fitting model is the red line.The residuals of the fit are shown in the thin panels.

Figure 10 .
Figure10.Top: Transits of the secondary star by the tertiary star in KIC 7668648.The transit near day 700.5 is blended with a secondary eclipse and the event near day 900 is grazing.Bottom: Occultations of the third star by the secondary star in KIC 7668648.The occultation near day 979 occurred near a secondary eclipse, and consequently, the occultation event is almost 4 days long.In each case, the observations are shown with the black points and the best-fitting model is the red line.The residuals of the fit are shown in the thin panels.

Figure 11 .
Figure 11.Top:The Common-Period Observed Minus Computed (CPOC) diagram for the KIC 7668648 primary eclipses (black points) and secondary eclipses (red points), with the respective best-fitting models (black and red lines).Bottom: The residuals of the fit.Note the different scales on the y-axis for the TESS measurements in the lower panel on the right.

Figure 12 .
Figure 12.Top: The normalized Kepler light curve of KIC 5255552 (points) with the best-fitting model (red line).Module 3 failed on 9 January 2010, and consequently, no data were obtained for KIC 5255552 during Quarters 5, 9, 13, and 17.Bottom: The residuals of the fit.

Figure 14 .
Figure 14.Top:A close-up of the normalized Kepler light curve of KIC 5255552 (points) with the best-fitting model (red line).Each eclipse is marked with a fraction of the form M/N, where M is the label of the body in front and N is the label of the body in the back during the event.For example, "2/1" means "body #2 is in front of body #1", which, in this case, is a primary eclipse of the first binary.The first two of the four extra eclipses are caused by the primary in the second (non-eclipsing) binary passing in front of the secondary in the first (eclipsing) binary.The other two extra eclipses are caused secondary in the second binary passing in front of the secondary in the first binary.Bottom: The residuals of the fit.

Figure 15 .
Figure 15.Similar to Figure 14, but for the second set of extra eclipses seen in KIC 5255552.Here, the first (eclipsing) binary is passing in front of the second (non-eclipsing) binary.

Figure 16 .
Figure16.Similar to Figure14, but for the third set of extra eclipses seen in KIC 5255552.As was the case for Figure14, the second (non-eclipsing) binary is passing in front of the first (eclipsing) binary.

Figure 18 .
Figure 18.Top: The radial velocity curve of the KIC 5255552 primary star (black points) with the best-fitting model (black line) and the radial velocity curve of the secondary star star (red points) with the best-fitting model (red line) Bottom: The residuals of the fit.

Figure 19 .
Figure 19.Top: The Common-Period Observed Minus Computed (CPOC) diagram for the KIC 5255552 primary eclipses (black points) and secondary eclipses (red triangles), with the respective best-fitting models (black and red lines).Bottom: The residuals of the fit.

Figure 20 .
Figure 20.Top: The normalized K2 light curve of EPIC 2202 (points) with the best-fitting model (red line).

Figure 21 .
Figure 21.Top: The radial velocity curve of the EPIC 2202 binary #1 primary star (black points) with the best-fitting model (black line) and the radial velocity curve of the binary #1 secondary star (red points) with the best-fitting model (red line).Second: The residuals of the fit.Third: Similar to the Top above, but for the EPIC 2202 binary #2.Bottom: The residuals of the fit.

Figure 22 .
Figure 22.Top: Posterior samples for EPIC 2202 showing, from left to right, the period of binary #1, the period of binary #2, and the period of the outer orbit.Bottom: Posterior samples for EPIC 2202 showing, from left to right, the time of primary conjunction of binary #1, the time of conjunction of binary #2, and the time of conjunction of the outer orbit.The units of the times are BJD minus the offset given in each panel.In each panel, the vertical dashed lines indicate the sample median, the two vertical dotted lines delineate the central region that contains 68% of the sample.

Figure 26 .
Figure 26.Similar to Figure 22, but for the stellar masses (top) and stellar radii (bottom).

Figure 27 .
Figure 27.Similar to Figure 22, but for the mutual inclination between the two binary orbital planes (left) and for the mutual inclination between the orbital plane of binary #1 and the orbital plane of the outer orbit (right).

Figure 28 .
Figure 28.MIST isochrone models in the mass-radius plane for KIC 5255552 (upper left), KIC 7668648 (upper right), KIC 10319590 (lower left) and EPIC 2202 (lower right).For each panel, the observed masses and radii and their uncertainties for each star are indicated with the red crosses, and the isochrone ages are indicated by the various colored lines, with black being the youngest age, then blue, then green, and then finally red being the oldest age.The solar metallicity isochrone ([Fe/H] = 0) with an age of log(τ/yr) = 9.6 (≈4 Gyr) provides a close match for all four stars in KIC 5255552.KIC 7668648 appears to be old (τ ≈ 10 Gyr) and somewhat metal poor, where the metallicities greater than [Fe/H] − 0.5 give ages older than the galaxy.Solar metallicity isochrones with ages of log(τ/yr) = 9.75 (≈5.6 Gyr) and log(τ/yr) = 9.85 (≈7.1 Gyr) provide plausible matches in the mass-radius plane for KIC 10319590.The uncertainties on the stellar masses and radii for EPIC 2202 are too large to give meaningful constraints.

Figure 29 .
Figure 29.The long-term evolution of key orbital parameters for KIC 7668648 (left panels) and KIC 10319590 (right panels).From top to bottom, the quantities shown are e cos ω, the inclination i in

Figure 30 .
Figure 30.The long-term evolution of key orbital parameters for KIC 5255552 (left panels) and EPIC 2202 (right panels).From top to bottom, the quantities shown are e cos ω, the inclination i in degrees, the nodal angle Ω in degrees, and the mutual inclination between the two orbital planes I bin−out in degrees.The black lines indicate the quantities for binary #1, and the red lines indicate the quantities for binary #2, and the blue lines indicate the quantities for the outer orbit.

Table 1 .
Observational Data Summary for Target Systems.
Figure 17.TESS data for KID 5255552 with the best-fitting models.See TableA1in the Appendix A for the eclipse times.The vertical scale on the panels with the residuals is ±0.1.