BG Ind: the nearest doubly eclipsing, compact hierarchical quadruple system

BG Ind is a well studied, bright, nearby binary consisting of a pair of F stars in a 1.46-day orbit. We have discovered in the TESS lightcurve for TIC 229804573 (aka BG Ind) a second eclipsing binary in the system with a 0.53-day. Our subsequent analyses of the recent TESS and archival ground-based photometric and radial velocity data, reveal that the two binaries are gravitationally bound in a 721-day period, moderately eccentric orbit. We present the results of a joint spectro-photodynamical analysis of the eclipse timing variation curves of both binaries based on TESS and ground-based archival data, the TESS lightcurve, archival radial velocity data and the spectral energy distribution, coupled with the use of PARSEC stellar isochrones. We confirm prior studies of BG Ind which found that the brighter binary A consists of slightly evolved F-type stars with refined masses of 1.32 and 1.43 $M_\odot$, and radii of 1.59 and 2.34 $R_\odot$. The previously unknown binary B has two less massive stars of 0.69 and 0.64 $M_\odot$ and radii of 0.64 and 0.61 $R_\odot$. Based on a number of different arguments which we discuss, we conclude that the three orbital planes are likely aligned to within 17$^\circ$.


INTRODUCTION
BG Ind (κ1 Ind; HD 208496; TIC 229804573) is a bright, sixth magnitude eclipsing binary (EB) formed by two F-type stars. Its variability was reported first by Strohmeier et al. (1964), and its EB nature was found by Manfroid & Mathys (1984). At the same time, Andersen et al. (1984) obtained the first two spectrograms, and concluded that BG Ind is also a double-lined spectroscopic binary and calculated stellar masses and radii for the first time. The first photometric lightcurve analysis was carried out by van Hamme & Manfroid (1988).
In the forthcoming decades several new photometric and E-mail: borko@electra.bajaobs.hu spectroscopic observations were carried out. They are nicely summarized in Rozyczka et al. (2011), and therefore we do not repeat them here.
The most recent thorough spectroscopic and photometric analysis was carried out by Rozyczka et al. (2011). These authors analysed all the available lightcurves and radial velocity (RV) data including their own measurements. They performed extensive spectroscopic analyses to obtain accurate stellar temperatures, system abundances and then age and evolutionary status. We will compare their results with our findings later in Sect. 4, and therefore, here we highlight only a few noteworthy details. First, they found that the more massive and larger star has the lower temperature 1 , thereby indicating clearly that this component has already evolved away from the main sequence and is moving toward the subgiant regime. Second, they made attempts to resolve some problems with both the photometric phasing (already first noted in van Hamme & Manfroid 1988) and discrepancies in the systemic γ velocities obtained in the solutions of the RV curves measured during three highly different epochs by Andersen et al. (1984); Bakış et al. (2010); Rozyczka et al. (2011). However, they were not able arrive at any definitive conclusions regarding these inconsistencies.
We further note that BG Ind was included in the catalog of those detached eclipsing binaries for which the constituent masses and radii are known to at least 2% precision (Southworth 2015). And BG Ind was also selected for inclusion in the sample of 156 detached eclipsing binaries, which can be used as benchmarks for trigonometric parallaxes in the Gaia era (Stassun & Torres 2016). Finally, turning to the Gaia era, with the use of Gaia DR2 (Gaia collaboration 2018) and HIPPARCOS (van Leeuwen 2007) data, a significant proper motion anomaly was found that might indicate the presence of further, gravitationally bound components in the system (Brandt 2018;Kervella et al. 2019). At this point it should also be noted, that there is a remarkable discrepancy between the revised HIPPARCOS and Gaia EDR3 (Gaia collaboration 2020) parallaxes of BG Ind (πHIP = 14.90 ± 0.59 mas vs πEDR3 = 19.44±0.52 mas), which might be a further indicator of additional multiplicity in the system.
In this paper we confirm the-at least-quadruple nature of BG Ind. Using the high-precision TESS photometry with 2-min cadence we have discovered an obvious second eclipsing EB in the lightcurve of BG Ind with a period of 0.53 days 2 . Our comprehensive investigation of the TESS photometry, archival ground-based photometry and radial velocity curves, as well as the eclipse timing variations (ETV) data demonstrates that the two eclipsing binaries form a close 2+2 quadruple stellar system with a remarkably short outer period of ∼ 2 years.
In Sect. 2 we describe all the available observational data and their preparation for the complex, joint photodynamical analysis which is discussed in Sect. 3. Then, the results are discussed and finally, summarized in Sects. 4 and 5.
2 OBSERVATIONAL DATA

Catalog data
In Table 1, in addition to other catalog data, we collected the photometric passband magnitudes of the system from different surveys, e.g. Tycho-2 (Høg et al. 2000), 2MASS 1 In most binaries with unevolved stars the primary star is the more massive and the hotter star. In the A binary of this system (the dominant binary), the more massive star turns out to be the cooler of the two due to its evolution. We continue to refer to the more massive star as the 'primary' even though it is cooler. However, we will still refer to the deeper eclipse (i.e., when the primary eclipses the less massive (but hotter star) as the 'primary eclipse'. The usual naming convention still holds for the fainter binary B. 2 This second eclipsing binary was also found independently by Eisner et al. 2021. 51.0 ± 0.5 9 References.
(1) Gaia EDR3 (Gaia collaboration 2020); (2) HIPPARCOS (revised) (van Leeuwen 2007); (3) TIC-8 catalog (Stassun et al. 2018); (4) Tycho-2 catalog (Høg et al. 2000); (5) 2MASS All-Sky Catalog of Point Sources (Skrutskie et al. 2006); (6) AllWISE catalog (Cutri et al. 2013); (7) GALEX-DR5 (GR5) (Bianchi et al. 2011 (Skrutskie et al. 2006), AllWISE (Cutri et al. 2013), GALEX (Bianchi et al. 2011) and Gaia (Gaia collaboration 2020). These will be used to construct the spectral energy distribution ('SED') of the system. In turn, the SED along with theoretical isochrones and the photodynamical model of the system provide an opportunity to determine the masses of the components in an astrophysical model-dependent way (see Sect. 3 for details). Together with the passband magnitudes given in Table 1, we list their uncertainties as tabulated in the given catalogs. For the SED analysis, however, we used a minimum uncertainty of 0.03 mag to avoid the overdominance of the extremely precise Gaia magnitudes and also to counterbalance the uncertainties inherent in our interpolataion method during the calculations of theoretical passband magnitudes that are part of the fitting process. Furthermore, similar to the approach followed by Stassun & Torres (2016) we omitted the GALEX near-UV magnitude from our analysis as a distinct outlier. K. Stassun (private communication) kindly called our attention to the fact that even the largest available NUV aperture is missing flux.

TESS photometry
The TESS space telescope (Ricker et al. 2015) Figure 1. Two four-day sections from the beginning of Sector 1 and the end of Sector 28 SAP lightcurves of BG Ind (blue circles). The red and grey curves are spectro-photodynamical model solutions (see later, in Sect. 3). In the case of the red solution the small extra fluctuations of the lightcurve are probably due to the chromospheric/photospheric activities of the stars and were modelled mathematically with Fourier-harmonics simultaneously with the two-binary model, while the grey curve represent the pure twobinary part of the same solution. The residuals to the models are also shown below the lightcurves. SAP) lightcurves from the MAST portal 3 . We used the SAP lightcurves for our study. Because the presence of the small extra dips belonging to the eclipses of the previously unknown binary B (see Fig. 1) was discovered shortly after the release of the data of the first four TESS sectors, our analyses were carried out mostly with the use of Sector 1 data. We also did use Sector 27 and 28 data, but mainly for the purpose of extending the interval of the ETV study. Since, in the case of the faint binary B, the only sources of ETV data are the three sectors of high-quality TESS data, the inclusion of these new observations into our analysis significantly improved the accuracy of the outer orbit solution (including the dynamically inferred mass of binary B).

WASP photometry
BG Ind is one of millions of stars that have been observed as part of the WASP survey. The survey is described in Pollacco et al. (2006) andCollier Cameron et al (2006 July 2012 the WASP-South instrument was operated using 85-mm, f/1.2 lenses and an r filter. With these lenses the image scale is 33 arcsec/pixel. Observations of BG Ind were obtained simultaneously in two cameras on WASP-South over three observing seasons, from 2012 July 3 to 2014 December 6. Fluxes are measured in an aperture with a radius of 132 arcsec for the 85-mm data and instrumental trends are removed using the SYSRem algorithm (Tamuz et al. 2005). Data points more than 5 standard deviations from a phasebinned version of the light curve were rejected and the entire night of data was rejected if more than 1/4 of the observations were identified as outliers based on this criterion.
An eight-day section of the WASP measurements is shown in Fig. 2. Note, we converted the original HJD(UTC) times of the WASP observations to BJD(TDT) for the forthcoming analyses. We also applied the same transformations for all the archival data that we describe below.

Other ground-based archive photometric data used for our analysis
We downloaded publicly available Strömgren u, v, b, y photometric observations from the ESO archive (Manfroid et al. 1991;Sterken et al. 1993). These observations were carried out between JDs 2 446 581 and 2 447 069. The lightcurves in each bandpass contain 175 measurements. We used these data primarily to determine additional times of eclipses in binary A. Unfortunately, however, the majority of the nightly observations contain only a few measurements, and we were therefore able to determine the mid-times of only two primary eclipses (see below, in Sect. 2.6) from this data set. BG Ind was also observed by Jens Viggo Clausen and collaborators as part of their long-running observing programme to measure absolute dimensions for solar-type stars in eclipsing binaries, carried out since 1994 at the Strömgren Automatic Telescope at ESO, La Silla. Unpublished data and manuscripts that were in preparation from this observing programme have been made available to one of us (PM), from which we extracted 5 more eclipse times (see again, in Sect. 2.6). The observing procedures and data reduction for these observations are similar to those described in Clausen et al. (2001).

Disentanglement of the lightcurves
For the combined analysis of all the observational data we used the original TESS timeseries, i.e., the net lightcurve of the two binaries together. However, at the start of the analysis, we found it worthwhile to disentangle the lightcurves of the two binaries, so we could examine each one separately. We have described this process in substantial detail in Powell et al. (2021). Therefore we review only the highlights here. First we folded and binned the Sector 1 (i.e. Year 1) TESS SAP lightcurve with the period of binary A into 1000 equal phase cells. However, while producing the fold for binary A, we excluded those data points which were recorded during the eclipses of binary B. Then the mean flux of each cell was rendered to the mid-phase value of that cell. In such a way we obtained a folded, binned, and averaged lightcurve of binary A (see upper panel of Fig. 3). Then this lightcurve was removed from the original sector 1 SAP lightcurve point by point in such a manner that the flux to be removed at the actual phase of any given data point was calculated with a three-point local Lagrange-interpolation from the folded, binned, and averaged lightcurve. As the result of this removal, we have obtained a new, residual time series which now mainly contains the light variations of binary B, 4 without the eclipses and ellipsoidal variations of binary A. Therefore, this lightcurve can be used for determining the mid-eclipse times of binary B.
In the next step we folded, binned, and averaged this residual lightcurve with the period of binary B (see middle panel of Fig. 3). Finally, we subtracted this folded, disentangled lightcurve of binary B from the original Sector 1 TESS SAP lightcurve, thereby obtaining a time series of binary A without the small distortions caused by binary B. We applied the same process to the Sector 27 and 28 (Year 3) SAP lightcurves, as well (see bottom panel of Fig. 3).
Regarding the WASP observations, we carried out a very similar process with the slight modification that, in this case, for the much smaller number of individual data points we applied binnings of 200 and 500 cells instead of 1000. We carried out the whole process separately for the three seasons of the WASP observations. Though the eclipsing signal of the faint binary B is not readily detected in the original WASP timeseries, we were able to see it clearly in our disentangled version (see Fig. 4).
Finally in regard to disentangling the lightcurves, we have also used a second method which fits simultaneously for 50 harmonics of each of binaries A and B given their established periods. This technique, which is also described in detail in Powell et al. (2021), involves inverting a 201 × 201 matrix to solve for the linear coefficients to the 50 sines and 50 cosines for each of the two binaries. We find nearly perfect agreement for the disentangled TESS lightcurves from the two independent methods, and thus we do not show those results here. 4 Note, for practical reasons, we added a constant flux to these time series in such a way that the flux of the very first data point retained the same value as in the original time series. In this manner, we replaced the varying light of the extracted binary with a constant extra light.  Figure 3. Folded, binned, and averaged TESS lightcurves of the two binaries of the quadruple system BG Ind. Upper panel: Sector 1 lightcurve of binary A (blue circles), together with the folded, binned, and averaged combined spectro-photodynamical model lightcurve (red curve, see later, in Sect. 3). Middle and lower panels: Year 1 (Sector 1) and Year 3 (Sector 27 and 28) lightcurves of binary B, respectively. As in the case of the binary A lightcurve in Fig. 1, the red solution curve exhibits some small extra fluctuations that are probably due to the chromospheric/photospheric activities of the stars (see Sect. 3 for details). These were modelled mathematically with Fourier-harmonics simultaneously with the two-binary model, while the thin grey curves represent the pure two-binary part of the same solution. The fold of the residuals to the models are also shown below the folded lightcurves.  Figure 4. Folded, binned, averaged WASP lightcurves of binary A and B for observing season 2012/2013. For illustrative purposes we phased both curves with the ephemeris calculated for Sector 1 TESS data. In such a manner, the shift of the primary and secondary eclipses from phases 0. p 0 and 0. p 5, respectively, i.e., the phasing problem mentioned in the Introduction, is clearly visible. In the case of disentangling the WASP data, the results for binary B are actually somewhat improved using the Fourier approach and we show that lightcurve in Fig. 5 as well for comparison.

TESS ETV results
In order to calculate accurate eclipse times from the TESS lightcurve we used the disentangled time series (see above in Sect. 2.5). The 91 eclipse times of binary A (including both primary and secondary eclipses) from Sectors 1, 27, and 28 are presented in Table 2. In Table 3 we list the eclipse times for binary B including a combined 259 primary and secondary eclipses.

Ground-based ETV results
We also utilized WASP and ESO (Manfroid et al. 1991;Sterken et al. 1993) data, including the unpublished observations of J. V. Clausen, as well, to calculate 85 additional eclipse times for binary A. Furthermore, a primary and a secondary eclipse of BG Ind were observed by one of us (MB) using a DSLR camera. Images were recorded in RAW format and the green, blue and red channels were extracted into separate images. The times of minimum were measured from each colour filter with the Peranso software 5 using a 5th order polynomial fit. The average of the mid-eclipse times were converted into BJD. Finally, we collected one other eclipse time from the the paper of van Hamme & Manfroid (1988) and converted it into BJD. All these eclipse times are tabulated in Table 4. The eclipses from binary B in the archival data were too weak to derive meaningful eclipse times.

BG Ind ETV Results
The overall ETV curves for BG Ind A and B are plotted in Fig. 6 along with the best-fit spectro-photodynamical model that is described in the next section. The ETV curve of BG Ind A exhibits a clearly cyclic pattern with a period of ∼ 2 years. Even in the absence of any other indications of additional stars in the system, the most plausible explanation of this ETV behavior would be the light-travel time effect (LTTE) caused by a gravitationally bound, distant, third component. Therefore, we carried out, a preliminary, 'traditional' analysis of the ETV curves of binary A by fitting the LTTE-term with our analytic ETV-solver (Borkovits et al. 2015). We found that the very first eight ETV points deviate systematically from the LTTE solution. Therefore, we added a quadratic term to the analysis and obtained the following quadratic ephemeris: We also tabulate the parameters of this preliminary LTTE solution in Table 5, and plot this simple model together with the spectro-photodynamical model, in Fig. 6. Turning to the ETV points of binary B, they appear to be moving with the opposite phase to that of the bright binary A which makes it very likely that the two binaries form a bound, quite tight quadruple system. As we will discuss below in Sects. 3 and 4 our detailed analysis robustly confirms this hypothesis.

Radial Velocity data
We used three sets of radial velocity (RV) data for our analysis. These are as follows: (i) Bakış et al. (2010)   Notes. Integer and half-integer cycle numbers, as above, refer to primary and secondary eclipses, respectively. Left panel shows all the available observations, while in the right panel we zoom in on the regions of the better-covered WASP and TESS data. Larger red circles represent ETV points calculated from the observed eclipse events of binary A, while the smaller blue circles stand for the ETVs of binary B. Note, for simplicity, we do not separate primary and secondary eclipses. (The validity of this can easily be verified since both binaries have circular orbits and, furthermore, the primary and secondary eclipses within each binary can be calculated with the same accuracy due to their similar depths.) Black and grey lines stand for the combined spectro-photodynamical model ETV solution (Sect. 3) for binary A and B, respectively, while the green line denotes the preliminary, 'classic', analytic LTTE+quadratic ETV solution discussed in Sect. 2.6.3. (Note, for clarity, in the left panel, the ETV solution of binary B, i.e., the grey curve, is plotted only for the narrow interval around the TESS observations.) The residuals of the observed vs. modelled ETVs are plotted in the bottom panel. Here, as above, red and blue dots represent the residuals of binary A and B ETV points against the spectro-photodynamical model, while green dots stand for the residuals of binary A data against the analytic LTTE+quadratic ETV model.  Rozyczka et al. (2011) with the fibre-fed Giraffe spectrograph on the 1.9-m Radcliffe telescope at the South African Astronomical Observatory; and (iii) finally, we found in the ESO publicly accessible archive 6 a number of spectra taken with the HARPS and FEROS spectrographs between JDs 2 453 191 and 2 456 910. From these data we determined an additional 54 RV points of both components of binary A.
To measure the RVs from these latter spectra we used the 6 http://archive.eso.org/cms.html broadening function method (Rucinski 1992) as implemented in the software package RaveSpan (Pilecki et al. 2017). The template used for the analysis was a synthetic spectrum for a star with T eff = 7000 K, log g = 3.5. We used a simultaneous least-squares fit of two rotationally-broadened profiles to measure the radial velocities of the two stars from the broadening profile. The RV points obtained in this way are tabulated in Table 6.
The phase-folded RV points (after the correcting for the orbital motion around the center of mass of the whole quadruple system) together with the best-fit photodynamical solution (see below, in Sect. 3) are plotted in Fig. 7.

JOINT ANALYSIS OF THE AVAILABLE DATA
We used the software package Lightcurvefactory (see Borkovits et al. 2019aBorkovits et al. , 2020, and further references therein) to carry out a complex spectro-photodynamical modeling of the system based on the data collected in Sect. 2. Lightcurvefactory calculates stellar positions and velocities for each object and emulates the light-, RV and ETV curves of any arbitrary quadruple system (having either 2+2 or 2+1+1 hierarchies), including mutual eclipses amongst any two (or more) components. Moreover, the software may (optionally) use built-in, pre-calculated PARSEC isochrone tables 7 (Bressan et al. 2012) to constrain the stellar parameters theoretically Our combined analysis is primarily based on the following observational inputs: (i) the high-quality TESS lightcurve (see Sect. 2.2), (ii) the RV data available in the literature (see Sect. 2.7), (iii) the ETV data calculated from all the available photometric observations (see Sect. 2.6 and Tables 2-4) and, (iv) the observed passband magnitudes of the target taken from standard catalogs (see Sect. 2.1 and Table 1).
Note that, in its present form, Lightcurvefactory is unable to handle period variations caused by non-few-body perturbations. Therefore, we did not model the small linear period variation of binary A which manifests itself in the form of    11.5 ± 0.6 P A /P A [10 −8 yr −1 ] 6.1 ± 0.4 Notes. a AB sin iout denotes the line-of-sight projected semi-major axis of the outer orbit of binary A around the center of mass of the quadruple system, while the other orbital elements and associated parameters are noted in their usual manner. Moreover, we tabulate two derived parameters: f (m B ), and K A which are the mass function and the amplitude of the radial velocity curve of the center of mass of binary A on its outer orbit. Finally, in the last row, we give the rate of the continuous period variation of binary A, which is derived from eqn. (1).
quadratic deviations of the first few ETV points (see above, in Sect. 2.6.3). Considering the fact that our analysis depends primarily on the TESS measurements obtained over the last 2.5 years, and on the radial velocity data gathered within a relatively narrow eight-year-long interval about a decade ago, we do not expect that such a small, long-term effect will have any significant influence on our results. Nevertheless, we will return to this question in Sect. 4. Regarding the archival photometric observations, we decided not to use the photometric fits to these lightcurves themselves for the analysis. This decision was based mainly on the fact that their large scatter was found to be of the same order as, or even higher than, the eclipse depths of the faint binary B. We did, however, utilize in our analysis the most relevant information that could be mined from these observations, namely the best of the mid-eclipse times that could be derived from these data. Furthermore, the other benefit of these data is that they reveal that the eclipse depths of binary A have remained constant during the last ∼ 40 years, the relevance of which will be discussed later.
Before the analysis, we also took further preparatory steps on the TESS lightcurve. In order to save computational time we binned the 2-min cadence data into 30-min bins, and in the following analyses we worked with the binned data. While carrying out some preliminary fitting runs on this 30-min binned Sector 1 lightcurve, we realized that the residual curve exhibits small amplitude, quasi-periodic variations. We had also found very similar patterns in the residual lightcurve at the end of the prior lightcurve disentangling process (Sect. 2.5), i.e., after removing both binaries from the original time series. We therefore concluded that these small quasi-cyclic variations cannot be the consequence of some misadjustements of the lightcurve parameters during our analysis, but should be real effects. In order to find the dominant frequencies of these fluctuations we calculated the power spectrum of the residual curve with a discrete Fourier-transform. We found two independent sets of frequency peaks. The frequencies of one set were close to the orbital frequency of binary A (and its multiples), while the other set was clearly related to the orbital frequency of binary B.
Similar to what we have done in some of our previous work (see e.g. Borkovits et al. 2018) we modelled these fluctuations during our analysis in the following manner. We found that the use of the two most dominant frequencies of both sets of frequencies resulted in a significant improvement in the solutions. The process itself works as follows: in each trial step, after the removal of the blended EB lightcurves from the observed data, the residual lightcurve is modelled with harmonic functions of the four fixed frequencies, of which the eight (plus one) coefficients are obtained via matrix inversion. Then, this mathematical model of the residual lightcurve is added to the double binary model lightcurve and the actual χ 2 value is calculated for this mixed model lightcurve. Finally note, since the fluctuations were found to be quasi-periodic instead of strictly periodic, we found that our process is the most effective if we use only a short section of the TESS lightcurve. Therefore, for the main portion of our analysis we used only a seven-day-long section of the Sector 1 TESS lightcurve. 8 The combined analyses were carried out in two different stages. In the first stage we worked only with astrophysical model-independent parameters. Therefore, we fitted simultaneously only the TESS lightcurve, and the RV and ETV curves, but did not include SED data and theoretical stellar isochrones. During this phase the 20 adjusted parameters were as follows: (i) Seven lightcurve related parameters: the temperature ratios of (T2/T1)A,B and TBa/TAa, i.e. the secondary over primary temperature ratios of both binaries, and the ratio of the temperatures of the two primaries; the durations of the two primary eclipses (∆tpri)A,B; the ratios of the radii in both pairs (R2/R1)A,B; and the gravity darkening coefficients of the two stars of binary A (β Aa,Ab ). (ii) One parameter for each inner binary orbit, i.e. the observed inclinations iA,B of the orbital planes of binary A and binary B, and five orbital parameters of the outer orbit: period (Pout), time of periastron passage τout, eccentricity and argument of periastron (e cos ω)out and (e sin ω)out, and the inclination iout.
(iii) Four mass-related parameters: the masses of the two primaries (mAa,Ba), and the mass ratios of the two binaries (qA,B).
Regarding the other orbital parameters of the inner binaries, the periods (PA,B) of these EBs, as well as their orbital phase (through the time of an arbitrary primary eclipse -T pri A,B ) were constrained internally through the ETV data. Furthermore, the eccentricities of both inner orbits were set to zero. Moreover, for the large Pout/PA,B ratios we found that all three orbits (two inner binary orbits and the outer orbit) can be considered as pure, unperturbed Keplerian motion. Due to this latter consideration our dataset does not contain any information about the positions of the orbital nodes relative to each other. Therefore, the sixth orbital element, the longitude of the node of each orbit (ΩA,B;out) was fixed at zero. Finally, we note that the systemic radial velocity of the whole quadruple system (γ) was also constrained internally by minimizing the χ 2 RV contribution a posteriori in each trial step.
Turning to the atmospheric parameters of the four stars, in contrast to our previous analyses, we now adjust the gravity darkening coefficients (β) of the strongly non-spheroidal components of the bright binary A. The reason is that, in contrast to the widely used classic model of Lucy (1967) which predicts a unique gravity darkening coefficient of β = 0.32 for all convective stars, recently Claret & Bloemen (2011) have shown that the true relations are much more complicated. This is especially true for stars close to the transition region between convective and radiative envelopes, where the components of binary A are located. On the other hand, in the case of binary B, we kept fixed the usual value of β = 0.32 prescribed in Lucy's model. Other atmospheric parameters, such as the logarithmic limb-darkening coefficients (x, y)TESS were interpolated in each trial step with the use of passbanddependent tables downloaded from the Phoebe 1.0 Legacy page 9 . These tables are based on Castelli & Kurucz (2004) atmospheric models and are primarily used for the original version of the Phoebe software (Prša & Zwitter 2005). Furthermore, for the components of the bright binary A (Aa and Ab), we include the reflection/irradiation effect into the lightcurve model and, therefore, we take into account the bolometric limb-darkening coefficients (x, y) bol , interpolating them in each trial step in the same manner as was done with the passband-dependent coefficients.
At this stage of the analysis we required only one further parameter that is undetermined by the model and, therefore, has to be set externally. This was the effective temperature of the primary of binary A, which was set (and kept fixed) according to the findings of Rozyczka et al. (2011).
At the end of this stage of the analysis we obtained accurate dynamical masses not only for the two members of the bright binary A but, in addition, we obtained the total dynamical mass of binary B. 10 Furthermore, the temperature ratio of the two primaries provides reliable information about the characteristics of the two stars forming the faint binary B. Finally at this stage, the orbital elements of the three orbits were also accurately determined.
In the next and final stage of the analysis, we included the SED information into the analysis as well as the built-in PARSEC tables. Now the seven lightcurve-related parameters described above were no longer adjusted but, instead, the radii and temperatures of all the four stars were constrained, i.e., recalculated at the beginning of each trial step by interpolating their values from the three-dimensional (mass, metallicity, age) grids of the PARSEC tables 11 . During this phase of the analysis three additional adjustable quantities were introduced, including (i) the metallicity ([M/H]) and (ii) the (logarithmic) age of the quadruple. These two parameters, together with the mass of the given components, determined the position of each star within the PARSEC grids and, therefore, determine the interpolated fundamental stellar parameters and theoretical passband magnitudes. The third parameter (iii) was the stellar extinction E(B − V ) for the SED fitting. Moreover, while fitting the model SED to the dereddened observed SED points, the distance of the system comes in as an additional parameter. The software constrains this parameter a posteriori in each trial step by minimizing the value of χ 2 SED . After some initial trials, however, we found it necessary to introduce a fourth, extra parameter to adjust, namely the age of the evolved component of binary A, in order to obtain model lightcurves which yield similarly low χ 2 LC values to the ones obtained in the previous astrophysical modelindependent stage. This procedure requires some further explanation. It is generally expected that the components of a close binary (multiple) system are coeval. Theories, however, allow for small departures from exact coevality (see, e.g. Tokovinin 2018b), which during some critical rapid stages of stellar evolution might be significant. Furthermore, even in the case of exact coevality, the approximative nature of our interpolation method certainly carries with it inherent inaccuracies which might lead to modest discrepancies in the derived stellar parameters, especially during the very rapid sensitive evolutionary stage of the evolved star in binary A. Therefore, as a counterbalance to these uncertainties, we allowed for the age of the evolved component to be set independently from the other three stars.
In Table 7 we tabulate the median values and the 1σ statistical uncertainties of the parameters obtained during the last stage of our analysis. The synthetic model light curves derived from the best-fit joint solution are displayed in Figs. 1 and 2. The corresponding ETV and RV curves are presented in Figs. 6 and 7, respectively. Finally, in the two panels of Fig. 8, we illustrate the SED-fitting part of the combined solution both in the flux and the passband magnitude domain.

Orbital configuration
Our analysis confirms the hierarchical 2+2 type quadruple star nature of BG Ind. Thanks to the available high-quality TESS photometry and the long-term ground-based photometric and spectroscopic observations, BG Ind now takes its place as (i) one of the most compact 2+2 quadruples known, as well as (ii) the quadruple system with the most accurately Notes. a: Time of the inferior conjunction of the secondary component (i.e. mid-time of a primary eeclipse); b: Interpolated (or derived) from the PARSEC isochrones; c: interpolated linear (x) and logarithmic (y) limb-darkening coefficients. Note, bolometric coefficients used only during the calculation of the reflection effect, therefore they were not set for binary B; d: gravity darkening coefficients; e: the age of the evolved primary component of binary A was allowed to vary independently of the other three stars -see text for details.
known stellar masses and other stellar parameters. The outer period of the system is found to be Pout = 721 ± 3 days, which is the shortest amongst doubly eclipsing quadruple systems with an accurately known outer period. 12 Note, how- 12 We emphasize that this holds only for doubly eclipsing 2+2 quadruples. The shortest period known 2+2 quadruple system, VW LMi has a much shorter outer period of Pout = 355 days (Pribulla et al. 2008(Pribulla et al. , 2020. Furthermore, in the case of the doubly eclipsing quadruple star EPIC 220204960 Rappaport et al. (2017) ever, that despite the relatively short outer period, both outer to inner period ratios are large enough (Pout/PA ≈ 492, Pout/PB ≈ 1365) so that we do not expect readily measurable short-term mutual three-(four-) body perturbations.
In other words, all three orbits can be considered as essentially purely Keplerian. The outer orbit is moderately eccenreported that the outer period is very likely between 300 and 500 d, but an accurate value for that system is unknown.  Table 1) versus the model passband magnitudes derived from the absolute passband magnitudes interpolated with the use of the PARSEC tables (blue filled circles). In the right panel the dereddened observed magnitudes are converted into the flux domain (red filled circles), and overplotted with the quasi-continuous summed SED for the quadruple star system (thick black line). This SED is computed from the Castelli & Kurucz (2004) ATLAS9 stellar atmospheres models (http://wwwuser.oats.inaf.it/castelli/grids/gridp00k2odfnew/fp00k2tab.html). The separate SEDs of the four stars are also shown with thin green, black and purple lines, respectively.
tric with eout = 0.21 ± 0.05, and is seen nearly along the direction of the minor axis (ωout = 2 • ± 9 • ). These relatively small uncertainties in the BG Ind quadruple system, however, should be treated with some caution. The two main reasons for this caveat can nicely be seen in the ETV plots in Figure 6. First, due to the unlucky fact that the outer period is nearly exactly equal to two years (Pout ≈ 1.973 yr), the annual observing seasons of the target can, and do, miss the most informative two parts of the ETV curve, i.e. its two extrema. Second, as was discussed above, the very first eight pre-WASP ETV points show clear deviations from the pure LTTE solution, and might indicate a continuous, constant increase in the orbital period of binary A, i.e.,ṖA.
Though the modelling ofṖA was not included in the comprehensive spectro-photodynamical approach, its effect can be quantified by comparing the orbital parameters of the outer orbit obtained through the classic, analytic LTTE+quadratic solution of the ETV of binary A (Table 5) with the detailed spectro-photodynamical model ( Table 7). As one can see, the outer period and eccentricity match well within their estimated uncertainties, while the argument of pericenter, the periastron passage time and the RV amplitudes are discrepant at the 2 − 3σ level. Therefore, we can conclude that, as was expected, the omission of the quadratic ETV term in the complex spectro-photodynamical analysis did not influence our basic solutions, but suggests that the actual uncertainties in the orbital elements should be somewhat larger than cited in Table 7.
The inclination of the outer orbit is found to be iout = 86 • ± 5 • . On the other hand, the inclinations of the two inner eclipsing binaries are found to be iA = 73. • 1 ± 0. • 1 and iB = 84. • 3 ± 0. • 9. From these values, and in the absence any information on the longitude of the nodes of the three orbits (ΩA,B,out), the only thing one can say is that the whole quadruple system is certainly not perfectly flat. Since the mutual inclination of two planes cannot be smaller than the difference between the two observed inclinations of the planes considered (and, cannot be larger than their sum), the inclination of the bright binary A relative to the outer orbit must surely exceed ≈ 13 ± 5 • , but may even reach 90 • . (Similarly, the mutual inclination between the orbital plane of binary B and the outer orbital plane may be anywhere between coplanar and perpendicular.) As a consequence of such misalignments, one may expect the binary's orbital plane to precess. In that case, eclipse depth variations should be observed, or even the disappearance of the eclipses on a longer time-scale. The period of forced precession of a binary orbital plane in a hierarchical triple system can be well approximated with the expression (see e.g. Söderhjelm 1975) (Pprec)A = 4 3 1 + qout qout where C represents the total orbital angular momentum of the quadruple, while G2 is the orbital angular momentum stored in the outer orbit. In the present situation it can be readily seen that C/G2 ≈ 1, i.e., the majority of the orbital angular momentum of the quadruple is stored in the outer orbit. Therefore, with qout ≡ MB/MA 0.5 one can easily show that (Pprec)A 4200 yr. From this we conclude that there is no chance of detecting eclipse-depth variations given the available span of the observations.

Astrophysical properties and evolutionary status of the four stars
Turning to the fundamental astrophysical parameters of the bright components of binary A, we compare our results to those of the former thorough analysis of Rozyczka et al. (2011). The masses of the stars found in the two analyses agree quite well: 1.432 ± 0.020 vs. 1.428 ± 0.008 for the more evolved component, and 1.315 ± 0.025 vs. 1.293 ± 0.008 for its less evolved companion, where the first of each pair are from the current work. This agreement is good to 1 σ, in units of our error bars. The uncertainties given in Rozyczka et al. (2011) are smaller than ours by factors of 2-3. One should keep in mind, however, that Rozyczka et al. (2011) estimated their uncertainties from an analysis of their own set of RV data which have smaller rms residuals than the combined set of RVs that we used. Furthermore, during their final RV analysis they corrected the RV values for the distortions of the stellar components with the use of the Wilson-Devinney code (Wilson & Devinney 1971;Wilson 1979). In contrast to this, in our study, the effects of the distortions of the stars on the RV data are automatically taken into account within Lightcurvfactory. Therefore, we consider our somewhat larger uncertainties to be more realistic.
The radii of the two stars in binary A exhibit slightly larger differences between the two studies. Our analysis has yielded RAa = 2.34 ± 0.02 R for the evolved component, while Rozyczka et al. (2011) obtained the somewhat smaller value of 2.29 ± 0.02 R . 13 For the other less evolved star we found R Ab = 1.59 ± 0.04 R , in contrast to their somewhat larger value of 1.68 ± 0.04 R . Keeping in mind, however, that the sum of the (fractional) radii of the two stars (i.e. RAa/aA + R Ab /aA) is one of the most robustly determined parameters of an eclipsing binary's lightcurve (though, it does depend sensitively on the inclination), one can easily check, that this sum agrees well in the two solutions. 14 The relatively small discrepancies in the individual radii between our analysis and that of Rozyczka et al. (2011), can likely be explained by some combination of the following three effects. First, prior studies did not consider the small extra flux contribution (of ≈ 2.2%) coming from binary B, and the lightcurve distortions caused by its eclipses. Second, we used the high quality TESS photometry whose superiority over the former ground-based measurements is unquestionable. Third, we allowed for the reflection/irradiation effect which made our analysis more realistic, but this effect was not considered during the previous analyses. In conclusion, we emphasize again that the discrepancy in radii is fairly small.
We also found small departures in the effective temperatures of the two components of binary A compared with the previous results. Our results, which hinge to a large degree on the fit of the combined 4-star SED, resulted in slightly larger temperatures. We found TAa = 6442 ± 29 K and T Ab = 6816 ± 26 K in contrast to 6350 ± 260 K and 6650 ± 230 K (Rozyczka et al. 2011). Note. however, that our results are within the uncertainties of Rozyczka et al. (2011). On the other hand, by using the temperatures given by Rozyczka et al. (2011), Stassun & Torres (2016 found a consistent 13 Note, however, that they analysed six different lightcurves separately, and their results were scattered between 2.17 R and 2.40 R . 14 Note also that due to the significant tidal and rotational oblateness of the two stars, they will no longer be spherical; therefore, it should be clarified what is meant by 'radius'. We cite the volume equivalent radius and assume that Rozyczka et al. (2011) used the same definition. On the other hand, we note that in the above mentioned relation for the sum of the fractional radii, the so-called 'side' radius, (i.e. measured in the star's equatorial plane, in the direction perpendicular to the line joining the two stars) should be the relevant one during eclipses. The volume equivalent and side radii for our stars, however, agree to better than 1%.
SED solution for the binary. Of course, they did not consider the contribution of binary B, which might give a small excess at the red wing of the SED, and in turn which might force the fit toward slightly lower temperatures.
There is, however, an even more significant discrepancy between the system distance inferred from the SED solution and the trigonometric distance deduced from Gaia's measurements. Our solution has resulted in a photometric distance of d = 69.7 ± 0.8 pc, while Bailer-Jones et al. (2018) using the Gaia DR2 measurements have obtained dDR2 = 51.0 ± 0.5 pc. 15 The situation is more complicated than this seemingly straightforward discrepancy. First, the trigonometric distance that can be calculated from the new reduction of HIPPARCOS parallaxes (van Leeuwen 2007) is dHIP = 67.1 +2.8 −2.2 , which is within 1σ of our result. Furthermore, Stassun & Torres (2016) used a similar SED modeling analysis to ours, and found a photometric distance of dStassun+16 = 66.7 pc. Again, this is much closer to our result and that of HIPPARCOS than to the Gaia distance. However, since we know the fundamental stellar parameters of the dominant A binary quite accurately (including the bolometric luminosities) independent of the distance, and others have used partially different methods 16 to find a very similar distance, we tentatively conclude that the published Gaia DR2 and EDR3 parallaxes are probably subject to some systematic error. This discrepancy might have arisen from the fact that the period of the outer orbit in BG Ind is very close to two years (Pout = 1.973 yr) and, furthermore, the semi-major axis of BG Ind A's ellipse around the center of mass of the quadruple system is aout,A ≈ 0.82 AU. Therefore, the combination of the orbital motion of the photocenter of binary A along the outer orbit and a period near 2 years, may be responsible for causing some problems with Gaia's trigonometric parallax determination.
Turning to the newly discovered binary B, we find it to be a pair of two mid K-type dwarfs with a near unit mass ratio of qB = 0.93 ± 0.01. Despite the fact that the flux contribution from this binary is only about 2% in the TESS photometric band, due to the high quality TESS photometry of this rather bright quadruple, we were able to obtain quite good lightcurves (see Fig. 3) and a robust dynamical model for binary B. The out-of-eclipse sections of the lightcurve of binary B are distorted, and we explain that by chromospheric activity, i.e. stellar spots, which are quite usual for stars with thick convective envelopes. As one can see in the middle panel of Fig. 3, for Sector 1 data these variations can be well modelled (mathematically) with two Fourier-terms having frequencies equal to the orbital frequency and its first harmonic. On the other hand, in the case of Year 3 (i.e. Sectors 27 and 28 data) the same Fourier-representation was found to be less satisfac- 15 The parallaxes published in Gaia DR2 and EDR3 are well within 1σ of each other, and therefore we can assume that the distance derived from the EDR3 data will not differ significantly from the published DR2 distance. 16 Rozyczka et al. (2011) have determined the temperatures with a combination of spectroscopic analysis and lightcurve fitting. Stassun & Torres (2016) utilized SED modelling. And our study included a combination of lightcurve and SED fitting, where the consistency of the obtained stellar parameters were also probed by modelling the RV curves and identifying appropriate PARSEC isochrones tracks. tory, and we therefore added two additional harmonics of the binary B orbital frequency to the Fourier-representation (see the lowest panel of Fig. 3). However, even in this case, one can still notice some imperfections in the fit. We explain this fact with a possible rapid variation in the chromospheric activity which induces brightness fluctuations that cannot be well represented by a few smooth harmonics (even after averaging the lightcurve over a few weeks of the TESS observations). Note, however, that this discrepancy is less than ∼100 ppm, and therefore, it would remain under the detection limit for any ground-based photometric observations. The timings of the shallow eclipses from binary B are in accord with the 1.973-yr-periodicity in the ETV curve found from binary A, both of which are dominated by the light travel time effect (Fig. 6). Moreover, the dynamically deduced total mass of binary B, coupled with the dimensionless stellar parameters, lead to physical parameters of the stars in binary B. And these, according to the PARSEC tables we used, are fully consistent with the parameters of two mainsequence K-dwarfs, having the same age and metallicity as those of the members of the bright binary A. Therefore, there is no question that the two eclipsing binaries form a compact, gravitationally bound, hierarchical quadruple star system.
Regarding the global parameters of the quadruple, the combined solution prefers a slightly metal-deficient abundance of [M/H] = −0.19 ± 0.04 which, again, is in perfect agreement with the previous result of [F e/H] = −0.2 ± 0.1 by Rozyczka et al. (2011). As mentioned before, we did not enforce strict coevality among the stellar components during our analyses and found an age of τAa = 2.51 ± 0.12 Gyr for the evolved primary, and τ Ab,B = 2.14 ± 0.10 Gyr for the three main-sequence components. These two ages differ by 370 ± 150 Myr, or ∼7% of the age, with a significance of only 2.5 σ. We consider this discrepancy to be not a 'small' departure from the coevality (though it does not have a high statistical significance). Our impression is that it might arise from the rapid rotation as well as the tidal distortion of the evolved star. Therefore, it is possible that the spherical stellar radius given by the PARSEC tables would not strictly equal the volume-equivalent radius of a strongly tidally distorted star.
Finally, we briefly discuss the question of the likely continuous period increase in the binary A period, detected through the systematic deviations of the very first ETV points from a simple linearly sloped LTTE model. Such period variations have been observed in a fair number of EBs. In the case of semi-detached and contact systems, the most common explanation is some kind of mass exchange between the stellar components. However, since all the previous studies have found that BG Ind A is a detached system, and our detailed analysis confirms this scenario, this period increase cannot be explained via mass transfer. (And this is not to mention the fact that, in this case, the increasing period would imply that the less evolved, lower mass secondary star would be the mass-donor, which is an unphysical scenario.) On the other hand, however, as was shown, e.g. by Pringle (1985) and Demircan et al. (2006) mass loss from a close binary star (e.g., due to a stellar wind) always leads to an increasing period. Therefore, the quadratic ETV-term in BG Ind A might imply an enhanced stellar wind from the surface of the evolved component. For a quantitative study of this possibil-ity, further observations of high-quality eclipse times over a longer time interval would be extremely useful.
We also note that naturally, an LTTE effect driven by a more distant, low-mass fifth stellar component may also be the correct explanation. Obviously, the confirmation or refutation of this scenario also requires further eclipse follow up observations.

SUMMARY AND CONCLUSIONS
In this paper we report the discovery of the doubly eclipsing quadruple nature of the previously known, bright, southern eclipsing binary BG Ind. We present the first comprehensive analysis of BG Ind in its entirety. TESS observations provided high-precision photometry covering two intervals of one and two months, respectively, and separated by two years. Even though these high-quality TESS observations covered only short segments of the outer orbit with Pout = 1.973 yr, we were able to use ground-based archival lightcurve and RV data to determine accurately the orbital and dynamical parameters of the system. BG Ind is found to have one of the shortest outer periods among all quadruple systems having a 2+2 hierarchy. According to the recent version of the Multiple Stellar Catalog (Tokovinin 2018a) there are only five such systems (including BG Ind) with outer periods shorter than 3 years. These are tabulated in Table 8.
The remarkably small number of such compact 2+2 quadruples that are known probably arises from an observational selection effect rather than for astrophysical reasonsspecifically, they are quite difficult to discover. In contrast with the discovery of a third companion of a known binary star, which can be made by, e.g., astrometry, long-term RV measurements (or, in the case of an EB) ETV studies, or in exceptional cases, observing serendipitous extra eclipses, the binary nature of such a third component would remain hidden in most cases. The only rare exceptions are when the second binary happens to be also an EB 17 (as is the case in four of the five systems listed in Table 8), or it is bright enough to be observable as a second spectroscopic binary (as in the case of the fifth tabulated system, VW LMi). Furthermore, another possibility for discovering the binarity of a component in a binary or multiple star system might arise from its faintness in comparison to its dynamically determined mass (as in the cases of κ For, Tokovinin 2013 and ζ Cnc C, Tokovinin 2017). 18 17 On the other hand, despite the fact that such large photometric surveys such as, e.g., the gound-based Optical Gravitational Lensing Experiment (OGLE; Udalski et al. 2015), or the TESS mission have observed hundreds of lightcurves exhibiting blends of at least two EBs, the gravitationally bound nature of the blended EBs have been proven definitively for only a relatively small fraction of these objects (see, e. g. Zasche et al. 2019). 18 In the case of the compact hierarchical triple system IU Aur Drechsel et al. (1994) and Özdemir et al. (2003) have also concluded that the large third mass vs. small third light and weak spectroscopic signal discrepancies could be resolved by postulating that the third companion is a binary itself. If this assumption was true, IU Aur would be the shortest outer period 2+2 quadruple with Pout = 294 d, but the system needs further investigations. Moreover, note that most recently Marcadon et al. (2020) have Notes. * The 7.931-day-period binary in VW LMi does not exhibit eclipses. References: (1) Pribulla et al. (2008); (2) Rappaport et al. (2017); (3) present paper; (4) Rowden et al. (2020); (5) Zasche & Uhlař (2016) BG Ind is also one of only a very few compact quadruple systems where the key parameters of all fours stars are known with an accuracy of better than ∼3%, including masses, radii, and T eff values. Likewise, the three sets of orbital parameters, such as periods, semi-major axes, eccentricities, and inclination angles are all known rather precisely. The one notable exception is that we do not have any way of determining the mutual inclination angles among the three orbital planes. The observational inclination angles are all close to edge on (i.e., 90 • ), and we might surmise from statistical arguments that the most likely configuration is nearly coplanar for all three orbits. But we cannot be certain that this is indeed the case.
Future interferometric and astrometric observations may help to solve this problem. The semimajor axis of the outer orbit is 36 mas, so in principle it is resolvable by speckle interferometry or adaptive optics. However, the high contrast between the A and B binaries makes resolution with singledish telescopes a challenging task. A much better prospect is offered by long-baseline interferometers, e.g. the GRAVITY instrument at VLTI 19 (Gravity collaboration 2017). The contrast in the K band is more favorable compared to the visible (∆K = 2.98 mag vs ∆V = 4.73 mag), so the visibility and phase modulation caused by the outer pair can be well measured. Furthermore, the spectral resolution of R ∼ 4000 offered by GRAVITY will allow detection of opposite phase shifts in the spectral lines of Aa and Ab at times near their maximum separation, thus enabling one to measure the orientation of the Aa,Ab orbit (its semimajor axis is 0.5 mas) and, perhaps, even the orbit of Ba,Bb.
The orientation of the outer orbit on the sky will also be known from future Gaia data releases because the amplitude of the photocentric orbit is quite large, ∼12 mas. The proper motion anomaly (difference between the short-term proper motion measured by Gaia and the long-term proper motion deduced from the HIPPARCOS and Gaia positions) is quite large, (+19.4, −13.7) mas/yr (Brandt 2018). Moreover, the long-term proper motion deduced from the HIP-PARCOS and Gaia positions, (+5.4, +29.0) mas/yr, is close to (+5.0, +30.2) mas/yr measured by Gaia EDR3 on a 2yr time base that effectively averages the outer orbit, while proposed that the Pout = 180.4 days period outer component of V1200 Cen might also be a binary, forming a more tight half-yearlong period quadruple system. 19 https://www.eso.org/sci/facilities/paranal/ instruments/gravity.html Gaia DR2 proper motion measured on a 1.5 yr baseline is substantially different.
Finally since BG Ind is so bright, and the eclipses of binary A are relatively deep (at ∼15%), we encourage amateurs to continue the eclipse timing. The historical archival data are very helpful in this regard, but not as accurate as targeted observations of this star would be with small to modest sized telescopes. Furthermore, a secure verification of the suspected continuous period change of binary A also needs long-term, continuous follow up timing observations.