BU Canis Minoris – the most compact known flat doubly eclipsing quadruple system

ABSTRACT

We have found that the 2+2 quadruple star system BU CMi is currently the most compact quadruple system known, with an extremely short outer period of only 121 d. The previous record holder was TIC 219006972 (Kostov et al.), with a period of 168 d. The quadruple nature of BU CMi was established by Volkov, Kravtsova & Chochol, but they misidentified the outer period as 6.6 yr. BU CMi contains two eclipsing binaries (EBs), each with a period near 3 d, and a substantial eccentricity of ≃0.22. All four stars are within ∼0.1 M_⊙ of 2.4 M_⊙. Both binaries exhibit dynamically driven apsidal motion with fairly short apsidal periods of ≃30 yr, due to the short outer orbital period. The outer period of 121 d is found both from the dynamical perturbations, with this period imprinted on the eclipse timing variations curve of each EB by the other binary, and by modelling the complex line profiles in a collection of spectra. We find that the three orbital planes are all mutually aligned to within 1 deg, but the overall system has an inclination angle near 83.5°. We utilize a complex spectro-photodynamical analysis to compute and tabulate all the interesting stellar and orbital parameters of the system. Finally, we also find an unexpected dynamical perturbation on a time-scale of several years whose origin we explore. This latter effect was misinterpreted by Volkov et al. and led them to conclude that the outer period was 6.6 yr rather than the 121 d that we establish here.

1 INTRODUCTION

There are currently more than 300 known 2+2 quadruples consisting of an orbiting pair of eclipsing binaries (EBs), most of which have been found with Transiting Exoplanet Survey Satellite (TESS; Kostov et al. 2022; Zasche, Henzl & Mašek 2022; Kostov et al. 2023). The criteria for accepting these as quadruples are: (i) there are two EBs that are (ii) unresolved at the pixel level with TESS, and (iii) which show only one dominant star in Gaia within the TESS pixel. However, given that Gaia does not often distinguish stars that are ≲1/2 arcsec apart, and that these objects are typically a kpc away, this implies only that the physical separation of the EBs is ≲500 au or so. The corresponding outer orbital periods are only constrained to an order of ≲5000 yr.

At the largest of these orbital separations, the quadruples would be too wide for easy-to-measure dynamical interactions that could lead to a determination of the outer orbital period. However, a small percentage of these quadruples have much closer separations of less than a few au, and these have led to measurable outer orbits as well as other interesting dynamical interactions, such as TIC 454 140 642 (432 d; Kostov et al. 2021), TIC 219 006 972 (168 d; Kostov et al. 2023), and VW LMi (355 d; Pribulla et al. 2008).

With quadruples having outer orbital periods as short as approximately a couple of years, interesting and informative dynamical interactions to look for include: (i) dynamical delays resulting from a changing period of either EB due to the varying distance to the other EB; (ii) dynamically forced orbital precession in the EBs; and (iii) forced precession of the orbital planes leading to eclipse depth variations.

Every once in a while, one of these quadruples turns out to have a dramatically short outer period (e.g. TIC 219006972; Kostov et al. 2023) and it becomes quite feasible to use eclipse timing variations (ETVs) data as well as radial velocity (RV) data to completely diagnose most of the important stellar and orbital parameters. In this work, we report on BU CMi, a 2+2 quadruple which we find to have the shortest known outer orbital period of 121 d, and a very interesting array of dynamical effects.

This work is organized as follows. In the remainder of Section 1, we discuss how the quadruple nature of BU CMi came to be known. Section 2 briefly reviews the observations that we bring to bear on the analysis of BU CMi, including with TESS and follow-up ground-based photometry and spectroscopy. In Section 3, we briefly describe our use of broadening functions to extract RVs from the spectra. Section 4 details our production of a long-term ETV curve from archival as well as TESS data, and what we can learn from a visual inspection of it. In Section 5, we explain in some detail how we fit a model of stellar and orbital parameters directly to the complex and overlapping line profiles in each spectrum. Section 6 is devoted to a discussion of disentangling the four spectra, after we determined the stellar and orbital parameters. We review our spectro-photodynamical model for ascertaining the system parameters in Section 7. A discussion of our work and its implications follows in Section 8. We give some concluding remarks in Section 9.

1.1 Prior work on multiple stars

The TESS (Ricker et al. 2015) has been instrumental in the discovery of multiple star systems. There has been significant progress on the problem of identifying triple, quadruple, quintuple, and even sextuple star systems from TESS data; for instance, the Planet Hunters citizen science project uncovered tens of multiple star system candidates (Eisner et al. 2021). In parallel, a combination of machine learning techniques and human vetting led to the creation of an extensive catalogue of quadruple stellar systems (Kostov et al. 2022), in addition to the identification of the first-ever sextuply eclipsing six star system (Powell et al. 2021).

1.2 BU Canis Minoris

BU CMi (HD 65241, HIP 38945, TIC 271204362) is listed as a suspected Algol-type eclipsing binary in the 74th special namelist of new variable stars (Kazarovets et al. 1999). In spite of its brightness (V = 6.42), the object was not subject to any detailed study until Volkov et al. (2021) found it to be a quadruple system composed of two EBs with P_A = 2.94 d and P_B = 3.26 d. Although the authors had observed BU CMi with a photoelectric photometer in 2012, they first noticed two systems of eclipses in MASCARA photometric data in 2020. According to their analysis, all four stars in the quadruple are similar (A0 spectral type), having masses in the range 3.1–3.4 M_⊙. The orbital eccentricities of the inner binaries were found to be relatively high for close binaries, with e_A = 0.20 and e_B = 0.22. The authors interpreted the observed ETVs as being due to (i) relatively rapid apsidal motion in the inner (i.e. binary) orbits with the periods U_A = 25.4 yr and U_B = 26.3 yr, and (ii) light travel time effects (LTTE) from the outer orbit between the two binaries having a period of 6.62-yr and a high-eccentricity (e_AB = 0.7).

Gaia DR3 (Gaia Collaboration 2023) lists BU CMi as a single object astrometrically and spectrosopically with a parallax of π = 4.0143 ± 0.0335 mas and the following atmospheric parameters |$T_{\rm eff} = 10173^{+43}_{-39}$| K, |$\log g = 3.727^{+7}_{-6}$|⁠, and [M/H] = 0.778|$^{+17}_{-40}$|⁠.

The main catalogued photometric and kinematic data for BU CMi are given in Table 1.

1.3 Our two independent discoveries of the quadruple nature of BU CMi

The Visual Survey Group (VSG), in its search for compact multistellar systems, continues to visually survey large numbers of light curves from the TESS mission (Kristiansen et al. 2022). Some of its findings, including in the area of multistellar systems, are given in Kristiansen et al. (2022) and Rappaport et al. (2019, 2022). In 2021 March, the group spotted BU CMi in the Sector 34 light curves and immediately identified it as a potential quadruple star system.

In parallel to the VSG’s survey of the light curves from the TESS full-frame images (FFIs), RJ and SR have also been searching the 2-min and 20-s cadence light curves for strongly periodic or time-varying stellar phenomena. Every 2-min-cadence light curve is passed through an algorithm that flags it for further review if there is a detection of at least a 12σ peak in its periodogram. When searching the Sector 34 short-cadence light curves through this algorithmic process, BU CMi triggered on this algorithm; upon human review, this light curve was found to be a bona fide quadruple system and flagged for further follow-up.

Since that time, we have continued to collect information on, and model, the BU CMi system. After 2 yr of study, we now report here on the discovery of a 121 d outer period for the system.

2 OBSERVATIONS

2.1 TESS observations

BU CMi (TIC 271204362) was observed during TESS Sectors 7, 34, and 61 (from 2019 January 7 to 2019 February 2, 2021 January 13 to 2021 February 9, and 2023 January 18 to 2023 February 12, respectively). In Sector 7, it was observed at 30-min cadence in the FFIs, while in Sectors 34 and 61, besides the 600 and 200-s cadence FFIs, it was also observed at 2-min cadence. For most of our study we used the FFI observations, from which the light curves were processed using the fitsh pipeline (Pál 2012). We used the 2-min cadence SAP-FLUX light curves, which were downloaded directly from the Barbara A. Mikulski Archive for Space Telescopes (MAST) website, only for the determination of mid-eclipse times over Sectors 34 and 61 data. Naturally, in the case of the Sector 7 observations, we could determine the mid-eclipse times only from the sparsely cadenced FFI data at 30 min.

We show some illustrative segments of the three TESS sectors in Fig. 1. The superposed model curves will be discussed in Section 7.

2.2 Ground-based photometry

BU CMi is bright enough (M_V = 6.42) that it can be reliably observed using small-aperture telescopes from the ground. We used archival data obtained from the Kilodegree Extremely Little Telescope (KELT–2600 points; Pepper et al. 2007, 2012), the Multi-site All-Sky CAmeRA (MASCARA–10800 points; see Talens et al. 2017 and Burggraaff et al. 2018), and the Hungarian Automated Telescope Network (HATnet–3200 points; Bakos et al. 2004). BU CMi was heavily saturated in the HATNet observations, so we performed a custom analysis of these data that differs from the standard methods used by the survey. We extracted aperture photometry for BU CMi and 200 other comparably bright stars through an annulus excluding the saturated core of the target. We then performed an ensemble correction using a linear combination of the 200 bright neighbours to correct for instrumental and atmospheric variations in the resulting light curve.

Additionally, we carried out further photometric follow-up observations between 2021 March and 2023 April. These were mainly obtained with a small 34-mm telescope in a private observatory in Jílové u Prahy in the Czech Republic, as well as remotely in Northern Italy (both by one of the co-authors RU). Another observing site was in Argentina as part of the Pierre Auger Observatory (Aab et al. 2021), and obtained by MM. On three additional nights, the target was also observed by ZH from two observing sites in the Czech Republic. Finally, two further eclipses of BU CMi were observed with the RC80 telescope of Baja Astronomical Observatory in 2022 and 2023.

2.3 Spectroscopy from Skalnaté Pleso observatory

High-dispersion spectroscopy was obtained with a 1.3 m, f/8.36 Nasmyth–Cassegrain telescope equipped with a fibre-fed échelle spectrograph at the Skalnaté Pleso (SP) Observatory, Slovakia. Its layout follows the MUSICOS design (Baudrand & Bohm 1992). The spectra were recorded by an Andor iKon-L DZ936N-BV CCD camera with a 2048 × 2048 array, 13.5 μm square pixels, 2.9 e⁻ readout noise, and a gain close to unity. The spectral range of the instrument is 4250–7375 Å (56 échelle orders), with a maximum resolution of R = 38 000. Because of the relatively long orbital period of both inner binaries, three 600-s exposures were combined to increase the SNR and to clean cosmic ray hits. The raw spectroscopic data were reduced as in Pribulla et al. (2015) using iraf package tasks, linux shell scripts, and fortran programs. In the first step, master dark and flat-field frames were produced, based on the spectra from the tungsten lamp and blue LED. In the second step, the photometric calibration of the frames was performed using dark and flat-field frames. Bad pixels were cleaned using a bad-pixel mask, and cosmic ray hits were removed using the program of Pych (2004). Order positions were defined by fitting sixth-order Chebyshev polynomials to tungsten-lamp and blue LED spectra. Subsequently, scattered light was modelled and subtracted, and then aperture spectra were extracted for both BU CMi and the ThAr lamp. The resulting two-dimensional spectra were then dispersion solved and combined to one-dimensional spectra. Finally, all spectra were continuum normalized. The typical RV stability of the spectrograph is 200 m s⁻¹.

In total, 60 spectra with per-pixel SNR ranging from 45 to 175 (at 5500 Å) were obtained from 2020 March 17 through 2023 March 3.

2.4 CTIO spectroscopy

Additional spectroscopy was obtained with the CHIRON fibre-fed échelle spectrograph at the 1.5m telescope of the Cerro-Tololo Interamerican Observatory, Chile. A detailed description of the spectrograph can be found in Tokovinin et al. (2013). All spectra were taken in the slicer mode providing R = 80 000. The spectrum was extracted from 70 échelle orders covering 4080–7000 Å. Similar to the case of the SP spectroscopy, three consecutive 600-s exposures were combined. Even without an iodine cell, the spectrograph stability is a few m s⁻¹. Because of a reflection causing a blemish in the blue part of the spectrum, an independent pipeline based on the iraf scripts has been developed to use the full spectral range of the instrument. In addition to the reduction steps taken for the SP spectra, the raw frames were trimmed and corrected for the overscan. The relative responses of four quadrants as read out by different amplifiers were taken into account. In total, BU CMi was observed on 16 nights from 2022 November 11, until 2023 March 29. The per-pixel SNR ranged from 70 to 290 (at 5500 Å).

2.5 Konkoly and Rozhen spectroscopic observations

Additional spectroscopic measurements were obtained with the R ∼ 20 000 échelle spectrograph mounted on the 1m RCC telescope at Konkoly Observatory, Hungary between 2021 December and 2022 March. The spectrograph is capable of covering the 3890–8670 Å wavelength range in a set of 33 échelle orders. The images were taken with a back-illuminated FLI ML1109 CCD camera having an array of 2048 × 506 pixels of size 12 μm, 10 e⁻ readout noise, and a gain close to unity. A total of 11 spectra were observed, with exposure times varying between 900 and 3600 s. Two additional spectra were taken with the 2 m RCC telescope at NAO Rozhen, Bulgaria equipped with the R ∼ 30 000 ESpeRo spectrograph in 2022 January. The details of the instruments can be found in Bonev et al. (2017). The spectra were reduced completely in the same manner as in Mitnyan et al. (2020) including the steps of bias, dark and flat corrections, wavelength calibration, continuum normalization, telluric line removal, and barycentric correction. The per-pixel SNR of these additional observed spectra is between 70 and 140.

3 RADIAL VELOCITIES

We extracted the broadening functions from all spectra of BU CMi using the publicly available software bf-rvplotter¹ in a ±300 km s⁻¹ velocity range using 3 km s⁻¹ RV bins, and a Gaussian smoothing using a rolling window of five datapoints. We tried templates of different (G-, F-, A- spectral type) stars, and the spectrum of Vega with A0V spectral type produced the best quality BFs over the 4400–4600 Å  and 4900–5300 Å wavelength range that we used for the calculations. The sum of four Gaussian functions was then fitted to the resulting broadening functions to find the peak positions of the profiles of the component stars in order to determine their RVs. Although the use of Gaussian functions is not the most appropriate choice for fitting the BFs of rapidly rotating stars, we found it to work well in the case of the BFs of BU CMi, as we were able to constrain satisfactory model fits to them using the sum of four Gaussian functions. An example of a typical BF along with the corresponding model fit is displayed in Fig. 2. The resulting RVs and their 1σ uncertainties are listed in Table 2.

4 PERIOD STUDY

4.1 Determination of times of minima

In order to carry out a preliminary period study with the usual method of analysing the ETVs, we determined accurate mid-eclipse times from the various photometric light curves. Due to the presence of frequently and strongly overlapping eclipses in the two EBs, this process required some extra care.

In the case of the quasi-continuous, 25–26-d-long TESS observations, we disentangled the light curves of the two EBs, and mutually removed the signals of the other EB from the light curve in the same manner as was also done for our analyses of V994 Her (Zasche et al. 2023). This needs to be done in order to avoid losing the overlapping eclipse events and, moreover, to enhance the accuracies of those times where the eclipses of the two EBs do not overlap each other, but other light-curve features (e.g. the periastron bumps) of one of the binaries distorts the eclipsing signal of the other EB. Due to the rapid and large amplitude ETVs, this method was done separately for the three sectors and, even then, in the case of a few fully overlapping eclipsing events, the disentanglement was not perfect. As a result, we had to drop out two to three outlier eclipse times from our analysis.

In the absence of well-covered out-of-eclipse light curves for our 2021–2023 targeted ground-based follow-up eclipse observations, we were unable to carry out the same disentanglement process, and had no choice but to simply exclude a few blended eclipses from our ETV analysis.

In the case of the MASCARA, KELT, and HAT data, after their conversion from the original Heliocentric Julian Date (HJD) times to Barycentric Julian Date (BJD), we again disentangled the two binaries’ signals from the light curves. We then took those eclipse observations where both the ingress and egress portions of the same eclipses were measured, and we determined mid-eclipse times for these individual eclipses. The resulting mid-eclipse times as determined from the TESS, KELT, and MASCARA archival data, as well as our new ground-based observations, are listed in Tables 3 and 4.

Finally, we utilized some additional published eclipse times from the papers of Ogloza et al. (2008) and Volkov et al. (2021). Regarding this latter work, Volkov et al. (2021) tabulate their own eclipse times determined from the MASCARA data and the TESS Sectors 7 and 34 observations. We do not use those particular times, as we have determined our own eclipse times from the same data sets using our own method. We do use, however, most of the eclipse times that Volkov et al. (2021) determined from their own observations, and to which we do not have access. Unfortunately, they do not provide uncertainties for these times, so we arbitrarily take an estimated error of 0.0001 d for each of these data points.

The combination of all these data has provided 102 and 81 mid-eclipse times for binaries A and B, respectively.

4.2 Preliminary ETV study

The ETV curves formed from all the data listed in Tables 3 and 4 are shown in Fig. 3. The top two panels show the overall ETV curves spanning all the data. The left and right panels are for binary A and binary B, respectively. The bottom panels are zoom-ins around the three TESS sectors and one ground-based segment.

A first visual inspection of the ETV curves of both binaries shows two very remarkable features. First, each EB exhibits a large sinusoid with an amplitude of ∼±0.2 d, and periods of 25–30 yr. Note that the corresponding curves for the primary and secondary eclipses are anticorrelated. This clearly indicates apsidal motion of the eccentric EBs. We will show in Section 7 that this apsidal motion is driven by the presence of the ‘other’ binary.

Secondly, the high-quality modern ETV points determined from the TESS and ground-based follow-up observations exhibit clear ∼120-d periodic variations, with amplitudes of 30–40 min. The shapes of these short-period (∼120-d) ETVs clearly resemble the dynamically driven ETVs of several other tight and compact triple and quadruple star systems (see e.g. Borkovits et al. 2015, 2016, 2020; Kostov et al. 2021, 2023). We thereby associate the ∼120-d ETV features with the outer orbital period of the quadruple in which binary A orbits binary B. Moreover, the apsidal motion periods are in perfect agreement with periods which are theoretically expected for tight triple systems with the observed inner and outer periods. It was therefore clear to us that BU CMi is an exceptionally compact doubly eclipsing 2+2 hierarchical quadruple system.

As a first attempt to understand these ETVs quantitatively, we made analytic ETV fits using the code and method described in Borkovits et al. (2015). As this code was developed for tight triple systems, unfortunately, it is unable to simultaneously fit the ETVs of both binaries. On the other hand, however, the analytical descriptions of the P_out-time-scale, as well as the apse-node perturbation terms, are included in the software package and, hence, it is very useful for getting reliable third-body solutions for dynamically dominated ETV curves. When used to analyse the BU CMi ETV curves, the solutions confirmed the P_out ≈ 121.5 d outer orbital period and, moreover, the fact that the observed rapid apsidal motions are driven by the third-body perturbations of the other EB.² While our analytic ETV fits confirmed our initial hypothesis about the system configuration, the obtained fit was surprisingly poor for both binaries. This inadequate fit is due to an ‘extra feature’ in the ETV curve with a period of ∼1000–1200 d that can be seen by casual inspection of the curves.

We were able to improve our analytical ETV fit substantially, however, when we ‘added’ an additional, more distant fourth (i.e. in the present situation, fifth) stellar component to the system. A recent improvement of our analytic ETV software package makes it possible to fit simultaneously a second LTTE. We found that adding such a second LTTE component with a period of P_outermost ≈ 1100 d results in statistically excellent fits for the ETV curves of both binaries. We emphasize, however, that this was a solution that was useful in the context of the fitting code, rather than having any physical significance. Such an extremely tight 2+2+1 quintuple system would be at best marginally stable; furthermore, the strong gravitational perturbations of the outermost component should also have to be taken into account. Later, we demonstrate, when discussing the fully consistent, photodynamical treatment (see Section 7), that the explanation for the origin of the extra cyclic terms in the ETVs with a ≈1000–1200-d period does not require the presence of an additional, more distant stellar component in the system. Rather, we show that the extra cyclic feature in the ETV curve arises naturally from the mutual gravitational perturbations between the two binaries.

Finally, we note that the results of this preliminary, analytic ETV analysis were used only for finding reliable input parameters for the complex, photodynamical fitting procedure (Section 7). Thus, a more detailed, quantitative discussion of the ETV results will be given later, in Section 8, in the context of the photodynamical results.

5 DIRECT FITTING TO THE SPECTRA AND TESS LIGHT CURVES

The ETV analysis presented in Section 4.2 strongly indicates that the outer orbit of the system is very tight, with P_out ≈ 121.5 d. This means that the mutual orbit of the EBs is relatively fast, resulting in EB centre of mass RV semi-amplitudes of the order of tens of km s⁻¹. Such a large RV variability superposed on EB orbits of a few days should be easy to detect.

The analysis of the spectra of BU CMi is, however, significantly complicated by the fact that all its stellar components possess an early spectral type. The spectra are therefore dominated by strong, wide hydrogen Balmer lines; these lines are very difficult to work with, as we cannot reliably detect their splitting and robustly calculate RVs of the individual stars in the system. On the other hand, the metal lines of the system are rather shallow; as a result, a relatively high signal-to-noise ratio is necessary to use these in our analysis. The line splitting due to the orbital motion in both pairs is clearly visible in the strongest metal line (Mg ii 4481 Å). Because the rotational velocities of the components (⁠|$v\, \sin \, i\sim$| 50–80 km s⁻¹) are approximately half the RV semi-amplitudes (K ∼ 120–130 km s⁻¹), the components’ profiles are almost always somewhat blended. Thus, a direct measurement of the RV is rather difficult and may yield inaccurate results (see Section 3).

Consequently, we attempted to directly model the observed spectra, under the assumption of two binary stars in a Keplerian orbit about each other. Two wavelength ranges were found to be appropriate: a blue region around the Mg ii 4481 Å line, covering 4380–4605.4 Å, and a green region around the Mg i triplet, covering 4900–5379.6 Å. The spectra were weighted using the signal-to-noise ratios provided by the reduction pipeline at both 4400 and 5500 Å. The standard deviation at either of the spectral ranges was assumed to be 1/SNR.³ Before being used for modelling, the original spectra were rebinned to a logarithmic wavelength grid with an RV step of 3.5 km s⁻¹ in both spectral ranges. Due to the small number of observations during epochs in which either of the systems was in an eclipse, no spin-axis orbital-plane misalignment for either of the components was assumed; thus, we set the parameter λ = 0 deg.

The sum of synthetic spectra for all stellar components was fitted to the observed spectrum at each epoch. Prior to summing, the synthetic spectra of individual components were convolved with a theoretical limb-darkened rotational profile. Mutual eclipses of the components were taken into account. Keplerian motion in both the inner and outer orbits was assumed. Synthetic spectra were computed using ispec (Blanco-Cuaresma et al. 2014; Blanco-Cuaresma 2019), a program based on the spectrum code (Gray & Corbally 1994). For the spectra, we assumed a solar metallicity (log [m/X] = 0.0) and that 9000 K ≤ T_eff ≤ 12 000 K. We also allowed for apsidal motion for the inner orbits. No other gravitational perturbations were assumed in this relatively simple model. We found that the principal atmospheric parameters corresponding to the best model for either of the spectral ranges are log g = 4.5 for all four stars; T_eff = 11 000 K for components 2, 3, and 4; and T_eff = 11 500 K for component 1.⁴

We modelled the blue and red spectral ranges separately. In order to constrain the fit parameters – specifically the inclination angle i, orbital periods of inner binaries, P_A, P_B, and the apsidal motion rates of the binaries P_aps – we included the TESS photometry from Sectors 7 and 34 in the modelling. The ellipsoidal light variations were modelled using the analytical approximations of Esteves, De Mooij & Jayawardhana (2013); these were found to be sufficiently accurate for both inner binaries. While the semi-amplitudes of the RV changes K_i are free parameters, the semi-amplitude of the radial-velocity variations of the centre of masses of the binaries depend on the masses of the individual components (M_i, i = 1, 2, 3, 4), the inclination angle i_out, eccentricity e_out, and the orbital period P_out of the outer orbit. The semi-amplitudes of these RV changes, K₁₂ and K₃₄, were hence computed using these parameters.

The simultaneous modelling of the TESS light curve and the spectra was performed via gradient-based optimization starting from multiple sets of trial parameters in order to arrive at the global minimum of χ². Initial optimization runs showed that the inclination of the outer orbit is close to 90°, but attempts to adjust it led to divergent results. Hence, the outer orbit was fixed at edge-on and the inclination was not adjusted. We were able to constrain the outer orbit and demonstrate that it is extremely tight, with P_out ∼ 121 d. This differs significantly from the estimate provided by Volkov et al. (2021), who cite a value over 20 times ours (P_out = 2420 ± 40 d or 6.62 yr). In Section 8, we discuss the differences in analyses, and lay out a case for why 121 d is definitely the period of the outer orbit.

The parameters corresponding to the best models for the blue and red regions are listed in Table 5, and the best-fitting models to the first 20 observed spectra of BU CMi in the blue wavelength range, close to the Mg ii 4481 Å line, are plotted in Fig. 4. The best-fitting stellar parameters for the component stars are listed in Table 6. For further analysis, we will equate the blue and yellow spectral ranges with the Johnson B and V passbands, and use B and V magnitudes where relevant.

6 SPECTRA DISENTANGLING AND ATMOSPHERIC PARAMETERS

After obtaining a set of model RVs for both stellar components at all epochs of the spectroscopic observations, we then attempted to disentangle the composite spectrum in order to obtain individual spectra for each star (1, 2, 3, and 4). The iterative disentangling approach of Bagnuolo & Gies (1991) was used, and the spectra were disentangled in both the blue and red ranges. As previously, a radial-velocity step size of 3.5 km s⁻¹ was used for both spectral ranges. The process converged to a solution for both spectral ranges, and the resulting composite spectra are plotted in Fig. 5. These spectra show small imperfections, perhaps from issues with the continuum rectification and artefacts from the disentangling.

However, our technique does not yield the flux ratio of the components (see e.g. Pavlovski & Hensberge 2010; Lehmann & Tkachenko 2012). As a result, we corrected the component spectra for the contribution of the other components before proceeding. Properly normalizing the spectra will have an impact on the line depths and is possible using the component passband-specific luminosities derived from simultaneously modelling the spectra and the light curve. Computed brightnesses for the B and V Johnson bands were assumed for the blue and red spectral ranges, respectively. As part of our modelling, we used the radii derived in Table 5 for the correction step and to determine the atmospheric parameters T_eff, log g, metallicity log [m/X]. The projected rotational velocities were taken from the results of the spectral fitting. As before, we used the ispec code.

The temperatures and surface gravities for binary A are T_pri = 10 760 ± 360 K, log g_pri = 4.23 ± 0.29, T_sec = 10 180 ± 430 K, log g_pri = 4.11 ± 0.30. The corresponding parameters for binary B are T_pri = 10 820 ± 380 K, log g_pri = 4.63 ± 0.26, T_sec = 10 040 ± 350 K, log g_pri = 4.18 ± 0.23. The metallicity was fixed at the solar value. The component temperatures are slightly lower than those of the best-fitting templates. We also highlight our finding that component 4, in binary B (P = 3.26 d), is cooler than the remaining components. As part of our presentation of these results, we add the caveat that the atmospheric parameters strongly depend on the technique used to normalize the individual disentangled spectra. Using different flux ratios would change the line depth and result in different temperatures and/or metallicities for each of the components in the system.

7 SPECTRO-PHOTODYNAMICAL ANALYSIS

Independent of the spectroscopic analyses described in Section 5, we carried out a simultaneous, joint analysis of the TESS light curve, the ETV curves calculated from the mid-times of all the observed eclipses (i.e. from TESS and ground-based measurements; see Section 4.1), RVs (determined by Gaussian-profile fitting to the BFs; see Section 3), as well as all the available multiwavelength SED data for BU CMi using the lightcurvefactory software package (Borkovits et al. 2019, 2020). During the analysis we followed the same steps as for TIC 454 140 642 and TIC 219 006 972 (Kostov et al. 2021, 2023), two other compact quadruple systems discovered with TESS. For further details of the photodynamical analysis, we refer the reader to section 5.1 of Kostov et al. (2021). Here, we describe only the steps in the data preparation that are specific to this system.

In order to retain similar sampling and to reduce computational costs, we binned the 2-min cadence data sets from Sectors 34 and 61 to 1800 s, similar to the Sector 7 FFI-cadence light curve. When comparing some preliminary fits to the short- and long-cadence light curves, we found no statistically significant discrepancies in the resultant parameters and, hence, we decided to use exclusively the 1800-s cadence light curves for all three available TESS sectors.

During our earlier initial analysis runs (made in 2021 and 2022), we also used a second set of light curves compiled from the ground-based photometric follow-up observations. However, when the Sector 61 observations of the TESS spacecraft became available and, hence, the high-quality and homogeneous satellite photometry then covered more than a 4-yr-long interval, we decided to no longer use the more inhomogeneous and less accurate ground-based photometric data for direct light curve fitting. Naturally, the eclipse times derived from these ground-based observations did continue to be used for the ETV analysis part of the simultaneous, joint spectro-photodynamical fitting process.

Because our software package lightcurvefactory is unable to handle direct fitting of the spectral lines (see Section 5) and fits only direct RV data, we used a second, different approach. We mined out model-independent RV data from all available spectra in the manner described briefly in Section 2.5. Due to the highly blended spectral lines, the correct identifications of the individual stellar components were quite problematic in the case of several spectra. Hence, after obtaining the first model fits, we had to change the labelling of the stars in the RV data in several cases. The RV data used for the spectro-photodynamical analysis are listed in Table 2.

In addition to the TESS light curves and the RV data, as was mentioned before, we used the four ETV curves (primary and secondary ETV data for both binaries; see Tables 3 and 4).

Finally, we also used the catalogue passband magnitudes listed in Table 1 for the SED analysis. For this analysis, similar to our previous works, we used a minimum uncertainty of 0.03 mag for most of the observed passband magnitudes. This was done in order to avoid an outsized contribution from the extremely precise Gaia magnitudes, as well as to counterbalance the uncertainties inherent in our interpolation method during the calculations of theoretical passband magnitudes that are part of the fitting process. The only exception is the WISE W4 magnitude, for which the uncertainty was set to 0.3 mag.

Table 7 lists the median values of the stellar and orbital parameters of the BU CMi quadruple system that have been either adjusted, internally constrained, or derived from the MCMC posteriors, together with the corresponding 1σ statistical uncertainties. Sections of the light curves, the ETV and RV curves of the lowest |$\chi ^2_\mathrm{global}$| solution are plotted in Figs 1, 3, and 6, respectively.

One caveat should be noted regarding the proper interpretation of the orbital parameters listed in Table 7. The tabulated orbital elements, with a few exceptions, are so-called ‘instantaneous osculating elements’, which are valid at the moment of the cited epoch t₀. Thus, they cannot be simply compared with those orbital elements that are deduced directly from photometric and/or spectroscopic observations. These latter orbital elements can be considered to be some kind of long-term averaged orbital elements and connected to such Keplerian orbits that represent the approximations of the time-averaged envelopes of the true (non-Keplerian) motions. This question was discussed in more detail in section 5.1 of Kostov et al. (2021).

Regarding the exceptions, in the (second) row |$\overline{P}_\mathrm{obs}$| we give the average or, ‘observable’, periods for the three orbits (A, B, AB) which were obtained from a longer-term numerical integration initiated with the parameters of the best-fitting spectro-photodynamical model. The values given for the inner binaries (A, B) stand for the long-term average of their eclipsing periods (which technically means that calculating the ETV curves with these periods, the averages of the ETVs over full apsidal motion cycles remain constant). In the case of the outer orbit, however, |$\overline{P}_\mathrm{obs}$| was determined as the time average of consecutive periastron passages, i.e. this is the average anomalistic period of the outer orbit.

The other exception is the first row of the apsidal motion related parameters, |$P_\mathrm{apse}^\mathrm{obs}$|⁠. Here, the duration of an apsidal motion cycle, i.e. the time needed for the complete, 360° variation of the observable arguments of periastron of the three orbits, ω_{A, B, AB}, are given. The other apsidal motion parameters are calculated internally by the software package lightcurvefactory with the use of the usual analytic formulae (see Kostov et al. 2021, section 6.2, for details).

8 DISCUSSION AND IMPLICATIONS

8.1 System parameters

Our two independent analyses (Sections 5–7) have resulted in very similar results, at least qualitatively; however, quantitatively, a number of the discrepancies exceed the estimated uncertainties. For example, both models agree that BU CMi consists of four very similar hot stars. Moreover, the mass ratio of the two primaries (i.e. m_Ba/m_Aa) are ∼0.993 and ∼0.992, according to the direct spectral fitting and the spectro-photodynamical solutions, respectively (i.e. they agree well within 0.1 per cent). The direct spectral fitting approach, however, resulted in systematically more massive (by ∼4–5 per cent) and hotter (by ∼3–8 per cent) components than the spectro-photodynamical analysis. At the moment we are unable to resolve completely these discrepancies. One should keep in mind, however, that the individual masses are primarily controlled by the RV data which, in the current situation, are quite uncertain for the reasons discussed in Sections 3 and 5. Regarding the temperatures, our attempts to obtain reliable spectroscopic solutions with the use of slightly lower temperature templates have failed. On the other hand, if one compares the photometric distance calculated from the spectro-photodynamical solution d_phot = 245 ± 8 pc with the Gaia DR3-derived one d_GaiaDR3 = 247 ± 2 pc (Bailer-Jones et al. 2021) the agreement looks perfect. In contrast to this, in the case of the direct spectral-fitting solution which gives hotter and more massive (and, hence, larger) stellar components, the total luminosity of the quadruple was found to be larger by ∼66 per cent (compare the derived individual luminosities in Tables 6 and 7) resulting in a |${\sim}29~{{\ \rm per\ cent}}$| larger photometric distance. However, one should keep in mind that in the Gaia DR3 astrometric solutions, stellar multiplicity was not taken into account and hence the Gaia-derived distance might be subject to some systematic errors. Thus, in conclusion, we believe that we are not currently in a position to prefer one or the other solution; hence, we conclude that the true uncertainties of our solutions, at least in the masses and temperature, must be around 5–7 per cent. The quantitatively similar results of the two independent approaches clearly show that the parameters are robust and BU CMi is indeed a very tight quadruple with ∼120 d outer orbital period. The differences in the orbital and component parameters between the two analyses give us, however, an independent check of the parameter uncertainties. The similar brightness of the inner binaries makes the astrometric signal of the system’s photocentre tiny; thus we cannot expect that Gaia would provide independent orbital parameters. Because the photodynamical analysis is more complex and in addition to RVs it takes into account all available photometric data (including observed colours and brightness), we prefer the corresponding parameters.

Turning to the orbital geometry of the system, we find that BU CMi is not only the most compact but also an extremely flat, 2+2 type quadruple system. All three orbital planes (that of the two inner EBs and also the outer orbital plane) are well aligned within 1°. Note, this is a common feature of the previous two TESS-discovered compact doubly eclipsing 2+2 type quadruple systems: TIC 454 140 642 (Kostov et al. 2021) and TIC 219 006 972 (Kostov et al. 2023). In drawing any general conclusions from the common flatness of these new compact 2+2 quadruples, however, one should keep in mind that these findings might be biased by observational selection effects, at least, in part. This is so because, in the case of substantially inclined orbital planes, the three planes would precess with different amplitudes and periods on time-scales of decades or centuries (depending on the mutual inclination angles), and this would make it very unlikely that both inner binaries would show eclipses at the same time. And, in the absence of eclipses in either or both of the binaries, there would be no (or only a very minor) chance of detecting the given target as a compact 2+2 system. On the other hand, however, such a selection effect does not explain the extreme flatness (i.e. orbital alignments within 1–2°). For example, the current binaries A and B would already produce eclipses for inclinations i_A ≳ 72.9° and i_B ≳ 74.3°, from which it follows that, keeping the inclination of the system’s invariable plane at the currently derived value of i₀ = 83.8° ± 0.2°, both inner binaries would continuously produce well-detectable eclipses in the case of mutual inclinations of, let us say, i_{A − AB; B − AB} = 5°. So, in conclusion, in our view the very strong flatness is likely a consequence of the formation processes of the most compact 2+2 systems, but, due to the strong observational bias described above, this conclusion is not highly robust.

Besides the similarities of BU CMi to the previously mentioned two quadruples, there are strong differences as well. For example, in the case of the other two quadruples, the inner binaries are in nearly circular orbits, while the inner binaries of BU CMi have remarkably substantial eccentricities (e_A = 0.2191 ± 0.0008 and e_B = 0.2257 ± 0.0009, respectively). As BU CMi has the smallest P_{A, B}/P_out ratio of these three quadruples and, moreover, contains the closest inner binaries amongst these three systems, these findings at first sight appear to be quite unexpected for a number of reasons. First, the less tight a hierarchical system is, the less effective are the gravitational perturbations of the outer component(s) acting upon the Keplerian motion(s) of the inner pair(s). Second, due to the compactness of the inner binaries, the fractional radii of the constituent stars in BU CMi exceed the same quantities of the stars in the other two quadruples by factors of 4–5. Hence, one can expect that tidal forces and tidal dissipation, which are related to the fifth (for equilibrium tides) and eighth order (for tidal dissipation) of the fractional radii are much more effective in BU CMi than in the other two systems.

These discrepancies can be resolved by considering the facts that (i) in contrast to TICs   454140642 and 219006972, BU CMi is formed by four early-type, massive, radiative stars, in which tidal dissipation is much less effective; and moreover, (ii) this system is much younger (i.e. ∼300-Myr-old, in contrast to the other two quadruples which are several Gyr old). Thus, one may conclude that there has been insufficient time for the circularizaton of the inner orbits since the formation of the quadruple system BU CMi.

The significant inner eccentricities in BU CMi, as well as the non-edge-on view of the two inner orbital planes (i_A = 83.4° ± 0.1°; i_B = 83.9° ± 0.1°), together with the very similar surface brightnesses of the binary stars, have an interesting observational consequence. Namely, depending on the orientations of the orbital ellipses relative to the Earth, the depths of the two consecutive eclipses in each EB alter each others. Regarding binary A, currently the periastron passage, i.e. the smallest separation between the two stars, is much closer to the inferior conjunction of the more massive (primary) component. Hence, due to the larger fraction of the occulted stellar disc of the secondary, the deeper light minimum occurs during this event, i.e. when the more massive, hotter star eclipses the less massive and cooler component. In contrast to this, during the HAT observations, nearly half an apsidal cycle earlier, the currently shallower eclipse was actually the deeper eclipse. In the case of binary B, the tendency is just the opposite, as can be seen in Fig. 7.⁵

Turning to the apsidal motion cycles mentioned above, our direct spectral fits give a ∼25 yr-period for both inner binaries. Allowing for the fact that the duration of the full spectroscopic data set is ∼2.7 yr, i.e. about 10 per cent of the full cycles, we feel that these results are in essentially perfect agreement with the findings of the spectro-photodynamical analysis (P_{apseA, B} = 28.7 ± 0.1 and 25.1 ± 0.1 yr, respectively). The latter results are based on the ETV data, which extend to more than the half of the full apsidal cycles. In addition to the apsidal motion periods of the inner binaries, the photodynamical solution also gives the apsidal motion for the outer orbit, and this is also found to be remarkably rapid at P_apseAB = 145.3 ± 0.2 yr. These values, which were ‘measured’ from the numerical integration of the best-fitting photodynamical solution, are in good accord with the theoretically calculated apsidal motion periods – of which the medians were found to be P_apse = 31.8 ± 0.1, 29.0 ± 0.1, and 151.6 ± 0.8 yr for orbits A, B, and AB, respectively.⁶ Besides the ‘measured’ and theoretically calculated apsidal motion periods, we also tabulate the contribution of the dynamical (third-body), relativistic, and classic tidal effects to the apsidal advance rates (Δω_{3b, GR, tide}, respectively). It is evident that the dynamical effects substantially dominate over the tidal and relativistic ones, and this clearly supports our previous statement in Section 4.2 that the rapid apsidal motion must have a dynamical origin.

This latter statement leads us back to the previously mentioned contradiction between our results and those of Volkov et al. (2021) which will be discussed below.

8.2 Comparison with the results of Volkov et al. (2021)

As mentioned in the Introduction, Volkov et al. (2021) also published a detailed analysis on this quadruple system. Their conclusions are remarkably different from our findings. Here, we discuss the origins of these discrepancies and make attempts to resolve them.

First, the most substantial difference is that they report an outer period of about 6.6 yr (2420 ± 20 d), instead of the much shorter outer period of 121 d that we report in this work. They arrived at this conclusion via an analysis of the ETV curves, where besides the evident large-amplitude sinusoidal variations characteristic of dynamically driven apsidal motion with a ∼25 yr period, they detected additional cyclic variations with a period of ∼2400 d. They interpret the latter variations as the LTTE originating from the orbit of the two EBs around their common centre of mass. Moreover, they also detect another variation of small amplitude in the high-precision TESS Sector 7 eclipse times, which they claim might be due to a small ∼60-d-period libration (which they call ‘nutation’) in the lines of the apsides of both binaries, caused by the other pair. They do not investigate, however, how such a relatively wide outer orbital separation would be able to produce such a short-period effect, nor how it could account for the rapid apsidal motion.

In contrast to their results, we found that these small-amplitude, so-called ‘nutations’ are the principal third-body effects. These variations are, in fact, the same thing as we have seen in the ETVs of several tight, hierarchical triple stellar systems where the ETVs are dominated by the P_out time-scale third-body perturbations due to the tertiary (which, in the present case, is also a binary itself). These variations can be modelled both analytically (see e.g. Borkovits et al. 2015, 2016) and numerically, as in this work. Our photodynamical analysis shows clearly that these structures in the ETVs of both EB pairs can be explained with the mutual P_out time-scale third-body perturbations due to the other EB, and consequently, the two binaries orbit around each other with a period of P_out = 121.50 ± 0.02 d (see the lower panels of Fig. 3). This conclusion makes this 2+2 quadruple momentarily the one with the shortest known outer period.

Regarding the longest period – the largest amplitude cycles on the ETVs (of ∼25 yr) of both binaries – we agree with Volkov et al. (2021) that these arise from apsidal motion. We note, however, that this relatively short period of the apsidal motions in both binaries can be explained neither by combinations of the classic tidal effects and the general relativistic contribution in the EBs, nor by the gravitational perturbations of a relatively distant companion in an ≈6 yr-period. The characteristic time-scale of dynamically driven apsidal motion is proportional to |$P_\mathrm{apse}\propto P_\mathrm{out}^2/P_\mathrm{in}$| which, in the case of P_out = 2420 d yields ≈4923 and ≈5457 yr for binaries A and B, respectively. In contrast to this, as one can see e.g. in the upper panels of Fig. 3, and also in Table 7, where we tabulated the tidal, relativistic and third-body contributions to the apsidal advance rate separately, our photodynamical solution indeed produces the observed apsidal motions.

The most pressing issue regarding the ETV curves is the origin of the extra cyclic variations with a period of a few years. As was mentioned above, Volkov et al. (2021) find their period to be ∼2420 ± 40 d, while our complex photodynamical analysis (and the preliminary analytic ETV studies, as well) have resulted in a period close to half the value they found. Next, we attempt to explain the origin of these ‘extra’ ETV variations that occur on a time-scale of ∼1200 d.

8.3 Origin of the ETV variations on a time-scale between P_out and P_apse-node

As mentioned above, the most interesting question about BU CMi is the origin of the extra, ≈3 or 6 yr-period cyclic variations of the ETVs. This effect, at first sight, looks quite mysterious because these cycles appeared naturally in the numerically generated ETVs during our modelling processes (see e.g. the upper panels of Fig. 3). However, they do not appear even in our most detailed analytical ETV model, described in Borkovits et al. (2015). In other words, they are not present in our formulae based on the once-averaged and doubly averaged octupole-order analytic perturbation theories of hierarchical triple stellar systems.⁷

In order to identify the origin of these effects, we checked the variations of the instantaneous osculating orbital elements during each numerical integration step. In Fig. 8, we show that we found periodic variations in the semimajor axes of the two binaries’ orbits that are similar to what we see in the observed ETV curves. We display in the upper right panel of Fig. 8 that these variations have exactly the same periods and opposite phases in the two binaries, and hence, their effect on the ETVs may really mimic the signs of mutual LTTE. The connection between these small variations of the semimajor axes and the ETVs are also well illustrated in the lower left panel of Fig. 8. We can easily infer this connection not only for the similar periods, but also for the fact that the larger the amplitude of the variations in the semimajor axes, the larger the amplitude of the extra ETV bumps. And, moreover, one can see that there is an exact 90°-phase shift in between the extrema of the semimajor axis and ETV variations. This is in perfect accordance with the fact that the ETV represent some cumulative or, integrated variations. Or, mathematically, a sine-like perturbation in the semimajor axis (or, in the mean motion) will result in a cosine-like ETV for the additional integration. This finding is not surprising, insofar as these slight variations in the semimajor axes naturally reflect the instantaneous mean motion of the given body, to which the variations of the ETVs are extremely sensitives. On the other hand, such kinds of variations in the semimajor axes are surprising in the sense that both the once and doubly averaged perturbation theories of hierarchical triple systems arrive at the conclusions that there are no medium-, long-period, and secular perturbations in the semimajor axes. In other words, according to the perturbation theories that are usually used in the hierarchical stellar three-body problem, one should not expect any variations in the semimajor axes for which the periods exceed the period of the inner binary. Note, however, that in our last theoretical work about the analytical description of the ETVs of tight triples (Borkovits et al. 2015), we introduced some further terms, which allow periodic perturbations in the inner semimajor axes on the time-scale of the outer period; here, however, we detected periodic perturbations in these elements with a factor of ∼10 longer period.

In order to test the physical origin of this behaviour, we made some additional numerical runs where all but one of the initial parameters were the same (i.e. they were taken from the best-fitting model), and only the outer period of the quadruple was slightly modified. As one can see in the upper panels of Fig. 8, only a small change in the outer period can result in substantial variations in the amplitude and the period of the effect, and this can also be seen clearly in the corresponding simulated ETVs (bottom panels of Fig. 8).

In generating and studying Fig. 8, we found it interesting to note that the closer the ratio of the outer to inner periods is to an integer, the larger the amplitude and longer the period are of these few-year-long cycles. For example, in the order of decreasing amplitudes, in binary A the period ratios are 40.99; 41.33; 41.73; while for binary B: 36.94; 37.24; 37.61.⁸ In our view, these facts suggest some similarities with the mean-motion resonances of the classical planetary perturbation theories. Namely, when the mean motions are nearly commensurable this may lead to large-amplitude and long-period perturbations in the given orbital elements. Normally, however, only low-order mean motion resonances will produce large amplitude perturbations, because the amplitudes of the resonant perturbations in general are multiplied by some power-law functions of the eccentricities, in which the powers are proportional to the order of the resonances. Here, we suspect some similar behaviour: Namely, due to the commensurability of the inner and outer periods, some very high order ‘resonant’ terms yield some contribution which cannot be averaged out perfectly. This actually results in a very low-amplitude cyclic variation in the semimajor axes (the relative variations are of the order of 10⁻⁴–10⁻⁵). But, due to the extreme sensitivity of the ETVs to the mean-motions, these tiny variations may produce observable effects in the occurrence times of the eclipses.

We note that this effect is worthy of a more detailed and quantitative investigation. But it should be carried out first for a simpler configuration, i.e. for an actual 2+1 triple system instead of a 2+2 quadruple. In this regard we note that the photodynamical analysis of the tight triple TIC 167 692 429 (Borkovits et al. 2020), which was based on the first year of TESS data revealed similar extra cycles in the ETVs; these have now been verified by the third and fifth years of new TESS observations. We will further investigate these extra cycles in the context of this triple system in the near future.

9 CONCLUSIONS

The bright variable star BU CMi is composed of two short-period EBs. Our work has found it to be the tightest quadruple system known, with P_out ≃ 121.5 d. Although the quadruple nature of the system was first pointed out by Volkov et al. (2021) from the MASCARA photometry, they found a much longer outer orbital period of P = 6.62 yr, attributing the 121-d ETV variations to a ‘nutation’ effect. Simultaneous modelling of high-dispersion spectroscopy and TESS light curves conclusively demonstrates that the line profile variations cannot be explained by the 6.62-yr periodicity. The spectroscopic data clearly show large, ∼50 km s⁻¹, variations in the systemic velocity of the inner subsystems (i.e. the two EBs). For P_out = 6.62 yr, the amplitude of the RV changes of the binary centres of mass would be only about ∼18 km s⁻¹.

Even more conclusive proof of the extremely tight outer orbit comes from the detailed photodynamical modelling, which takes into account satellite photometry, all available eclipse timing data, and the RVs. The numerical integration of the orbits explains the ETV changes with the 121-d period (the outer orbit), establishes the longer-term ETV variability on the ∼1200-d period, and perfectly predicts the rapid apsidal motions in the EBs on ∼25-yr time-scales. This apsidal motion is very rapid, and it is dynamically driven by the mutual gravitational perturbation of the binaries. For the much longer outer orbit determined by Volkov et al. (2021), the apsidal motion rate would be substantially slower.

Although the orbital and stellar component parameters are well constrained by the complex photodynamical modelling, there remain a few open questions that require further observations. For example, we seek to constrain any spin-orbit misalignment; this possibility, however, appears unlikely due to the near-coplanarity of the inner and outer orbits. Determination of the spin-orbit misalignment would require dedicated spectroscopy of the system during the eclipses of the stellar components. The small observed range of the RVs of the components and their high rotation rates would complicate the modelling of the Rossiter–McLaughlin effect.

Assuming that the stellar rotational rates are quasi-synchronous (rotation rate equals the Keplerian orbital rate at periastron) and that the spin axes are perpendicular to the inner orbits, we can further constrain the stellar radii from their measured projected rotational velocities vsin i.

Although both of the inner binaries are eclipsing, the relatively low inclination of the outer orbit, i_{A − B} ∼ 83.8°, and a relatively large semimajor axis, a_{A − B} ∼ 221 R_⊙, preclude outer-orbit eclipses. A shallow outer eclipse could occur if the outer inclination angle were larger than about 88.5°.

Finally, we note that the BU CMi system could have had an outer period as short as ≈32 d and still be dynamically stable (see e.g. Mardling & Aarseth 2001; Mikkola 2008; Rappaport et al. 2013). This assumes the same EB periods and masses, the same outer eccentricity, and the same orbital coplanarity. Thus, there is much room in phase space to find even tighter quadruples. Whether evolution scenarios will permit such short-period quadruples is another matter. It is therefore worth trying to observationally push the boundaries to ever shorter outer periods.

ACKNOWLEDGEMENTS

AP acknowledges the financial support of the Hungarian National Research, Development and Innovation Office – NKFIH Grant K-138962. TP and RK acknowledge support from the Slovak Research and Development Agency – contract No. APVV-20-0148 and the VEGA grant of the Slovak Academy of Sciences No. 2/0031/22. GB, ZC, and JH acknowledge funding from NASA grant 80NSSC22K0315. We would also like to thank the Pierre Auger Collaboration for the use of its facilities. The operation of the robotic telescope FRAM is supported by the grant of the Ministry of Education of the Czech Republic LM2023032. The data calibration and analysis related to the FRAM telescope is supported by the Ministry of Education of the Czech Republic MSMT-CR LTT18004, MSMT/EU funds CZ.02.1.01/0.0/0.0/16|$\_$|013/0001402, CZ.02.1.01/0.0/0.0/18|$\_$|046/0016010, and CZ.02.1.01/0.0/0.0/18|$\_$|046/0016007.

This paper includes data collected by the TESS mission. Funding for TESS is provided by NASA’s Science Mission Directorate.

DATA AVAILABILITY

All photometric and spectroscopic data used in this paper and the codes used for the direct fitting of the spectra and photodynamical analysis will be shared upon a reasonable request to the corresponding author.

Footnotes

References