{'fraction_non_alphanumeric': 0.037523536268268626, 'fraction_numerical': 0.026248542992916703, 'mean_word_length': 4.620181405895692, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
TIC 454140642: A Compact, Coplanar, Quadruple-lined Quadruple Star System Consisting of Two Eclipsing Binaries

We report the discovery of a compact, coplanar, quadruply-lined, eclipsing quadruple star system from TESS data, TIC 454140642, also known as TYC 0074-01254-1. The target was first detected in Sector 5 with 30-min cadence in Full-Frame Images and then observed in Sector 32 with 2-min cadence. The light curve exhibits two sets of primary and secondary eclipses with periods of PA = 13.624 days (binary A) and PB = 10.393 days (binary B). Analysis of archival and follow-up data shows clear eclipse-timing variations and divergent radial velocities, indicating dynamical interactions between the two binaries and confirming that they form a gravitationally-bound quadruple system with a 2+2 hierarchy. The Aa+Ab binary, Ba+Bb binary, and A-B system are aligned with respect to each other within a fraction of a degree: the respective mutual orbital inclinations are 0.25 degrees (A vs B), 0.37 degrees (A vs A-B), and 0.47 degrees (B vs A-B). The A-B system has an orbital period of 432 days - the second shortest amongst confirmed quadruple systems - and an orbital eccentricity of 0.3.


INTRODUCTION
A few percent of F-and G-type stars are members of stellar quadruples, and the fraction likely increases with stellar mass (e.g Raghavan et al. 2010;De Rosa et al. 2014;Toonen et al. 2016;Moe & Di Stefano 2017;Tokovinin 2017). These systems are ideal laboratories to study stellar formation and evolution in a dynamicallycomplex environment where the constituent binary stars can be close enough to interact with each other. These interactions can occur multiple times over the lifetime of the system, involve e.g. Lidov-Kozai oscillations (Lidov 1962;Kozai 1962;Pejcha et al. 2013), mass-transfer, tidal and dynamical interactions, mergers (Perets & Fabrycky 2009, Hammers et al. 2021, or formation of tight binaries (e.g. Naoz & Fabrycky 2014, and references therein). Quadruple systems are important participants in stellar evolution as a potential source for Type I SNe (Fang et al. 2018), as well as interesting outcomes such as common-envelope evolution, interactions with circumbinary disks, etc.
Detection and detailed characterization of a bona-fide quadruple-star system can be challenging. For example, it is not uncommon for close eclipsing binary stars to exhibit deviations from strict periodicity and recent studies show that at least ∼ 7% of Kepler's eclipsing binary stars show such deviations within the 4-yearlong observing window of the Kepler prime mission (e.g. Rappaport et al. 2013, Orosz 2015, Borkovits et al. 2016. Such eclipse timing variations (ETVs) can be indicative of dynamical interactions between the eclipsing binary and additional star(s) on wider orbit(s), and provide an important detection method for triple, quadruple, and higher-order stellar systems, as well as for circumbinary planets (e.g. Borkovits et al. 2016, Welsh & Orosz 2018. Other pathways for the detection of quadruple (and higher-order multiples) are visual and spectroscopic observations (e.g. Raghavan et al. 2010, Tokovinin 2017, Pribulla et al. 2006), but these methods generally do not guarantee that the potentially ad-ditional star(s) are indeed gravitationally-bound to the known binary.
A target star exhibiting multiple sets of eclipses, ETVs, and per-EB divergent radial velocities provides a unique opportunity to not only confirm the multiplicity of the system but also measure the orbital and physical parameters of its constituents with exquisite precision (e.g. Carter et al. 2011;Borkovits et al. 2018Borkovits et al. , 2020aBorkovits et al. , 2021. We note that because the outer orbital period of quadruple systems is much longer than the periods of the inner binaries (for otherwise the system would become unstable), detecting such features requires continuous and long-baseline observations of the target. A combination of space-based photometric surveys, extensive ground-based surveys photometry, and dedicated spectroscopy follow-up is ideally-suited to provide such observations. This was demonstrated by the highlysuccessful Kepler mission which detected a number of EB systems exhibiting additional features indicative of either circumbinary planets or triple and higher-order stellar systems (e.g. Welsh & Orosz 2018;Kirk et al. 2015, Borkovits et al. 2016. The Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015) is equally capable of discovering multiply-eclipsing stellar systems, as already demonstrated by Borkovits et al. (2020bBorkovits et al. ( , 2021 for stellar triples and quadruples, and Powell et al. (2021) for sextuples. The target presented here, TIC 454140642, exhibits the typical 2+2 hierarchy for quadruple systems, illustrated in Fig. 1, and joins the small number of such confirmed, well-characterized systems -VW LMi (Pribulla et al. 2008), V994 Her (Zasche & Uhlař 2016), V482 Per (Torres et al. 2017), EPIC 220204960 (Rappaport et al. 2017, EPIC 219217635 (Borkovits et al. 2018), CzeV1731 (Zasche et al. 2020), TIC 278956474 (Rowden et al. 2020), and BG Ind (Borkovits et al. 2021).
We note that while quadruple star systems have lower occurrence rates compared to triple systems, the larger outer orbit necessary for stability in the latter systems is generally well beyond the 27-day TESS observation time in a single sector. Specifically, the long outer orbital period substantially reduces the probability that tertiary eclipses or occultations of a stellar triple will be identified in one sector of TESS data. While such triples can (and do) produce multiple tertiary eclipses or occultations during one conjunction of the outer orbit in a singe sector of TESS data (e.g. Borkovits et al. 2020b), detecting multiple such conjunctions in one sector of data is highly unlikely. In turn, this makes determining the outer period of such triples -and thus the overall configuration of the respective system -challenging. In contrast, quadruple systems do not need to be at a conjunction of the outer orbit to reveal their nature, and eclipses/occultations of their constituent binaries can (and do) easily occur in a single sector of TESS data.
Here we present the discovery of the quadruple star TIC 454140642 that exhibits four sets of eclipses as well as prominent, dynamically induced apsidal motion, and short-term four-body perturbations. At the time of writing, TIC 454140642 is the closest to coplanarity amongst the quadruple systems mentioned above, and has the second shortest outer period (after VW LMi). This paper is organized as follows. Section 2 outlines the detection of the target in TESS data. In Sections 3 and 4 we describe the analysis of the available archival data and our follow-up observations, respectively. In Section 5 we present our comprehensive spectro-photodynamical analysis of the system's parameters. We discuss the properties of the system in Section 6 and draw our conclusions in Section 7.

DETECTION
Our detection methodology is described in detail by Powell et al. (2021). Briefly, we found hundreds of thousands of eclipsing binaries using a neural network classifier on the TESS Full-Frame Image (FFI) light curves.
Through manual examination of the lightcurves identified by the neural network as EBs, we found many that showed eclipses belonging to multiple periods in the same light curve. Through photocenter analysis, we determined that the large majority of these lightcurves are consistent with a superposition of two EBs originating from unrelated target stars within the same TESS aperture. A fraction of these, however, demonstrated on-target photocenter for all detected eclipses, suggesting they originate from the same stellar system. TIC 454140642, a previously unknown EB, was one of the systems that we found early in our examination of TESS light curves. The nature of the system was initially quite challenging to decipher because all the eclipses are similar in depth and duration, in addition to the fact we only had one sector of data with each eclipse present twice in the lightcurve. Preliminary analysis indicated that there are two sets of primary and secondary eclipses. The TESS eleanor lightcurve (Feinstein et al. 2019) for Sectors 5 and 32 is shown in Fig. 2, highlighting the four sets of eclipses; there are no visible out-ofeclipse variations. In accordance with the designation system developed by the IAU, these binaries will be referred to as A (the brighter one with a period of 13.62 days, and components Aa and Ab) and B (the fainter one, with a period of 10.39 days and components Ba and Bb). The contamination due to nearby sources is low according to the TESS Input Catalog (0.001), and photocenter analysis confirmed that both sets of eclipses originate from the target (see Fig. 3), indicating a potential quadruple system. This motivated us to submit the target for 2-min cadence observations in Sector 32 (TESS DDT program 020 PI: Kostov), and also coordinate follow-up efforts with the TESS Follow-up Observing Program (TFOP). The parameters of the system are listed in Table 1.

PHOTOMETRIC DATA
There is extensive archival data on TIC 454140642 from the WASP and ASAS-SN surveys (Pollacco et al. 2006, Shappee et al. 2014, Kochanek et al. 2017, altogether covering more than 4000 days. The time span of the available photometric data is shown in Figure 4, which includes the TESS data as well as our follow-up photometry. The ASAS-SN data of TIC 454140642 provide coverage of about 2000 days and overlap with TESS data. The target was observed by WASP with both the 85-mm and the 200-mm lenses, although the SNR of the former is poor and thus for our analysis we use only the latter. All four sets of eclipses are clearly recovered in the ASAS-SN and WASP data as shown in Figure 5. Box     (Cutri et al. 2012), (7) GALEX-DR5 (GR5) (Bianchi et al. 2011) Least Squares" analysis (BLS, Kovács et al. 2002) analysis of the two datasets provides similar orbital periods for the two binaries, as indicated in the figure captions and listed in Table 2. We note that the periods listed in the table are derived from the photometry alone; the final orbital periods are derived from the comprehensive spectro-photodynamical model described in Section 5.1. Altogether, the available photometry phase-folded on the periods listed in Table 2 shows clear apsidal motion (see Fig. 6) and indicates dynamical interactions between the two binaries.

Spectroscopy
TIC 454140642 was monitored spectroscopically with the Tillinghast Reflector Echelle Spectrograph (TRES; Szentgyorgyi & Furész 2007;Fürész 2008) mounted on the 1.5m Tillinghast reflector at the Fred L. Whipple Observatory on Mount Hopkins (AZ). TRES is a benchmounted, fiber-fed instrument with a resolving power of R 44, 000, and a wavelength coverage of 3900-9100Å in 51 orders. In total we gathered 26 spectra between September of 2020 and February of 2021, with signal-tonoise ratios in the region of the Mg I b triplet (∼5185Å) ranging from 28 to 60 per resolution element of 6.8 km s −1 . Exposure times ranged between 450 and 2400 seconds. The spectra were extracted and reduced as per Buchhave et al. (2010), with wavelength solutions derived from bracketing thorium-argon lamp exposures.
The very first spectrum of TIC 454140642 showed four sets of lines, which we expected should represent the two components of two double-lined binaries in the system (see Figure 7). Radial velocities were measured using QUADCOR, a four-dimensional cross-correlation algorithm introduced by Torres et al. (2007), which uses four separate and possibly different templates matched to each star. In this case, however, the stars turn out to be quite similar to each other (see below), so we adopted the same template for the four components. We used a synthetic spectrum taken from a large library of calculated spectra based on model atmospheres by R. L. Kurucz (see Nordstroem et al. 1994;Latham et al. 2002), covering the region of the Mg I b triplet. The template parameters were T eff = 6250 K, log g = 4.5, [m/H] = 0.0, and no rotational broadening. The choice of zero rotational broadening was based on experiments with a range of values for v sin i, which indicated the broadening for all four stars is smaller than our spectral resolution. In addition to the velocities, we used QUADCOR to measure the flux ratios among the components at the mean wavelength of our observations. While approximate orbital periods were known from the photometry for the two binaries undergoing eclipses, at first it was not trivial to identify which set of lines corresponded to which component of which binary. In several of the spectra the lines were heavily blended, so the velocity determinations were not always possible for all four stars at this early stage, and were guided by the location of the peaks in plots of the one-dimensional crosscorrelation functions. The process was further compli-  cated by the degeneracy between the flux ratios and the velocity separations when the lines are blended, which made the velocities imprecise. Based on the more reliable flux ratios from the relatively few spectra with well separated lines, it was eventually realized that one star was brighter than the other three, another was the faintest, and the remaining two had similar brightness. Gathering the velocities for the brighter and fainter components, and sorting out the others based on the pre-liminary photometric ephemerides, we obtained initial orbital solutions for the two binaries based on the 18 observations through mid-December 2020. These fits showed considerable scatter, far in excess of the precision expected for the measurements. The primary and secondary residuals for the 13-day binary showed an obvious upward drift, and those of the 10-day binary a drift in the opposite direction (see Fig. 8). This was the first sure sign that the two binaries are gravitationally  Table 2), highlighting the prominent apsidal motion. The different datasets are vertically offset for clarity. Lower panels: Same as above but for Binary B.
bound to each other, and in hindsight, explains our initial difficulty in identifying the components of the two binaries.
Solving for this linear drift separately in each of the double-lined orbits improved the velocities significantly, and allowed us to determine approximate binary elements. In addition to the period (P ), the elements for each binary are the velocity semiamplitudes of the components (K), the orbital eccentricity (e) and argument of periastron of the secondary (ω), the common centerof-mass velocity of the primary and secondary (γ), a reference time of periastron passage (τ ), and a radial acceleration term (γ) for the binary center of mass. The orbital periods were better determined from the TESS photometry as well as the WASP and ASAS-SN archival data (see Section 3), and were held fixed. These elements enabled us to better predict the velocities for the four components at epochs of heavy line blending, which were then used as a starting point for a refined analysis of those 18 observations and subsequent ones with QUADCOR. We report our final velocities for all epochs in Table 3. Typical uncertainties for stars Aa, Ab, Ba, and Bb are 0.66, 0.73, 0.95, and 1.01 km s −1 , respectively, and are based on the rms scatter of the residuals from the orbits described above. The spectroscopic flux ratios we obtained for stars Ab, Ba, and Bb relative to the primary of the 13-day binary (Aa) are 0.85, 0.81, and 0.67 respectively, at the mean wavelength of our observations (∼5185Å).
The first-cut binary parameters determined from this analysis are given in Table 4, and served to set initial values for the global photodynamical analysis described in Section 5. The masses of the four stars in the quadruple system average 1.15 M (assuming sin i ≈ 1, as both binaries are eclipsing), and none of the stars has a mass that differs from this average by more than 0.05 M . Given that the total masses of the two binaries are rather similar, we note also that the straight average of the γ velocities, (γ A + γ B )/2 26 km s −1 , then reflects the approximate radial velocity of the system center of mass. Likewise, we expect that the radial accelerations for each binary will be nearly the same in mag- nitude, but of opposite sign -which they indeed are: Finally, we used the relative velocities and accelerations to derive a preliminary estimate of the outer orbital period and separation under the assumptions that the outer quadruple orbit (i) is being viewed nearly edge-on (i.e., i out 90 • ), and (ii) can be approximated as circular. The former assumption is based on the fact that the EB orbits are both viewed nearly edge-on (but we caution that this is not a proof of the overall 'flatness' of the system). As for the assumption of a circular outer orbit, it is actually e out 0.32, but we proceeded under the circular orbit assumption to obtain a quick and simple estimate of the scope of the outer orbit. There are two equations that can be solved to find the outer orbital separation, a out , and phase, φ out , at the time of the RV observations, where φ is the angle from inferior conjunction of binary A: (1) If we solve for a out and φ out , and then calculate the outer period, P out we find: where the uncertainties are due only to those associated with the measurements used in the calculations, and not to the actual eccentric nature of the outer orbit. These approximations for the outer orbit served to guide the full photodynamical modeling described later.

TESS Followup Observing Program
We acquired ground-based time-series follow-up photometry of TIC 454140642 as part of the TESS Followup Observing Program (TFOP) 1 . We used the TESS Transit Finder, which is a customized version of the Tapir software package (Jensen 2013), to schedule our eclipse observations. The photometric data were extracted using AstroImageJ (Collins et al. 2017). Observations were acquired using the Las Cumbres Observatory network (LCO; Brown et al. 2013) and the Hazelwood Private Observatory (Churchill, Victoria, Australia) as described in Table 5.

Other follow-up ground based photometry
In addition to the observations described above, we also used other photometric data for confirmation of the ETVs predicted from TESS data. These observations of eclipses of both pairs were obtained using the following instruments: • Ondřejov Observatory, Czech Republic, 65-cm reflecting telescope, using a G2-3200 CCD camera.
• Danish 1.54-m telescope located on La Silla in Chile, and using standard V photometric filter.

Speckle Imaging
The target star was placed on the Southern Astrophysical Research Telescope (SOAR) speckle program in 2020 December. It was observed on 2021 February 27. The star was unresolved and no companions with ∆I < 1 mag and having a separation greater than 0.048 were detected. The achieved detection limits are ∆I < 2.9 mag at 0.15 and ∆I < 4.4 mag at 1 . The details of these observations will be published jointly with other SOAR speckle results as part of the series of papers represented by .
The A-B pair, with an estimated semimajor axis of 5 mas, is not expected to be resolved at SOAR. However, the speckle imaging demonstrated the absence of other nearby faint stars that could otherwise be missed by the photometry or spectroscopy. To the best of our knowledge, then, TIC 454140642 is an isolated quadruple system.

Combined spectro-photodynamical analysis
We used the Lightcurvefactory software package (Borkovits et al. 2019(Borkovits et al. , 2020a to carry out a simultaneous, joint analysis of the lightcurves, ETV, RV curves and multi-passband SED data of the system. We used three separate lightcurves as follows: i) TESS Sectors 5 and 32 lightcurves. In the case of the Sector 32 data we used the 2-min cadence dataset for extracting the ETV data; however, for the joint analysis we binned these data into 1800 sec, resulting an effective exposure length of 1425sec, 2 to get the same sampling as was available for Sector 5 data. Furthermore, we excluded the flat out-of-eclipse lightcurve sections and kept only the narrow ±0. p 03-phase-domain regions around each eclipse.
ii) Similarly clipped regions of archival WASP data. After the removal of outlier data points, these data were averaged into 1800-sec bins for the analysis.
iii) The ground-based eclipse observations of our follow-up photometric campaign. These data were also binned at 1800 sec. South Africa 2020-11-09 z-short 6.8 ingress Notes. The z-short filter is the same as the Pan-STARRS Z-short pass-band, i.e. λcenter = 870 nm, λ width = 104 nm. The i filter is the same as the SDSS i pass-band, i.e. λcenter = 755 nm, λ width = 129 nm.     Table 3). The solid curves represent the lowest χ 2 spectro-photodynamical model solution (see in Sect. 5.1). The black curves show the RV contributions of the outer orbit. The residuals are also presented in the bottom panels. Note that the data obtained after JD 2 459 200 were not used for the analysis, although they agree perfectly with the solution.
Additionally, we used the four ETV curves (primary and secondary ETV data for both binaries; see Tables 6 and  7), and the first 18 epochs of the four TRES RV curves (see Table 3). Furthermore, the observed passband magnitudes tabulated in Table 1 were also used for the SED analysis. 3 The modelling runs proceeded as follows. Lightcurvefactory calculates the positions and velocities of each star at any requested instant via the numerical integration of the orbital motion. Then, taking into account the radii, temperatures and other atmospheric properties of each star, and their relative geometry with respect to the others and to the observer, the code predicts the net lightcurve of the whole system, enabling any kind of mutual (or multiple) eclipses. Note that in the present analysis, stellar radii, and temperatures were constrained through the stellar masses and the metallicity and age of the system with the use of built-in PARSEC stellar isochrones (Bressan et al. 2012), as was described in Borkovits et al. (2020a) in detail. Furthermore, the same PARSEC tables are used for generating theoretical passband magnitudes for the SED fitting part of the analysis. To solve the inverse problem, Lightcurvefactory employs a Markov Chain Monte Carlo (MCMC)-based parameter search, implementing the generic Metropolis-Hastings algorithm (see e.g. Ford 2005).
In most of the MCMC runs we adjusted the following 20 parameters: (i) 4 orbital elements for binary A, i.e., the eccentricity and argument of periastron in the combinations of (e sin ω) A and (e cos ω) A , its inclination, i A and the longitude of the ascending node Ω A relative to the node of binary B; (ii) 3 orbital elements for binary B, i.e., the same as above with the exception of Ω B , for which the value was kept at zero; (iii) 6 orbital parameters for the outer orbit, as its anomalistic period, P AB , eccentricity and argument of periastron, inclination, periastron passage time, τ AB and, the longitude of the ascending node relative to the node of binary B; 3 Similar to our previous works, for the SED analysis we used a minimum uncertainty of 0.03 mag for most of the observed passband magnitudes, in order to avoid the outsized contribution of the extremely precise Gaia magnitudes and also to counterbalance the uncertainties inherent in our interpolation method during the calculations of theoretical passband magnitudes that are part of the fitting process. The two exceptions are the WISE W 4 and GALEX F U V magnitudes, for which the uncertainties were inflated to 0.3 mag.
(iv) 4 mass-related parameters, as the masses of the two primaries, m Aa,Ba ), and the (inner) mass ratios of the two binaries, q A,B ; (v) 3 global parameters, as the (logarithmic) age of the system, its metallicity, [m/H], and the interstellar reddening E(B − V ).
Moreover, the following parameters were internally constrained: (i) 4 remaining orbital parameters of the inner orbits, i.e., their periods (P A,B ) and inferior conjunction times (T inf A,B ) at epoch t 0 , which were constrained through the ETV curves (see Borkovits et al. 2019); (ii) The systemic radial velocity of the entire quadruple system (γ) was obtained at each step by minimizing the χ 2 RV value a posteriori; (iii) 8 fundamental stellar parameters, i.e., the radii and temperatures of the four stars, which were interpolated from the PARSEC grids with the use of {mass, metallicity, age} triplets in the manner described by Borkovits et al. (2020a); (iv) The (photometric) distance of the system was also obtained through an a posteriori minimization of χ 2 SED ; (v) Finally, the (logarithmic) limb-darkening coefficients were calculated at each trial step internally from the pre-computed passband-dependent tables downloaded from the Phoebe 1.0 Legacy page 4 . These tables were originally used in former versions of the Phoebe software (Prša & Zwitter 2005). Table 8 lists the median values of the stellar and orbital parameters of the quadruple systems that are being either adjusted, internally constrained, or derived from the MCMC posteriors, together with their 1σ statistical uncertainties. The RV curves, lightcurves, ETV curves, and SED model of the lowest χ 2 global solution are plotted in Figs. 8-11.
In addition to the usual orbital elements tabulated in the first section of Table 8, we list the mutual (i.e. relative to each other) inclinations of the three orbital planes as well (i mutA−B;A−AB;B−AB ). We also provide the angular orbital elements in the dynamical frame of reference, i.e. the inclinations of the orbits relative to the invariable plane of the quadruple system (i dyn ), and the longitude of the pericenter ( dyn = Ω dyn +ω dyn ) of each orbit defined as the sum of the dynamical longitude of the node 5 and the dynamical argument of pericenter 6 . We note that due to the almost perfect coplanarity of the entire system, the intersections of the three orbits with the invariable plane (and therefore all the three pairs of Ω dyn and ω dyn ) become ill-defined; however, their sums, i.e. dyn = Ω dyn + ω dyn remain well-defined and we tabulate the latter quantities in Table 8. The relevance of these dynamical orbital elements lies in the fact that they, along with the mutual inclinations, appear directly in the analytic perturbation theories of multiple star systems (see, e.g. Borkovits et al. 2015, for further discussions). Furthermore, in this section we tabulate also some apsidal motion related parameters which will be discussed in Sect. 6.3. Finally, for completeness, we also provide the inclination and nodal angle of the invariable plane (i inv , Ω inv , respectively).
We note that the values listed in Table 8 represent the (instantaneous) osculating orbital elements at the (arbitrarily chosen) epoch t 0 = 2 458437.5. Therefore, the respective periods listed in the first row of the table should not be used for eclipse timing prediction for future follow-up observations. This is illustrated in Fig. 12 where we show the various relevant periods for the system. The osculating orbital elements (grey lines in the figure) are calculated at each instant from the Cartesian (Jacobian) coordinates and velocities of each orbit as if they represented pure, unperturbed, Keplerian twobody motions. In other words, these orbital elements describe the (osculating) two-body orbit on which the motion would continue if the perturbing forces vanished at that instant. The true motion occurs as the envelope of the osculating orbits. Naturally, in the presence of time-varying external forces such as those due to the gravitational perturbation from the second binary 7 , the osculating elements vary continuously, with a characteristic period which is half of the orbital period of the given binary. Furthermore, if the outer orbit happens to be eccentric (as in the present case), the perturbing forces -as well as the amplitudes of the variations of the osculating elements -naturally increase around the periastron passage of the outer orbit. As a result, the value of a given osculating element depends strongly on the 5 The angular distance of the intersection of the given orbital plane and the invariable plane -the dynamical nodal line -from the intersection of the plane of the sky and the invariable plane, measured along the invariable plane. 6 The angular distance of the pericenter from the dynamical ascending node along the orbital plane 7 We note that tidal forces may (and in several case should) be considered as well.
inner and outer orbital phase of the sampling. This is demonstrated in Fig. 12 by the green curve, which connects the anomalistic periods that were sampled during all the consecutive inner periastron passages of the given binary and, similarly, by the magenta curve which shows the osculating anomalistic periods during each primary eclipse. 8 However, none of the periods discussed thus far can be observed. The periods that can be observed are, e.g. the time between two consecutive primary (or, secondary) eclipses or, in case of the RV observations, the repeating time of the RV curve. As both the eclipses and the RV curve are periodic according to the true longitude-like quantity u = v + ω, where v denotes the true anomaly and ω stands for the argument of periastron in the observational frame of reference 9 , the measurable quantity is the 2π increment of u, instead of v, i.e. some kind of averaged sidereal (or, eclipsing) period, instead of an averaged anomalistic period. Therefore, in Fig. 12 we plot the time intervals between two consecutive primary (red) and secondary (blue) eclipses. One should be aware, however, that these observable eclipsing or, sidereal periods are systematically longer than the one-sidereal-period average of the instantaneous osculating anomalistic periods (brown symbols on the figure). The reason is that in a triple (and quadruple) system the clock of the inner binary is slowed as, on average, the perturbing forces from the outer bodies act as if the central mass of the binary were reduced. Finally, the black curve in Fig. 12 represents the average value of the eclipsing period which is practically the average of the red and blue curves, and can be determined most easily from the ETV curves as it is simply the linear term of the ETV ephemeris (as listed in the last row of Table 2, and in Fig. 9, as well).

PHOEBE analysis
For completeness, we also analyzed the TESS data via a classical light-curve modelling approach using the Wilson-Devinney method (Wilson & Devinney 1971) as implemented in the PHOEBE code (Prša & Zwitter 2005). Both binaries were analysed independently, using several simplifications for the analysis. Both orbital periods were kept fixed, and the mass ratio in each binary was set to the corresponding value given in Table 8. Likewise, the temperature of the primary in both binaries was taken from Table 8. See Table 9 for the results of the PHOEBE fit. These are in good agreement with those derived with the more sophisticated approach using the Lightcurvefactory code. Note, however, that the uncertainties associated with these PHOEBE results are statistical only, and have often been found to be underestimated.
6. DISCUSSION 6.1. Origin This hierarchical system has several distinctive features pointing to its origin, namely (i) nearly equal masses of all four stars, (ii) a compact outer orbit, and (iii) a near-perfect coplanarity. It is well known that close binaries cannot form directly with their presentday small separations owing to the so-called 'opacity limit' to fragmentation on the order of 10 au (e.g. Larson 1972). Accretion of gas on to a nascent binary shrinks its orbit and at the same time increases the mass ratio because the secondary star accretes more; this mechanism can explain statistics of the close-binary population, including inner subsystems in hierarchies (Tokovinin & Moe 2020). Preferential growth of the secondary component produces a subset of binaries with almost equal masses (twins) which is more prominent at short periods. A compact quadruple such as TIC 454140642 can be explained by the accretion-induced migration of both inner and outer subsystems, assuming that an initially wide and low-mass seed quadruple gained most of its mass by accretion. The seed quadruple could be produced, for example, through an encounter of two protostars A and B, each surrounded by a gas envelope. An encounter leads to a burst of accretion on to each star that produced the inner subsystems by disk fragmentation. However, other formation channels are also feasible. The important factor is a substantial mass growth of this progenitor quadruple by accretion that eventually created four stars with similar masses and shortened both inner and outer orbits. Dissipative interaction of the orbits with gas also led to the coplanarity of the entire system. This formation scenario of compact 2+2 quadruples makes a strong prediction of near-equal masses and orbit alignment. Indeed, other known compact quadruples share these properties with TIC 45414062. The closest such system, VW LMi (11029+3025, HD 95660), has an outer period of 355 days and inner periods of 0.478 and 7.93 days (Pribulla et al. 2020). The masses of the stars in the closest eclipsing pair are not equal, but this could be a result of the later mass transfer; the other, 7.9-day subsystem is a twin with q = 0.98, and the outer mass ratio is 0.93; all orbits are coplanar. Another compact and doubly-eclipsing quadruple V994 Her (= HD 170314) has an outer period of 2.9 yr and inner periods of 2.1 and 1.4 days (Zasche & Uhlař 2016). The inner mass ratios are 0.86 and 0.95, and the outer mass ratio is 0.67. There is another, more distant star in this system. Finally, the doubly-eclipsing quadruple OGLE BLG-ECL-145467, found by Zasche et al. (2019), has an outer period of 4.2 yr, with inner periods of 3.3 and 4.9 days, and the masses of all four stars are comparable.

Gaia EDR3 Astrometry
The Gaia EDR3 astrometric excess noise for TIC 454140642 is 0.144 mas at a significance of σ = 24.476; the number of good observations is 257 (Gaia EDR3, Gaia Collaboration 2020). The physical displacement corresponding to the excess noise is ∆a measured = 0.052 au (Eqn. 3, Belokurov et al. 2020). Taking into account the mass and luminosity ratios between binary A and B (0.94 and 0.8, respectively), as well as the eccentricity and argument of periastron of the AB system, the expected astrometric wobble is ∆a expected = 0.053 au (Eqn. 5, Belokurov et al. 2020) -practically identical to ∆a measured . In addition, using the methodology of Stassun & Torres (2021) and the measured renormalized unit weight error (RUWE) of 1.299, the photocenter semi-major axis is 0.27 mas which, at the measured parallax of 2.786 mas, corresponds to ≈ 0.1 au -comparable to ∆a measured . Altogether, this indicates that Gaia has seen the astrometric wobble of the AB system.

Long-term Dynamical Stability and Orbital Evolution
The orbital elements listed in Table 8 represent a snapshot of the system at the reference epoch. Due to the strong dynamical interactions between the two binary stars, these orbital elements, in accordance with the perturbation theory of hierarchical stellar systems (see, e. g., Harrington 1968) vary over time with at least three characteristic timescales: (i) the inner periods, (ii) the outer period and, (iii) the ratio of P 2 out /P in (see, e.g., Figs. 9, 12). The variations on the latter two timescales can nicely be confirmed on the ETV curves (Fig. 9), where, besides the strictly P out -period cyclic variations (not to be confused with the much smaller amplitude Rømer-delay, shown also in Fig. 9), the dynamically forced apsidal motion is also manifested.
To quantify the properties of the apsidal motions for all of the three orbits we calculated theoretical apsidal motion rates, ∆ω (i.e. the variation of the argument of periastron, ω, during one revolution of the given binary). For this approximate calculation, we used hierarchical three-body models. So, for example, to calculate the dynamical apsidal motion rate of binary A, we replaced the other binary (B) by a single outer body having mass  m B and, vice versa to compute the dynamical apsidal motion of binary B. Besides the dynamically forced apsidal motion, we also calculated the general relativistic and tidal contributions. For the two inner orbits we followed strictly the formulae and method described in Borkovits et al. (2015), Appendix C. For the outer orbit we considered the tidal contribution to be negligible and, furthermore, in the case of the dynamical term we simply calculated the algebraic sum of the contributions from the two separate 'triple systems' (formed by stars Aa,Ab,B and Ba,Bb,A). We list these apsidal advance rates in the 'apsidal motion parameters' part of Table 8. As one can see, the dynamical apsidal advance rate (∆ω 3b ; where '3b' stands for 'third body') is higher by 4-6 orders of magnitude than the relativistic and tidal rates (∆ω GR , ∆ω tide , respectively) for all the three orbits, therefore, these latter contributions are certainly negligible. Assuming that these apsidal motion rates re-main nearly constant (which is a plausible assumption for such a nearly coplanar system) we estimated approximate theoretical apsidal motion periods for the three orbits. These periods were found to be U theo = 87 ± 1, 109 ± 1, and 480 ± 1 years for binary A, B and AB, respectively.
In order to check the plausibility of our assumptions, we determined the 'true' apsidal motion periods of the three orbits. For this we numerically integrated the bestfit four-body model for ∼ 3 000 years, and then determined the average apsidal motion periods for all three orbits. We obtained U num = 70, 97, and 472 years, respectively. Therefore, we can conclude that our simple, theoretical three-body model slightly overestimated the true periods (i.e., in other words, underestimated the amplitude of the perturbations a bit); however, it resulted in a correct estimation of their magnitude.  Finally we note that the rapid apsidal motion of binary A, in addition to the readily visible divergence of the primary and secondary ETV curves, will also be detectable within a few years through the variation of the shape of the RV curve of binary A. This is illustrated visually in Fig 13. (Note, despite the similarly rapid apsidal advance of binary B, due to its almost circular orbit, the same effect will remain below the uncertainties of the currently available RV data.) We investigated also the short-term dynamical stability of the system. To confirm this we integrated the 4-body equations of motion for 5 million days (corresponding to ≈ 10, 000 outer periods), using the parameters from Table 8. We used the high-order IAS15 integrator with adaptive timesteps (Rein & Spiegel 2015), part of the REBOUND N-body package (Rein & Liu 2012). We provide a REBOUND jupyter notebook demonstrating the initialization and orbital element computation at https://github.com/vbkostov/TIC 454140642 10 As shown in Figure 16, there are no indications of chaotic motion and the system remains stable for the duration of said integration. The semi-major axis of binary A (a A ) oscillates between 31.8 R and 32 R 10 The orbital elements for binaries A and B represent the corresponding osculating two-body orbits, while the elements for A-B represent the orbit of the centers of mass of binaries A and B around one another. Despite the near edge-on orientation of the system, there are no mutual eclipses or occultations between the constituents stars of binary A and binary B. This is highlighted in Figure 17 where we show the evolution of the respective impact parameters between the four stars. While there are periodic, precession-induced variations in the corresponding minimum impact parameters, the latter never reach below ≈ 4.

SUMMARY
We have presented the discovery of the quadruplelined, doubly eclipsing quadruple system TIC 454140642. The target was observed by TESS in Sector 5 (at 30-min cadence) and again in Sector 32 (in 10-min cadence), and exhibited four sets of eclipses associated with two binary stars (A and B). Additional spectroscopic and photometric measurements reveal divergent radial velocities, eclipse-time variations and apsidal motion for both binaries, confirming that they form a gravitationally-bound quadruple system with a 2+2 hierarchy. The two binary stars have orbital periods of P A = 13.62 days and P B = 10.39 days, respectively, and slightly eccentric orbits (e A = 0.07, e B = 0.03). The four stars are very similar in terms of mass (individual masses between 1.105M and 1.195M ), size (radii between 1.09R and 1.26R ), and effective temperature (between 6188 K and 6434 K). The quadruple has an orbital period of P AB = 432 days and an eccentricity of e AB = 0.32. The entire system is essentially coplanar -the orbits of the two binaries and the wide orbit are aligned to within a fraction of a degree. TIC 454140642 is the newest addition to the small family of confirmed, well-characterized, multiply-eclipsing quadruple systems. . Photodynamical model (red line) plotted with the TESS data (left, blue) and WASP data (right, blue). Dark blue dots represent the data points used for fitting the photodynamical model. In the case of TESS data (left) these are the observed data points within the ±0. p 03 phase-domain of each eclipse, while in the case of the WASP observations (right) these are the 1800-sec averages of the individual data points. The pale blue points represent the remaining, out-of-eclipse TESS data (in the left panel), and the original WASP data (to the right). Residual curves are also plotted in the lower panels. Figure 12. Illustration of the various periods discussed in the text for binaries A and B (upper and lower panels, respectively). The grey curves represent the (instantaneous), osculating anomalistic periods (as derived from the osculating true anomalies via Kepler's third law). The brown circles represent the corresponding averages for one sidereal period of the binary (i.e., the interval during which the true longitude changes by 2π). The dark-green squares represent the osculating anomalistic period at the time of each periastron passage of the corresponding binary, while magenta triangles denote the osculating anomalistic period at the time of the primary eclipses. The red/blue curves show the elapsed time between two primary/secondary eclipses, and can be considered to be the true, cycle-by-cycle sidereal (eclipsing) periods between two consecutive primary/secondary eclipses. Finally, the horizontal black lines represent the "average" sidereal (eclipsing) periods -i.e., the periods which can be obtained from the corresponding linear term in the long-term ETV fitting; these are the periods listed in Fig. 9.  Figure 13. The folded RV curve of binary A, after the removal of the contribution of the outer orbit. Red and orange lines represent the fold of primary and secondary RV curves for the forthcoming 9.5 years. The consequence of the rapid, dynamically forced apsidal motion is clearly manifested in the varying shape of the curves. Black points represent the measured RVs (see in Table 3), while the blue lines represent the folded model curves during the interval of the observations. As one can see, the half-year duration of the available observations is insufficient to detect the apsidal motion, but it will be certainly detectable within a few years. Figure 15. Upper panels: Evolution of the semimajor axis (left), eccentricity (middle), and inclination (right) of the two binaries and the quadruple over the course of 100,000 days (≈ 230 outer periods). Lower panels: Same as above but zoomed-in to show the details around the 435-day outer orbit. We note that the periodic extrema in the semimajor axes and eccentricities correspond to the periastron passages of the outer orbit.  Evolution of the impact parameters between each star in each binary with respect to both stars of the other binary (e.g. AaBa for the impact parameter between Aa and Ba, AaBb for the impact parameter between Aa and Bb, etc.) over the course of 5 millions days days. Lower panel: same as above but zoomed-in on a 300-days interval to highlight two conjunctions.