© The Authors 2025

1. Introduction

Classical eclipsing binaries are sometimes presented as the foundations of modern stellar astrophysics. This is due to their role in 20th-century astrophysical research. Eclipsing binaries (hereafter EBs) are routinely used to precisely measure stellar radii and masses, as well as to calibrate the distance ladder (nowadays even outside of our Galaxy), and stellar evolution theories via looking into stellar interiors (for a discussion of all of these aspects, see e.g. Southworth 2012; Claret & Giménez 2010).

In the current era of large-scale surveys and ultra-precise satellite data, it is not possible to perform a detailed analysis of every single EB in the sky, since their number now reaches more than 2 million (Mowlavi et al. 2023). Therefore, it is essential to focus detailed research on some very specific types of EBs that are limited in number and that have high astrophysical importance, i.e. they offer us some unique insight into otherwise hard-to-detect processes and effects. By this we mean objects for which there are enough observational data for analysis, but that are still rare enough in the sky.

We are convinced that doubly eclipsing systems of the 2+2 architecture are clearly such objects. In total, among the known more than 2 million EBs, there are currently only about 1000 known candidates for doubly eclipsing systems. By a doubly eclipsing system candidate, we mean one point source in the sky, for which two different eclipsing periods are detectable in the data. We still mark them as candidates only, since we do not know anything about their true connection and gravitational coupling. Only after a proof of their mutual movement around a common barycenter is presented can we definitely say that these are real 2+2 quadruples.

However, it seems like the rate of optical doubles and blends of two close-by stars with eclipsing periods incorrectly attributed to only one star (i.e. false detections) is relatively low. Among the thousand known published doubly eclipsing systems, there are only several dozen proven optical doubles and blends (see e.g. Fezenko et al. 2022 and Zasche 2024). Therefore, one can ask why there are only about 60 definitely confirmed 2+2 quadruples among the order-of-magnitude-larger set of candidates with two known periods. We can see two different explanations for this apparent discrepancy. The first one is that the detection of the mutual orbit of the two doubles around a common barycenter is complicated and not always straightforward, and needs much more effort than the pure detection of the two eclipsing periods in the photometric data. The second reason could be the fact that we are still dealing with quite limited data. Mutual periods can be decades, centuries, or even thousands of years long, and our available time span of data is usually much shorter. Moreover, the orientation of the long orbit plays a role in the successful detection of the orbit. On the other hand, the 2+2 quadruples are obviously more common in the stellar population compared to the 3+1 quadruples, and their parameters can give us useful hints concerning their formation mechanisms (see e.g. Tokovinin 2021).

Our presented work is a natural continuation of our long-term effort to detect and study these 2+2 quadruples. Without our contribution, the whole sample of the doubly eclipsing candidates would be 17% smaller. Even more significant is the result of our previous publications on the number of proved 2+2 quadruples with derived mutual orbits. The total number of such systems would be about 65% lower without our contribution. Therefore, we feel that the long-term effort of collecting the ground-based data and analysing them in combination with the satellite ones is still useful for the topic and brings us a non-negligible contribution of new systems.

2. Target selection

The targets selected for the current study were chosen from the set of candidates among the published systems. Most of them are from our recent paper on TESS discoveries (Zasche et al. 2022). One was taken from an older paper dealing with the Corot data (Hajdu et al. 2017). The selected stars are a mixture of northern and southern targets. Their spectral types or masses also did not play a role in our selection process. The only criterion used was the detection of two eclipsing periods in the available data and the non-existence of any detailed analysis of the star. From this point of view, our present paper can be considered as a first confirmation of these stars as proven quadruples, i.e. with their mutual orbit firstly derived. A summary of information about the six targets is presented in Table 1. One can see there that the stars are relatively faint targets, and we have only very limited information in our hands at present (i.e. no detailed spectroscopy).

The detection of the systems as doubly eclipsing was originally done by (Zasche et al. 2022), who easily plotted the TESS data where the additional eclipses are directly visible by the naked eye. However, the derivation of the proper period of the secondary pair was done by subtracting the dominant eclipsing pair, and on the residuals the period was derived more robustly over several TESS sectors of data.

The selected systems analysed in our study seem to be only a drop in the ocean, and not exceptional among the other 2+2 quadruples that have already been proven. A comparison with the other systems is summarised below in Section 5. What can be stated already is the fact that we are still dealing with quite a low number of contact binaries as the inner pairs among these doubly eclipsing systems (see Figure 2 in our recent paper Zasche et al. 2022). Only about one quarter of the currently proved doubly eclipsing 2+2 quadruples contain at least one inner binary showing a contact-shaped light curve (hereafter LC). Even in our present paper, only one system (the first one, CzeV1254) shows such a LC. However, this deficiency is still probably just a consequence of the problematic analysis of these contact LCs and their disentangling, rather than an indication that there is a lack of such systems in the stellar population. As was already stated above, our contribution to the topic with our six new systems represents 10% of the total number of known 2+2 doubly eclipsing quadruples. With hundreds of them in the future, we can hope to find some clues concerning their formation mechanism, for example (see e.g. Tokovinin 2021 discussing different formation channels between 2+2 and 3+1 quadruples).

3. Data handling and fitting

The TESS data (Ricker et al. 2015) were the most important source of photometric information for our systems. This is due to their accuracy as well as their cadence. We extracted the photometry from the TESS archive using the software LIGHTKURVE (Lightkurve Collaboration 2018). The TESS data were used for the analysis of LCs of both inner binaries (except for the last system).

For the fitting of the LCs, the software PHOEBE (Prša & Zwitter 2005) was used. It uses the standard Wilson-Devinney algorithm (Wilson & Devinney 1971) with its later modifications. When dealing only with photometry, several important input parameters have to be fixed or estimated elsewhere.

The orbital period of a binary has to be derived outside of PHOEBE. This preliminary period estimation was done on the best TESS sector of data. After that, the period was step-by-step fixed for some short time interval, for example for one TESS sector, as it is apparently changing continuously as the binary revolves around a common barycenter with the other binary in the system. Such an approach is well substantiated due to the fact that the duration of one TESS sector of data (i.e. 27 days) is always an order of magnitude shorter than the period variation on the long orbit. Therefore, it produces errors in the times of minima at the level of 10⁻⁵ to 10⁻⁴ day at maximum, i.e. an order of magnitude smaller than typical uncertainties of the times of minima derivation precision. For the mass ratio of the particular binary, we started with the assumption of equal masses, i.e. q = 1.0, but later this parameter was also released and fitted for systems with larger out-of-eclipse variations.

For the whole analysis, we followed our step-by-step iterative procedure. The same method was already used for our previous studies on the topic (see e.g. Zasche et al. 2024, 2023). It means that we first performed a preliminary fit of the most pronounced pair, A. Then this fit was subtracted from the complete (A+B) photometry and the resulting pair, B, was fitted as well. This preliminary fit of pair B was again subtracted from the original combined photometry and only pair A was fitted again. This procedure was repeated several times until we obtained the final acceptable results for both pairs.

Together with the fitting of the LCs, the periods of both pairs were computed with a higher precision as well. This step was also done iteratively several times. We used our method of AFP (automatic fitting procedure, introduced in Zasche et al. 2014), which deals with the LC template and provides us with the times of eclipses. From the times of eclipses, the better period for the particular pair can be derived. With such an improved orbital period, a better LC fit can be obtained. We implemented this step-by-step procedure for both inner eclipsing pairs, A and B, until we got the final results for both pairs. For this AFP method, the orbital period was also fixed at its actual value (a linear fit from the long-term period variations) for particularly short intervals of data, i.e. hundred(s) of days typically.

We also have to be aware of the fact that the TESS data may in principle have different levels of light contamination in different sectors. This aspect of the whole TESS photometry was taken into account when fitting the LCs of both pairs. This means that the third light values of both A and B pairs were fitted independently as well. Therefore, we consider the third light value to be a free parameter for both pairs, trying to estimate the light contributions of both A and B pairs. However, due to varying values in between different sectors, this can only give us a rough estimate. The values presented below in Table 2 give only one value for the TESS sector with the best data (the one with the shortest cadence of data). In principle, this value can be higher than 100% for pairs A and B in total, which can indicate additional light from some close-by components that are not necessarily connected with the system itself, for example. We assumed that the value of the third light is constant (i.e. comes from a non-variable star) during at least one TESS sector of data, for instance.

The proper orbit of the two doubles around a common barycenter (and hence the proof of their 2+2 quadruple nature) can be derived from the long-term ETVs (eclipse timing variations). Both ETV curves of the A and B pairs should behave in the opposite manner.

Both the older (and less precise) photometry and the new dedicated ground-based observations from our observers can help us greatly in this regard, since the TESS photometry itself covers only a small part of the long orbit. Including these observations, our time baselines improved to much longer time intervals of about 20 years in some cases. The new observations were obtained mainly on these observation sites:

• The Danish 1.54-m telescope on La Silla, Chile, remotely operated, equipped with a CCD camera, and standard filters,

• The 65-cm telescope at Ondřejov observatory in Czech Republic, using a MII G2-3200 CCD camera,

• The 30-cm telescope at a private observatory in Veltěže u Loun, Czech Republic, equipped with a MII G2-8300 CCD camera,

• The 20-cm telescope at a private observatory in Dolní Bousov, Czech Republic,

• Two FRAM telescopes of 30 and 25 cm, located on La Palma (Prouza et al. 2019), and in Argentina (Aab et al. 2021),

• The 30-cm ROST telescope of the Variable Star and Exoplanet Section of the Czech Astronomical Society located at Ondřejov observatory in the Czech Republic, equipped with a MII G4-16000 CCD camera.

For the analysis of the period variations of both pairs, we used a standard least-squares method of fitting the data with the theoretical curve. By this curve we mean a standard LITE (light-time effect, or also named Römer delay) curve describing the orbit around a barycenter with the distant body (see e.g. Mayer 1990, or Borkovits et al. 2016):

$$\begin{array}{c}\hfill {\mathrm{ETV}}_{\mathrm{LITE}}=\frac{A_{\mathrm{LITE}}}{\sqrt{1-e^2{cos}^2\omega }}\left[\frac{(1-e^2)·sin(\nu +\omega )}{1+ecos\nu }\right].\end{array}$$(1)

However, for some tighter systems where the outer orbit of the two binaries is close enough to the inner two orbits (i.e. for orbits with periods of p_A − B ≲ 1yr) one also has to take into account the dynamical interaction between the orbits. We used an approach that has been presented for this so-called ‘physical delay’, for example by Rappaport et al. (2013):

$$\begin{array}{c}\hfill ET{V}_{\mathrm{phys}}=A_{\mathrm{phys}}((2·{cos}^2{i}_m-2/3)·(\nu +esin\nu -M)+\\ \hfill (1-{cos}^2{i}_m)·(sin(2\nu -2v_m)+e·sin(\nu -2v_m)+\frac{e}{3}sin(3\nu -2v_m))).\end{array}$$(2)

This means that for the more compact systems, both of these effects should be taken into account together and added. In these two equations, (1) and (2), there are the following symbols: A stands for the amplitudes of both effects, e is the eccentricity of the outer orbit and ω its argument of periastron, and ν is the true anomaly on the long orbit, while M is the mean anomaly, v_m stands for the orientation of the periapse of the outer orbit to the inner binary orbit, and i_m is the mutual inclination of the inner binary and the outer orbit. The amplitudes of both these effects depend on the semi-major axes of the orbits, or the masses of the individual components. At least in principle, all these computable quantities should be derived from good data (which means data that are well distributed in time and that have an adequate accuracy).

4. Individual systems

In the following sections, we focus on each system from our sample one by one, emphasising their differences and individual aspects of the analysis.

4.1. CzeV1254

The very first system is CzeV1254 (= UCAC4 732-033136, V = 13.3 mag out of eclipse). The star was first discovered as a classical EB. Later, in Zasche et al. (2022), it was also classified as a doubly eclipsing candidate due to the discovery of two eclipsing periods (of EA+EW type): 0.715419 + 0.362593 days. Both EBs show rather deep eclipses of about 0.1 mag; however, no detailed analysis of the system has been published to date. According to the GAIA DR3 temperature estimate, the system is probably of G-dwarf type, being about 565 pc distant from the Sun (GAIA DR3).

Of all the available TESS sectors, we selected sector 59 for further analysis. The two overlapping photometric periods are clearly seen there. We were able to separate the two signals using the methods described in Section 3. Despite the fact that both pairs have rather short periods, we were able to disentangle them and perform the independent LC analysis for both pairs. The PHOEBE fitting led to the LC parameters given in Table 2 and the fit presented in Fig. 1. One can see there that the light contribution of both pairs is roughly the same. It means both binaries should also be of a similar mass (using the standard and easiest assumption of the main-sequence mass-luminosity relation). The slight asymmetry of a pair A light curve is probably caused by some photospheric spots, but we have not tried to fit this.

Fig. 1. Analysis of the CzeV1254 system. Left-hand side: Light curves of both disentangled pairs, A and B, based on the TESS data (residuals given in the bottom plots). Right figures: ETV diagrams of both pairs as they revolve around a common barycenter. Full dots stand for primary eclipses, and open circles for the secondary ones. Green symbols represent the TESS data, and red ones our new dedicated observations. The remaining symbols are the other photometric surveys.

For the study of period changes in both pairs, long-term evolution, and the detection of mutual movement of both pairs around a barycenter, we collected all available data for the star. This meant that apart from the new TESS data and our new dedicated observations, we also used older ASAS-SN data (Jayasinghe et al. 2018; Kochanek et al. 2017), SWASP (Pollacco et al. 2006), and ZTF (Masci et al. 2019). Thanks to these data, the interval covered with data spread to 20 years. Therefore, our fit presented in Fig. 1 with its periodicity of about 15.6 years (see Table 3) is now covered. It is still possible that, due to rapid period change near the periastron passage, the period is a bit longer. Therefore, new observations in the upcoming seasons would still be useful to obtain. A slight asymmetry of the LC as is presented in the residuals in Fig. 1 could also cause an improper derivation of the time of the mid-eclipse, or even additional apparent period variation in the ETV diagrams (see e.g. our study on the topic published in Zasche 2011). However, in this case, only a slightly worse time of minimum derivation is caused, shown by a larger scatter in the ETV diagrams. The overall picture of our fit remains practically the same, since the fit is derived through the opposite behaviour of A and B pairs. A similar evolution of asymmetries for both pairs in time is highly improbable.

The comparable amplitudes of the A and B pairs indicate similar masses of both pairs, which is in agreement with the third light ratios derived from the LC analysis as presented in Table 2. Due to its relatively small distance from the Sun, the value of the predicted angular separation of the double in the sky of about 30 mas is quite favourable. This is well within the limits for modern interferometers. However, more problematic is its relatively low brightness for such a technique.

4.2. ASASSN-V J070838.27-171952.9

The second system analysed was ASASSN-V J070838.27-171952.9 (= TIC 148942773, V = 13.41 mag out of eclipse). This eclipsing system in the Canis Majoris constellation was first discovered as an EB by the ASAS-SN survey, which detected its prominent 4.3-day eclipses. However, with precise TESS data we were able to detect a secondary periodicity of about 3.77 days. Both of these eclipse-shaped curves indicate well-detached systems for both A and B binaries. No other detailed analysis exists and no spectroscopy is available for this star.

We followed more or less the same approach as for the first system. The TESS data were analysed and the LC fit was derived (using TESS sector 33). The parameters of this PHOEBE fitting are given in Table 2. Despite the detached configuration, both inner pairs are circular (the small eccentricity value of pair A is still rather uncertain, see the residuals in Fig. 2). Due to only very small out-of-eclipse variation, we did not try to fit the mass ratio of either of the pairs and let this parameter stay fixed to its suggested value of equal mass, q = 1.0. Pair B seems to be the less luminous one (as was expected from a much shallower eclipse depth). The LC fits of pairs A and B are plotted in Fig. 2.

Fig. 2. Analysis of the ASASSN-V J070838.27-171952.9 system. The description of the plots and colour coding for individual symbols are the same as in Figure 1.

Collecting all older photometric data, we tried to detect the long-term period modulation of both pairs. We used the following sources for this system: ASAS-SN, ZTF, Atlas (Heinze et al. 2018), and our own data. The result of our fitting is plotted in Fig. 2, and the results are written in Table 3. Due to the rather poor quality of the data and the only very small amplitude of ETV variations in both pairs, we used only a simplified approach of a circular outer orbit with zero eccentricity. A much larger number of more precise observations would be needed for a better characterisation of this mutual orbit.

4.3. ASASSN-V J091951.17-593306.9

The third star analysed is ASASSN-V J091951.17-593306.9 (= TIC 386900081, V = 13.83 mag out of eclipse). This star was originally also classified as a standard EB with its orbital period of about 2.669 days, based on the ASAS-SN data (Jayasinghe et al. 2020). However, thanks to new precise TESS data we were able (Zasche et al. 2022) to detect an additional secondary eclipsing period of about 14 days. No other detailed analysis for this star was published to date.

As for the previous cases, we used the TESS photometry for detailed modelling of both inner pairs, A and B. The results are given in Table 2, while the fit is plotted in Fig. 3. It is evident that both pairs A and B are well detached. Both pairs have a similar brightness, as can be seen from the similar values of the third lights for both fittings. Pair B has a rather eccentric orbit and apsidal motion (e = 0.17, U = 37 yr). There is still a possibility that pair A also has a very small eccentricity, but only at a level of about 0.001. This is still barely detectable, even with precise TESS data.

Fig. 3. Analysis of the ASASSN-V J091951.17-593306.9 system. In the right-hand side figures are the ETV variations plotted twice for pair B. The upper plot is the complete ETV diagram (i.e. apsidal motion + Römer delay + physical delay). The lower figure is plotted after subtraction of the apsidal motion. The blue curve is the final combined fit, the red one stands for the physical delay, and the magenta is the Römer delay (LITE).

For this system, we had to use a combined approach of fitting LITE (Römer delay) together with the physical delay (i.e. both equations 1 + 2 above) due to the very short mutual orbit. The two doubles revolve around a common barycenter with a periodicity of only about 347 days (the ratio of the inner period, P_B, over P_AB is only 24). With such an approach, instead of deriving five LITE parameters (p_AB, e, ω, A, and T₀) we have to also include another three (i_m, v_m, and M_A/M_B) in our computations. This means that the dynamical interaction between the orbits is not negligible at all. One can even expect some interaction of the orbits, such as the orbital precession of either of the orbits in the case of a non-coplanar configuration of the system (see e.g. Borkovits et al. 2022). This should lead to eclipse depth variations over a longer timescale.

4.4. WISE J100820.07-731554.0

The star WISE J100820.07-731554.0 (= TIC 453459782, V = 14.9 mag out of eclipse) is the fourth system studied in our sample, and one of the faintest. It consists of two quite dissimilar binaries. The more prominent one is a contact-like binary with an orbital period of about 0.37 days, discovered as an EB even in WISE data (Chen et al. 2018). The other EB is well detached, with a period of about 7.9 days. Moreover, the system was also resolved as a double by GAIA (Gaia Collaboration 2023). The two components are about 0.8″ distant and are probably bound together, since both share a similar parallax and proper motion.

A similar approach as in previous cases was also used here. The LC fitting was based on the precise TESS data (from sector 65). The results of this fitting are presented in Table 2. The final fit of the inner pairs, A and B, is plotted in Figure 4. One can see that both pairs have a quite comparable brightness (and hence probably also masses), pair A being perhaps slightly more prominent. There is a slight asymmetry in the LC of pair A that is especially visible in regions outside of eclipses, probably due to spots on the surface of either of the stars.

Fig. 4. Analysis of the WISE J100820.07-731554.0 system. The description of the plots and colour coding for individual symbols are the same as in Figure 1. The ETV diagram of pair B was plotted after subtraction of the slow apsidal motion.

Our ETV analysis was performed with the collection of all available data. However, the variation is only very small, at the edge of detectability, even for the precise TESS data. Pair B shows very long apsidal motion in its eccentric orbit (e = 0.14, U > 1000 yr). The results of our fitting are plotted in Figure 4, and the parameters as well as ephemerides for both pair A and pair B are given in Table 3. Due to the low quality of the data and low amplitude of the ETV variation, we only used a simplified approach of a circular A-B orbit. The periodicity is barely visible, and needs to be confirmed with new dedicated observations of reasonable quality in the coming years, if possible. The predicted 7 mas separation of the two A&B binaries that results from our analysis indicates that this pair cannot be easily attributed to the detected components from GAIA (about 100× larger separation). This gives us a weak indication of its quintuple nature.

4.5. ASASSN-V J143536.01-721459.4

The next system, ASASSN-V J143536.01-721459.4 (= TIC 401924368, V = 12.71 mag out of eclipse), is the brightest star in our sample. It was originally detected as an EB in the ASAS-SN data (Jayasinghe et al. 2018) and later classified as doubly eclipsing by Zasche et al. (2022). Both of its inner binaries are detached with a rather long orbital periods. Pair A shows quite deep eclipses (about 0.3 mag), while pair B is about ten times shallower.

No detailed analysis for this system exists; therefore, we analysed the TESS data and performed the LC fitting. Using sector 65 we arrived at the LC parameters given in Table 2. The final plots for both inner pairs are given in Fig. 5. The two pairs have rather different brightnesses, with pair A dominating significantly. It is therefore not surprising that the ETV analysis presented in Figure 5 also shows that pair B has a higher amplitude of the ETV variation. It is still not a very robust solution due to the poor quality of the eclipsing times data for pair B, and its amplitude could be even higher. This is another system for which new data in the coming years would be beneficial to finally confirm our hypothesis with a higher conclusiveness.

Fig. 5. Analysis of the ASASSN-V J143536.01-721459.4 system. The description of the plots and colour coding for individual symbols are the same as in Figure 1.

4.6. Corot 310284765

The last system in our set of analysed stars is also the faintest one, named Corot 310284765 (= TIC 44874302, V = 15.28 mag out of eclipse). The doubly eclipsing nature was first detected by Hajdu et al. (2017) using the Corot data. However, since then no detailed analysis of the star has been published to date. According to its spectral classification (see Table 1) this is the system with the earliest spectral type; hence, it also has the highest temperature in our sample. According to GAIA DR3, there are also some close-by components about 2.7″ and 4.3″ away from Corot 310284765. Nonetheless, it is unlikely that these objects are gravitationally bound to the star, as is evidenced by their distinct parallaxes and proper motions.

This is the only system from our sample for which we used the Corot photometric data for LC modelling instead of TESS. This is due to their much better quality (the telescope has a significantly larger aperture, and hence the data obviously have a much lower scatter). Owing to the better angular resolution of the Corot data (about 10× better compared to TESS), also the two EBs appear to be deeper (approximately 2×) in the Corot data compared to the TESS data. The results are given again in our Figure 6, where one can see the fits for both pairs, and the final parameters are given in Table 2.

Fig. 6. Analysis of the Corot 310284765 system. The description of the plots and colour coding for individual symbols are the same as in Figure 1. The Corot data are shown as magenta symbols.

Both pair A and pair B show slightly eccentric orbits, which is clearly seen in the LCs. This can also be seen in our ETV analysis (plotted in Figure 6), in which the slow apsidal motion is evident for both pairs. Apart from that, periodic variation in the pairs as a consequence of movement around a common barycenter is also visible. The solution is still rather preliminary, since its coverage is poor, and new observations can still shift the solution significantly. Nevertheless, the bound nature is proven with our presented solution.

5. Conclusions

The doubly eclipsing quadruple systems still present a challenge when it comes to modelling their observational data. These are usually only photometric observations, because of their lower brightness. To date, only a select number of doubly eclipsing systems have been studied in detail, combining both photometry and spectroscopy (see e.g. Kostov et al. 2021, or Pribulla et al. 2023). These systems provide us with a unique insight into the star formation and consequent orbital evolution due to the possible derivation of overall architecture of the whole quadruple, derivation of mass distribution, and mutual inclination of the orbits, for example.

Quite recently (Zasche et al. 2023), we also discussed the possibility of mean motion resonances in the quadruple system (i.e. the A and B orbital periods) and the spectral type of the component stars. However, here in our sample the individual systems seem to be rather distant from the resonances. The closest ones are CzeV1254, which is about 1.4% distant from the 2:1 resonance, and Corot 310284765 about 1.1% distant from the 4:5 resonance. However, both values are too far away from any potential resonance configuration.

The two figures plotted in Fig. 7 show the proven 2+2 doubly eclipsing quadruples. They are compared to the sets of systems published by the group around T. Borkovits (using their photodynamical modelling), our group, and the other systems published by various other authors (unconnected to any previous groups of authors, see e.g. Ádám et al. 2023). As one can see, the systems discovered and analysed by T. Borkovits et al. are usually more compact (i.e. there is a smaller outer-to-inner period ratio), and sometimes have a mutual A-B period shorter than 100× the inner period. Only one of our systems is like that (here presented ASASSN-V J091951.17-593306.9), which is the reason why we had to use both the geometrical LITE as well as the dynamical term in our calculations for this particular system.

Fig. 7. Inner and outer periods for the 2+2 quadruples. Individual proven 2+2 systems: red dots published by T. Borkovits group, black dots by our group, and the blue ones are the rest of other authors. The newly presented systems in this paper are denoted by magenta asterisk signs. The sloped dashed lines denote the period ratios 100, 1000, and 10 000, respectively. One can clearly see that the red dots have in general a shorter A-B period than our systems plotted in black.

In general, we can say that the presented method of ETV detection and discovering multiples with periods of several years provides us with an excellent method of filling the incomplete statistics. When examining the distribution of periods as presented in Tokovinin (2008), one can see that the longer periods are usually discovered by astrometry (interferometry, CPM pairs, etc.), while the shorter ones are discovered by spectroscopy. The gap in between can be effectively filled with our presented ETV method, which preferably finds systems with periods of years to decades.

We also discussed the possibility of interferometric detections of the doubles based on the period of the A-B pair and the distance to the system. However, for most of the systems this is impractical because of their lower brightness and only a small angular separation. On the other hand, for the system CzeV1254, there is other indirect evidence that the system is not a single star, even in the GAIA DR3 catalogue. It provides us with information about the so-called ‘GAIA re-normalised unit weight error’ (RUWE) value, which indicates that this particular star also shows some non-negligible movement (RUWE exceeding 2.5). We can only speculate that this is caused by our detected 15.6-yr orbit. Future observations over the coming years will shed light on the precise orbital parameters for all of the systems, and are therefore desirable.

Data availability

Table A.1 is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/702/A234.

Acknowledgments

An anonymous referee is greatly acknowledged for his/her helpful and critical suggestions. We do thank the ZTF, ASAS-SN, SWASP, Atlas, OGLE, and TESS teams for making all of the observations easily public available. We are also grateful to the ESO team at the La Silla Observatory for their help in maintaining and operating the Danish 1.54m telescope. The research of P.Z., J.K., and J.M. was also supported by the project COOPERATIO – PHYSICS of Charles University in Prague. The research of JM was supported by the Czech Science Foundation (GACR) project no. 24-10608O. We would like to thank the Pierre Auger Collaboration for the use of its facilities. This work is supported by MEYS (Czech Republic) under the projects Czech Republic MEYS LM2023032, LM2023047, and CZ.02.01.01/00/22_008/0004632. The observations by ZH in Veltěže were obtained with a CCD camera kindly borrowed by the Variable Star and Exoplanet Section of the Czech Astronomical Society. This research made use of Lightkurve, a Python package for TESS data analysis (Lightkurve Collaboration 2018). This research has made use of the SIMBAD and VIZIER databases, operated at CDS, Strasbourg, France and of NASA Astrophysics Data System Bibliographic Services. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

References

Appendix A: Online data

All Tables

All Figures

Fig. 1. Analysis of the CzeV1254 system. Left-hand side: Light curves of both disentangled pairs, A and B, based on the TESS data (residuals given in the bottom plots). Right figures: ETV diagrams of both pairs as they revolve around a common barycenter. Full dots stand for primary eclipses, and open circles for the secondary ones. Green symbols represent the TESS data, and red ones our new dedicated observations. The remaining symbols are the other photometric surveys.

In the text

	Fig. 2. Analysis of the ASASSN-V J070838.27-171952.9 system. The description of the plots and colour coding for individual symbols are the same as in Figure 1.
In the text

Fig. 3. Analysis of the ASASSN-V J091951.17-593306.9 system. In the right-hand side figures are the ETV variations plotted twice for pair B. The upper plot is the complete ETV diagram (i.e. apsidal motion + Römer delay + physical delay). The lower figure is plotted after subtraction of the apsidal motion. The blue curve is the final combined fit, the red one stands for the physical delay, and the magenta is the Römer delay (LITE).

In the text

	Fig. 4. Analysis of the WISE J100820.07-731554.0 system. The description of the plots and colour coding for individual symbols are the same as in Figure 1. The ETV diagram of pair B was plotted after subtraction of the slow apsidal motion.
In the text

	Fig. 5. Analysis of the ASASSN-V J143536.01-721459.4 system. The description of the plots and colour coding for individual symbols are the same as in Figure 1.
In the text

	Fig. 6. Analysis of the Corot 310284765 system. The description of the plots and colour coding for individual symbols are the same as in Figure 1. The Corot data are shown as magenta symbols.
In the text

Fig. 7. Inner and outer periods for the 2+2 quadruples. Individual proven 2+2 systems: red dots published by T. Borkovits group, black dots by our group, and the blue ones are the rest of other authors. The newly presented systems in this paper are denoted by magenta asterisk signs. The sloped dashed lines denote the period ratios 100, 1000, and 10 000, respectively. One can clearly see that the red dots have in general a shorter A-B period than our systems plotted in black.

In the text