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

ABSTRACT

BG Ind is a well-studied, bright, nearby binary consisting of a pair of F stars in a 1.46-d orbit. We have discovered in the TESS light curve for TIC 229804573 (aka BG Ind), a second eclipsing binary in the system with a 0.53-d period. Our subsequent analyses of the recent TESS and archival ground-based photometric and radial velocity (RV) data reveal that the two binaries are gravitationally bound in a 721-d 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 light curve, archival RV data, and the spectral energy distribution, coupled with the use of PARSEC stellar isochrones. We confirm prior studies of BG Ind that found that the brighter binary A consists of slightly evolved F-type stars with refined masses of 1.32 and 1.43 M_⊙, and radii of 1.59 and 2.34R_⊙. The previously unknown binary B has two less massive stars of 0.69 and 0.64 M_⊙ and radii of 0.64 and 0.61R_⊙. Based on a number of different arguments that we discuss, we conclude that the three orbital planes are likely aligned to within 17^○.

1 INTRODUCTION

BG Ind (κ₁ 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, Knigge & Ott (1964), and its EB nature was found by Manfroid & Mathys (1984). At the same time, Andersen, Jensen & Nordstrom (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 light curve analysis was carried out by van Hamme & Manfroid (1988).

In the forthcoming decades, several new photometric and 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 light curves 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 Section 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,¹ thereby indicating clearly that this component has already evolved away from the main sequence and is moving toward the subgiant regime. Secondly, 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), and 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 catalogue of those detached eclipsing binaries for which the constituent masses and radii are known to at least 2 per cent 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 et al. 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 versus π_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 EB in the light curve of BG Ind with a period of 0.53 d.² Our comprehensive investigation of the TESS photometry, archival ground-based photometry and RV curves, as well as the eclipse timing variations (ETV) data demonstrate that the two EBs form a close 2+2 quadruple stellar system with a remarkably short outer period of ∼2 yr.

In Section 2, we describe all the available observational data and their preparation for the complex, joint photodynamical analysis that is discussed in Section 3. Then, the results are discussed and, finally, summarized in Sections 4 and 5.

2 OBSERVATIONAL DATA

2.1 Catalogue data

In Table 1, in addition to other catalogue data, we collected the photometric passband magnitudes of the system from different surveys, e.g. Tycho-2 (Høg et al. 2000), 2MASS (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 Section 3 for details). Together with the passband magnitudes given in Table 1, we list their uncertainties as tabulated in the given catalogues. 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 interpolation 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.

2.2 TESS photometry

The TESS space telescope (Ricker et al. 2015) has observed this target in 2-min cadence mode during Sectors 1, 27, and 28. We downloaded both the simple aperture photometry (SAP) and the pre-search data conditioning SAP (PDCSAP) light curves from the MAST portal.³ We used the SAP light curves 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).

2.3 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) and Collier Cameron et al (2006). From 2012 July, 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, Mazeh & Zucker 2005). Data points more than 5 standard deviations from a phase-binned version of the light curve were rejected and the entire night of data was rejected if more than one-fourth of the observations were identified as outliers based on this criterion.

An 8-d section of the WASP measurements is shown in Fig. 2. Note, we converted the original HJD(UTC) times of the WASP observations into BJD(TDT) for the forthcoming analyses. We also applied the same transformations for all the archival data that we describe below.

2.4 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 light curves 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 Section 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 five more eclipse times (see again, in Section 2.6). The observing procedures and data reduction for these observations are similar to those described in Clausen, Helt & Olsen (2001).

2.5 Disentanglement of the light curves

For the combined analysis of all the observational data, we used the original TESS time-series, i.e. the net light curve of the two binaries together. However, at the start of the analysis, we found it worthwhile to disentangle the light curves 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 light curve with the period of binary A into 1000  equal phase cells. However, while producing the fold for binary A, we excluded those data points that 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 light curve of binary A (see the upper panel of Fig. 3). Then this light curve was removed from the original Sector 1 SAP light-curve 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 light curve. As the result of this removal, we have obtained a new, residual time-series that now mainly contains the light variations of binary B,⁴ without the eclipses and ellipsoidal variations of binary A. Therefore, this light curve can be used for determining the mid-eclipse times of binary B.

In the next step, we folded, binned, and averaged this residual light curve with the period of binary B (see the middle panel of Fig. 3). Finally, we subtracted this folded, disentangled light curve of binary B from the original Sector 1 TESS SAP light curve, 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 light curves, as well (see the 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 time-series, we were able to see it clearly in our disentangled version (see Fig. 4).

Finally, in regard to disentangling the light curves, we have also used a second method that 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 light curves from the two independent methods, and thus we do not show those results here. 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 light curve in Fig. 5 as well for comparison.

2.6 ETV data

2.6.1 TESS ETV results

In order to calculate accurate eclipse times from the TESS light curve, we used the disentangled time-series (see above in Section 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.

2.6.2 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⁵ using a fifth-order polynomial fit. The average of the mid-eclipse times were converted into BJD. Finally, we collected one other eclipse time from 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.

2.6.3 BG Ind ETV results

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 Sections 3 and 4, our detailed analysis robustly confirms this hypothesis.

2.7 RV data

We used three sets of RV data for our analysis. These are as follows: (i) Bakış et al. (2010) have obtained 41 RVs between JDs 2 453 968 and 2 453 996 (i.e. in 2006 August) with the High Efficiency and Resolution Canterbury University Large Echelle Spectrograph (HERCULES) of the Department of Physics and Astronomy, New Zealand; (ii) an additional 23 RVs for both components of binary A between JDs 2 454 363 and 2 454 376 (2007 September/October) were collected by 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⁶ 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 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 centre of mass of the whole quadruple system) together with the best-fitting photodynamical solution (see below, in Section 3) are plotted in Fig. 7.

3 JOINT ANALYSIS OF THE AVAILABLE DATA

We used the software package lightcurvefactory (see Borkovits et al. 2019a, 2020, and further references therein) to carry out a complex spectro-photodynamical modeling of the system based on the data collected in Section 2. lightcurvefactory calculates stellar positions and velocities for each object and emulates the light curves, RV curves, 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⁷ (Bressan et al. 2012) to constrain the stellar parameters theoretically through their evolution tracks, and also to model the combined stellar energy distribution (SED) of the four stars. To solve the inverse problem, the code employs a Markov chain Monte Carlo (MCMC)-based parameter search with an implementation of the generic Metropolis-Hastings algorithm (see e.g. Ford 2005).

Our combined analysis is primarily based on the following observational inputs: (i) the high-quality TESS light curve (see Section 2.2), (ii) the RV data available in the literature (see Section 2.7), (iii) the ETV data calculated from all the available photometric observations (see Section 2.6 and Tables 2–4), and (iv) the observed passband magnitudes of the target taken from standard catalogues (see Section 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 that manifests itself in the form of quadratic deviations of the first few ETV points (see above, in Section 2.6.3). Considering the fact that our analysis depends primarily on the TESS measurements obtained over the last 2.5 yr, and on the RV 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 Section 4.

Regarding the archival photometric observations, we decided not to use the photometric fits to these light curves 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 yr, the relevance of which will be discussed later.

Before the analysis, we also took further preparatory steps on the TESS light curve. 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 light curve, we realized that the residual curve exhibits small amplitude, quasi-periodic variations. We had also found very similar patterns in the residual light curve at the end of the prior light curve disentangling process (Section 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 misadjustments of the light curve 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 light curves from the observed data, the residual light curve 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 light curve is added to the double binary model light curve and the actual χ² value is calculated for this mixed model light curve. Finally, 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 light curve. Therefore, for the main portion of our analysis we used only a 7-d-long section of the Sector 1 TESS light curve.⁸

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 light curve, 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:

• Seven light curve related parameters: the temperature ratios of (T₂/T₁)_{A, B} and T_Ba/T_Aa, 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 (Δt_pri)_{A, B}; the ratios of the radii in both pairs (R₂/R₁)_{A, B}; and the gravity darkening coefficients of the two stars of binary A (β_{Aa, Ab}).

• One parameter for each inner binary orbit, i.e. the observed inclinations i_{A, B} of the orbital planes of binary A and binary B, and five orbital parameters of the outer orbit: period (P_out), time of periastron passage τ_out, eccentricity and argument of periastron (ecos ω)_out and (esin ω)_out, and the inclination i_out.

• Four mass-related parameters: the masses of the two primaries (m_{Aa, Ba}), and the mass ratios of the two binaries (q_{A, B}).

Regarding the other orbital parameters of the inner binaries, the periods (P_{A, B}) of these EBs, as well as their orbital phase (through the time of an arbitrary primary eclipse – |$\mathcal {T}^\mathrm{pri}_\mathrm{A,B}$|⁠) were constrained internally through the ETV data. Furthermore, the eccentricities of both inner orbits were set to zero. Moreover, for the large P_out/P_{A, 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 data set 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 RV of the whole quadruple system (γ) was also constrained internally by minimizing the |$\chi ^2_\mathrm{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 passband-dependent tables downloaded from the phoebe 1.0 Legacy page.⁹ 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 light curve 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.¹⁰ 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 light curve-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.¹¹ 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, determined 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 |$\chi ^2_\mathrm{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 light curves that yield similarly low |$\chi ^2_\mathrm{LC}$| values to the ones obtained in the previous astrophysical model-independent 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 that 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-fitting 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.

4 RESULTS AND DISCUSSION

4.1 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 known stellar masses and other stellar parameters. The outer period of the system is found to be P_out = 721 ± 3 d, which is the shortest amongst doubly eclipsing quadruple systems with an accurately known outer period.¹² Note, however, that despite the relatively short outer period, both outer to inner period ratios are large enough (P_out/P_A ≈ 492, P_out/P_B ≈ 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 eccentric with e_out = 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 Fig. 6. First, due to the unlucky fact that the outer period is nearly exactly equal to 2 yr (P_out ≈ 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. Secondly, 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. |$\dot{P}_{\rm A}$|⁠.

Though the modelling of |$\dot{P}_{\rm 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 pericentre, 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 i_out = 86^○ ± 5^○. On the other hand, the inclinations of the two inner eclipsing binaries are found to be i_A = 73|${_{.}^{\circ}}$|1 ± 0|${_{.}^{\circ}}$|1 and i_B = 84|${_{.}^{\circ}}$|3 ± 0|${_{.}^{\circ}}$|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.)

4.2 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 versus 1.428 ± 0.008 for the more evolved component, and 1.315 ± 0.025 versus 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 |$\lesssim 1\, \sigma$|⁠, 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 that 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 |$R_\mathrm{Aa}=2.34\pm 0.02\, \mathrm{\it R}_\odot$| for the evolved component, while Rozyczka et al. (2011) obtained the somewhat smaller value of |$2.29\pm 0.02\, \mathrm{\it R}_\odot$|⁠.¹³ For the other less evolved star we found |$R_\mathrm{Ab}=1.59\pm 0.04\, \mathrm{\it R}_\odot$|⁠, in contrast to their somewhat larger value of |$1.68\pm 0.04\, \mathrm{\it R}_\odot$|⁠. Keeping in mind, however, that the sum of the (fractional) radii of the two stars (i.e. R_Aa/a_A + R_Ab/a_A) is one of the most robustly determined parameters of an eclipsing binary’s light curve (though it does depend sensitively on the inclination), one can easily check that this sum agrees well in the two solutions.¹⁴

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 |$\approx 2.2{{\ \rm per\ cent}}$|⁠) coming from binary B, and the light-curve distortions caused by its eclipses. Secondly, we used the high-quality TESS photometry whose superiority over the former ground-based measurements is unquestionable. Thirdly, we allowed for the reflection/irradiation effect that 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 four-star SED, resulted in slightly larger temperatures. We found T_Aa = 6442 ± 29 and T_Ab = 6816 ± 26 K in contrast to 6350 ± 260 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) and Stassun & Torres (2016) found a consistent 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 d_DR2 = 51.0 ± 0.5 pc.¹⁵ 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 |$d_\mathrm{HIP}=67.1_{-2.2}^{+2.8}$|⁠, 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 d_{Stassun + 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¹⁶ 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 2 yr (P_out = 1.973 yr) and, furthermore, the semi-major axis of BG Ind A’s ellipse around the centre of mass of the quadruple system is a_{out, A} ≈ 0.82 au. Therefore, the combination of the orbital motion of the photocentre of binary A along the outer orbit and a period near 2 yr 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 q_B = 0.93 ± 0.01. Despite the fact that the flux contribution from this binary is only about 2 per cent in the TESS photometric band, due to the high-quality TESS photometry of this rather bright quadruple, we were able to obtain quite good light curves (see Fig. 3) and a robust dynamical model for binary B. The out-of-eclipse sections of the light curve 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 satisfactory, 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 that induces brightness fluctuations that cannot be well represented by a few smooth harmonics (even after averaging the light curve 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 traveltime 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 main-sequence 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 [Fe/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 per cent 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 possibility, 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.

5 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 P_out = 1.973 yr, we were able to use ground-based archival light curve 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 Catalogue (Tokovinin 2018a), there are only five such systems (including BG Ind) with outer periods shorter than 3 yr. 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 reasons – specifically, 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¹⁷ (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, see Tokovinin 2013 and ζ Cnc C, see Tokovinin 2017).¹⁸

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 per cent, 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 single-dish telescopes a challenging task. A much better prospect is offered by long-baseline interferometers, e.g. the GRAVITY instrument at VLTI¹⁹ (Gravity collaboration 2017). The contrast in the K band is more favourable, compared to the visible (ΔK = 2.98 mag versus Δ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, i.e. ∼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 Hipparcos and Gaia positions, (+5.4, +29.0) mas yr⁻¹, is close to (+5.0, +30.2) mas/yr measured by Gaia EDR3 on a 2-yr time base that effectively averages the outer orbit, while 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 per cent), 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.

ACKNOWLEDGEMENTS

The authors are grateful to Drs Keivan Stassun and Andrei Tokovinin for very valuable discussions.

TB acknowledges the financial support of the Hungarian National Research, Development and Innovation Office – NKFIH Grant KH-130372.

PM is grateful for support from Science and Technology Facilities Council (STFC) research grant number ST/M001040/1.

We thank Allan R. Schmitt for making his light curve examining software ‘lctools’ freely available.

This paper includes data collected by the TESS mission. Funding for the TESS mission is provided by the National Aeronautics and Space Administration (NASA) Science Mission directorate. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts.

This work was partly based on data obtained from the ESO Science Archive Facility under request number 556483.

This work has made use of data from the European Space Agency (ESA) mission Gaia,²⁰ processed by the Gaia Data Processing and Analysis Consortium (DPAC).²¹ Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

We used the Simbad service operated by the Centre des Données Stellaires (Strasbourg, France) and the ESO Science Archive Facility services (data obtained under request number 396301).

DATA AVAILABILITY

The TESS data underlying this paper were accessed from MAST (Barbara A. Mikulski Archive for Space Telescopes) Portal (https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html). A part of the data were derived from sources in public domain as given in the respective footnotes. The derived data generated in this research and the code used for the photodynamical analysis will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES