A systematic search for double eclipsing binaries in Zwicky Transient Facility data

Context. Double eclipsing binaries are gravitationally bound quadruple systems in a 2+2 configuration where both of the binaries are eclipsing. These systems are interesting objects to better understand stellar formation, to investigate the dynamical interaction between the two binary systems or to study certain stages of stellar evolution. Aims. With this work, we aim to determine if double eclipsing binaries can be found using ZTF data and what the difficulties are in doing so. Secondly, we aim to significantly increase the number of known double eclipsing systems and determine how this sample differs from samples of double eclipsing binaries found with other telescopes. Methods. We develop a new method to systematically search for double eclipsing binaries in sparsely sampled light curves. For this we use box-least-squares (BLS) to search for the period of the first binary in the system. We then remove that signal from the light curves, and search the residual light curve again with BLS to find the second period. We applied this method to ZTF light curves of 575 526 eclipsing binaries known in the Gaia eclipsing binary catalogue. Results. We report the discovery of 198 new double eclipsing binary systems. Conclusions. We successfully implemented a method that systematically searches for double eclipsing binary systems in sparsely sampled data. In total 198 new double eclipsing binary systems have been found in 575 526 light curves. The ZTF sample typically contains more short period binaries compared to the TESS sample, but is also able to find systems with longer periods. We expect that at least three to four times more quadruples can be found by applying this method to all ZTF stellar light curves, by increasing the number of data points as a result of longer observations, and by implementing an automatic detection mechanism that replaces visual inspection.


Introduction
Most stars in the Universe can be found in gravitationally bound binaries or higher-order hierarchies (Tokovinin 2014).It is estimated that approximately 4% of all solar-type stars can be found in quadruple systems (Tokovinin 2021), either in a 3+1 or 2+2 configuration.A small fraction of 2+2 quadruple systems are double eclipsing binaries, in which two eclipsing binaries also orbit each other (Zasche et al. 2019).If the motion of the inner pair is not strongly perturbed by the outer companions, the motions of the stars are approximated by stable Keplerian orbits and can survive for a long time (Tokovinin 2014).
Compact double eclipsing binaries are interesting objects of study to better understand stellar formation or to investigate the dynamical interaction between the two binary systems.An important example of such a dynamical interaction is the Kozai-Lidov mechanism where, on long timescales, a periodic exchange between the binaries' eccentricities and inclination can take place (Kozai 1962;Lidov 1962).Furthermore, these objects serve as an astrophysical laboratory when studying certain stages ⋆ Table A.1 is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https: //cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/682/A164 of stellar evolution, such as common-envelope events or Type Ia Supernovae (Kostov et al. 2022).Therefore, the identification of more double eclipsing binaries leads to a larger sample that can be used to study the origin and life cycle of these objects.Quadruples are intrinsically rare and require a favourable alignment of both binaries with the observer in order to be detected.However, in the case of fortunate alignment, the objects can be identified by periodic eclipses in the light curve.Zasche et al. (2019) systematically searched the OGLE-LMC data by visually inspecting the light curves of eclipsing binary systems and found 72 systems of double eclipsing binaries.More recently, light curves of the Transient Exoplanet Survey Satellite (TESS, Ricker et al. 2015) have been used to identify double eclipsing binaries in this way.Zasche et al. (2022b) identified 116 double eclipsing binaries using TESS, and Kostov et al. (2022) detected in total 199 double eclipsing binaries, 18 of which overlapped with Zasche et al. (2022b).
In this paper, we present the search for double eclipsing binaries in Zwicky Transient Facility light curves.In Sect.2, we present the target preselection and the ZTF data.In Sect.3, we discuss the method we used to identify double eclipsing binaries.In Sect.4, we present the analysis of the data.The results are presented in Sect. 5 and we discuss them further in Sect.6.We end the paper with a summary and recommendations for future work in Sect.7.

Data
As part of the Zwicky Transient Facility (ZTF), the Palomar 48-inch (P48) telescope has been imaging the sky every night since 2018 (Graham et al. 2019;Bellm et al. 2019;Dekany et al. 2020).Most of the time ZTF uses the g and r bands, but a small fraction of the observations are also made in the i band.The exposure times are predominantly 30 s for g and r, and 60 or 90 s for i exposures.The median limiting magnitude, averaged over the lunar cycle, is ≈20.5 in all three bands.We used the PSF-photometry-based light curves which are automatically generated for all persistent sources detected in the ZTF reference images (for a full description see Masci et al. 2019) and are publicly available.
In order to avoid searching billions of ZTF light curves, we made a preselection of objects.We processed only light curves of objects that have been identified as eclipsing binaries by Gaia (Mowlavi et al. 2023).Although this preselection inevitably means that we miss any system that has not been identified as an eclipsing star by Gaia first, the advantage is that it strongly reduced the number of light curves we needed to search.
The Gaia eclipsing binary catalogue contains 2 184 477 sources distributed over the whole sky and down to magnitude G ≈ 20.5.While ZTF goes slightly deeper than Gaia, ZTF only covers the sky from the celestial north pole to a declination of −30 deg.This means there are 1 210 001 sources in the Gaia eclipsing binary catalogue that are also in the ZTF footprint and are fainter than G = 12 (the approximate ZTF saturation limit).To download ZTF light curves for these objects, we used ZTFquery1 (Rigault 2018).

Method
Light curves of double eclipsing binaries contain the signal of two eclipsing binaries with different periods.If the light curve is well sampled (e.g. the continuous light curves of TESS), multiple epochs sample each individual eclipse and the individual eclipses are easy to identify.This makes the identification of double eclipsing binaries relatively straightforward (Zasche et al. 2019(Zasche et al. , 2022a;;Kostov et al. 2022).However, if a light curve is sparsely sampled (the typical time between observations is much larger than the duration of an eclipse), individual eclipses are sampled by a few or just a single epoch only and the light curve needs to be phase-folded in order to identify eclipses.Therefore, searching for double eclipsing binaries requires a slightly different approach for sparsely sampled light curves like ZTF.Instead of identifying individual eclipses, we used the fact that both eclipse signals are periodic and used period-finding methods to identify both periods.
We implemented this method using box-least-squares (BLS; Kovács et al. 2002) which is available through the Python package Astrobase (Bhatti et al. 2017).First, the data from the different bands (g and r) were combined by normalising the data in each band to its corresponding mean.Then we performed the BLS period search on the light curve to find the first period, which we defined as period A. In our experience, this period was often half the orbital period and thus the primary and secondary eclipse overlapped.To avoid this overlap, we phase-folded the light curve to twice period A. To remove the signal of period A, we binned the phase-folded light curve using an empirically We also see that not all data points fall nicely on the curve, indicating the presence of another signal.Lower panel: phase folded light curve for binary B, after removal of the signal at period A, with P B = 0.404213 days.determined number of 100 bins.We then divided each flux measurement in the light curve by the mean of its corresponding bin.After the removal of the signal at period A, we applied the BLS method again to find the second period; period B. At each point in this process, we saved the phase-folded and binned light curves for visual inspection.An example of a (phase-folded) light curve can be seen in Fig. 1.

Analysis
To ensure that the method worked correctly, we used a test dataset consisting of 68 double eclipsing binaries found in prior research (Zasche et al. 2022b) for which ZTF data is available.Furthermore, this test sample also helped to establish a method for identifying double eclipsing binaries quantitatively.For this, the relative peak height (RPH) was introduced as RPH ≡ peak periodpower − periodpower σ periodpower (1) A164, page 2 of 10 where max periodpower is the BLS-periodogram peak associated with the best period, periodpower is the mean of powers, and σ periodpower is the standard deviation.In ASTROBASE these periodpowers are referred to as lspvals (Bhatti et al. 2017).
For readers not familiar with the BLS algorithm, an example of a periodogram is shown in the appendix in Fig. A.2.Here we see the periodogram power value of each orbital period.
Additionally, ASTROBASE offers a function to calculate the signal-to-noise ratio (S/N) for each periodic signal (Bhatti et al. 2017).From the test sample we empirically determined that a RPH > 10 and S /N > 10 may indicate the presence of a double eclipsing binary.As it was not feasible to visually inspect all objects in the sample, we selected only candidates with RPH > 10 and S /N > 10.This significantly reduced the number of light curves that needed visual inspection to identify double eclipsing binary candidates.These selection criteria are discussed further in Sect.6.Some objects were, based on a high relative peak height and S/N, initially identified as quadruple candidate but appeared to exhibit a sinusoidal rather than an eclipsing signal.These objects were associated with relatively short periods P B < 0.30 days, which could be an indication of a variable star in the binary system that lies on the instability strip, such as δ Scuti stars (Pietrukowicz et al. 2020).Another common false positive is associated with the imperfect removal of signal A. This would result in period B being the same as period A, or an integer multiple of period A. We therefore removed any candidate for which period B was a multiple of period A.
When testing the method, it appeared that light curves containing high cadence observations (nights with more than 30 epochs) had a negative impact on the ability to find the right periods.This is probably caused by correlated noise introduced by the calibration of the data.Therefore, observations with more than 30 epochs per night were removed.A similar effect was observed by outliers in high and low flux.For this reason, outliers in normalised flux were removed using the interquartile range with a cut-off value of 5. Other unreliable data, flagged by the ZTF data process pipeline, was also removed.After this data pre-processing and cleaning, a minimum threshold of 500 data points per light curve was set for analysis.This decreased the sample size from 1 210 001 to 575 526 objects.
Based on the objects in the test sample, we set the period search range for new objects between 0.1 and 150 days with a stepsize determined by the minimum transit duration and observational period.We chose this lower limit as it is not expected to find any objects with shorter orbital periods (Drake et al. 2014).The upper limit was based on the 5-yr observational period of the ZTF and the decreasing likelihood of finding longer periods.For the transit duration, the BLS performed a search for a periodic signal lasting between 0.01 and 0.2 in phase.
Finally, we applied the pipeline to the 575 526 ZTF light curves.We used the HELIOS computer cluster to process the light curves.The average processing time per object is approximately 10 s.
We used the BLS output statistics to select 5350 promising candidates.These were visually inspected in order to confirm the double eclipsing nature of the objects.At this point we also determined, to the best of our abilities, what the two orbital periods are since the BLS algorithm typically finds half the orbital period.Finally, we also checked if the photocenter position was correlated with the eclipse periods to identify chance alignments of two unrelated eclipsing binaries.The typical RMS positional accuracy is ≈0.1 ′′ .No chance alignments were found in this manner.Here P short is defined as the shortest period of A and B, and P long as the longest period of A and B. The figure shows that most objects have both periods shorter than a few days.Using a kernel density estimate, this distribution seems to be the same for the candidates and test objects.

Results
Of the 575 526 light curves analysed (500 data points or more), 203 (≈0.035%) are identified as double-eclipsing binaries candidates.They are summarised in the appendix in Table A.1, a full version of which is available at CDS.An example is shown in Fig. 1.Almost all candidates for double eclipsing binaries are not currently known in any scientific literature (i.e. they are not listed as such on SIMBAD; Wenger et al. 2000).Only four of these 203 objects have been found in prior research.One of these four objects, Gaia DR3 source ID 2003428628134264448, has previously been identified as a triple star system (TIC 388459317; Borkovits et al. 2021), but the ZTF data suggests that it is a double-eclipsing binary.The other three objects, Gaia DR3 source ID's 2060949991271681664 (Zasche et al. 2022a), 407849617690288128 (TIC 354314226;Kostov et al. 2022), and 158123726424621824 (TIC 150055835;Kostov et al. 2022) have been identified as double eclipsing binary systems by previous studies that used TESS data.Using the ZTF data, we find that the periods of these objects agree between both studies.
The periods of the double eclipsing binaries in the ZTF sample range from 0.11 days to 323.20 days, and the mean of periods A and B are 1.21 days and 10.52 days respectively, see Fig. 2. The figure shows that most objects have both periods shorter than a few days.Using a kernel density estimate, this distribution seems to be the same for the candidates and the test objects.
To understand if there is any bias in our search, all 1 210 001 binaries in the sample were plotted in an HR-diagram, see the left panel of Fig. 3.Here we see that the double eclipsing binaries, plotted separately as dots, are scattered throughout the diagram and do not seem to have any preferred location.
Since it is interesting to not only compare the double eclipsing binaries to the input sample but also to a general star population within the Milky Way, we refer to the HR-diagram of all stars in the Gaia DR3 archive within 200 pc of the sun.This can be seen in the right panel of Fig. 3.The quadruple candidates and test objects are systematically above the single stars in this diagram and are mostly Sun-like stars.We note, however, A164, page 3 of 10 that quadruple systems on the HR-diagram have a lower magnitude of 2.5 log 10 (4) = 1.505 compared to single stars, assuming all stars in the quadruple are similar and equally bright.

Confusion with triple eclipsing systems
As we already noted, Gaia DR3 source ID 2003428628134 264448 was previously identified as a triply eclipsing triple system (Borkovits et al. 2021) (a multistar system where all three stars eclipse each other).We initially identified this system as a quadruple because we detected a primary and secondary eclipse at a period of 88.8 days.However, as is shown by Borkovits et al. (2021) using well-sampled TESS and ground-based follow-up light curves, the eclipse shapes are complex and only consistent with a tertiary star that eclipses, and is being eclipsed by, a short period binary.Triple eclipsing stars are also extremely The fact that one of our double eclipsing binary candidates turned out to be a triply eclipsing triple system, prompted us to reconsider the long period systems that we recovered.We focus on long period systems since, for triple systems, the outer orbit has to be at least 5 times the inner orbit, while a factor of 10 is required for slightly eccentric orbits (Mardling & Aarseth 2001;Borkovits et al. 2021).With just the ZTF data, individual eclipses are not well-resolved, which makes an unambiguous identification of triples difficult.However, an eclipse of a triple eclipsing system is more noisy in a folded light curve because of the complex and changing shape of individual eclipses.We use this to determine whether some systems might be triple-eclipsing systems instead of double-eclipsing binaries.
A164, page 4 of 10 Objects 2020970648976206208, 2210495392378660864, and 2038100043701450112 show somewhat irregular, longduration eclipses in ZTF data and could be eclipsing triple systems.Objects 2058085351143518464, 2003428628134264448, 2201723866572302976, and 413930359370384384 show a primary and a shallow secondary eclipse at period B, with the primary eclipse somewhat irregular, suggesting that they could also be triply eclipsing systems.Object 2041659197183314304 shows well-behaved eclipses, and we judge that this system is unlikely to be a triply eclipsing system.
Objects 207943285475761280 and 2032118077650924672 are two systems for which the eclipse of the long period signal is much deeper in g than r and maybe a shallow secondary eclipse can be seen.This suggests that a large, cold star is eclipsed by a smaller, hotter star.We suspect that these are double eclipsing systems where one of the components has evolved off the mainsequence.For 260945106052343296 and 2012994482376100736 the period is too uncertain to draw any conclusions.
Since we cannot definitively classify these objects as triple eclipsing systems, we consider all these objects double eclipsing binary candidates in the rest of this paper.Follow-up light curve observations of the eclipses are needed for each object to definitively determine their nature.

Orbital period distribution
In general, binaries with shorter periods seem easier to detect than binaries with longer periods.One explanation is that stars in binaries with longer orbital periods are further apart from each other.As the orbital distance increases, the eclipse probability decreases since this requires a more precise alignment with the observer for the eclipse to be detected.In addition, longer orbital periods thus lead to shorter eclipse durations as fraction of the orbit.Some light curves showed that the eclipse duration of binary B was either very short (<0.01 in phase), or contained very few data points, which makes it difficult for the algorithm to identify the eclipse.This agrees with the results shown in Fig. 2, where most double eclipsing binary candidates have both periods shorter than a few days.The orbital periods of the candidates found in this research seem fairly similar to those found by Kostov et al. (2022), Zasche et al. (2022b), see Fig. 4.However, one difference is that we find some shorter periods in this research.Whereas previous research shows that most objects show a period ratio of P B /P A between 1 and 2 (Kostov et al. 2022), in Fig. 5 we see that here also much smaller ratios are found, indicating a larger difference between period A and period B. An explanation for this can be the data cadence of TESS, which is 26 days, meaning it most likely does not find any significant periods larger than approximately 13 days.Furthermore, short 30 min exposures by TESS may lead to the spreading out of short eclipses, which is why TESS might not find very short periods either.We also note that all stars in our sample are listed in the Gaia catalogue.Any bias of Gaia towards shorter periods may therefore reflect in our results.

Detection efficiency and completeness
In this section, we briefly discuss the detection efficiency and completeness of our search.Of the 575 526 analysed objects, at least 203 (≈0.035%) passed the visual inspection tests and  2022b), and the ZTF sample from this work.Here P short is defined as the shortest period of A and B, and P long as the longest period of A and B. Fig. 5. Period ratios P short /P long of each object.Here P short is defined as the shortest period of A and B, and P long as the longest period of A and B. The vertical lines indicate where potential harmonics would fall.We note that integer multiple harmonics were removed in the object search.However, noninteger harmonics were not removed.Here we see that most objects do not fall on any of the harmonics lines indicating that there does not seem to be any particular correlation between the periods.
showed to be a double eclipsing binary candidate rather than a binary system.This is approximately a factor 4.5 lower success rate than (Zasche et al. 2022b; 116 out of 70 000) and a factor 7 lower than (Pawlak et al. 2013 15 out of 6138).However, we note that OGLE-III has observed all stars in a dense area in the Small Magellanic Cloud for a longer period of time (Pawlak et al. 2013).Zasche et al. (2022b) have, like this research, only looked at eclipsing binaries, but a difference with this research is the data cadence, which for ZTF typically is one observation per A164, page 5 of 10 Fig. 6.Detection efficiency of double eclipsing binaries by number of data points and apparent magnitude.Left: cumulative histogram for the number of data points per light curve.After passing a threshold of approximately 1000 data points there is a significant increase in the quadruple detection efficiency.Note that the plateau and discontinuity at around 2000 data points is a result of the ZTF survey strategy.Right: cumulative histogram for the magnitude of stars for the complete sample and quadruple candidates in the magnitude range 13-20.5.We see that for stars fainter than G ≳ 18 the quadruple detection efficiency rapidly approaches zero.The sharp increase at G ≈ 15 is not fully understood.night only for 5 yr.TESS, on the other hand, provides continuous undisturbed photometry for 26 days (Zasche et al. 2022b).Lastly, the angular resolution of ZTF, OGLE, and TESS are very different.This, combined with the new method used here, makes it difficult to compare success rates, but could explain the difference.
In the rest of this section, we briefly discuss a few causes that limit the completeness of our search.First, we consider the efficiency of our quadruple detection method.Our method was able to retrieve the correct periods for most known objects in the test dataset.For periods that could not be recovered, further investigation revealed why this was not possible.In some cases, many data points were flagged as unreliable, leaving only very few data points in the light curve for period finding.In other instances, the eclipse depth of one binary was very small, in the order of a few percent, making it difficult for the algorithm to detect any periodic signal.Other reasons why the algorithm was not successful included that some light curves showed a very small trend in flux over time, resulting in the cadence of the telescope being found as a significant period; this typically being one or two days.Even light from the Moon coming in as background noise could lead to the period of the Moon being found as a significant period.For compact double eclipsing binaries, eclipse-timing variations (ETV) can smear periodic signals which may affect the period finding.However, when comparing our results with known systems with ETV, we found this not to be the case.We therefore expect this method to be suitable for finding compact double eclipsing binaries.
The completeness of our search for double eclipsing binaries is also limited by our use of the Gaia eclipsing binary catalogue.When we look at the results of other studies (Kostov et al. 2022;Zasche et al. 2022b), we see that approximately 30% of the found quadruples are in the Gaia catalogue.In this research, we have only looked at ZTF objects that are known in the Gaia catalogue.This means that potentially three times more quadruples can be found in the ZTF data that are not in the Gaia catalogue.
As it was not practical to look at 575 526 light curves, we used automatically calculated statistics to significantly decrease the number of light curves to inspect, which could also affect the completeness of our search of the ZTF data.We used the test objects to empirically determine cuts; relative peak height >10 and S /N > 10.As can be seen in Fig . A.3, the wide range of values of relative peak height and S/N cannot reliably indicate whether an object might be a quadruple candidate.This suggests that more double eclipsing binaries can be found for values of relative peak height and S/N below 10.It can also be seen in the figure that some objects did have a high relative peak height and S/N but were not identified as a quadruple candidate.This could have several reasons.In some cases, the algorithm was not able to completely remove the first periodic signal, which means that the signal, or an harmonic, was returned again for period B. Other explanations include, but are not limited to, an outlier in the light curve, the cadence of the telescope, or the influence of the Moon.Sometimes the algorithm did not find an eclipsing signal but a sinusoidal signal, possibly due to the pulsation of one of the stars in the eclipsing binary.Even though this is not a double eclipsing binary, a strong periodic signal was still found, resulting in a high relative peak height and S/N.
We suspected that the detection efficiency of double eclipsing binary stars is also strongly affected by the number of epochs in the light curve.If we look at the left panel in Fig. 6, we see a cumulative histogram for the number of data points per light curve.It shows that, only after passing a threshold of approximately 1000 data points, the quadruple detection efficiency significantly increases.After that, there are more double eclipsing binaries as the number of data points increases.The difference, however, is small.This means that, as ZTF keeps collecting data (or ZTF is combined with external data), we can expect the number of detectable double eclipsing binaries to double.
We also briefly consider to what brightness we can reliably find double eclipsing binaries in ZTF data.On the right panel in Fig. 6 we see a cumulative histogram for the magnitude of stars, both for the complete sample and for the quadruple candidates.The ZTF has looked at stars in the magnitude range 13-20.5.Since the curve of the quadruples is constantly above the overall sample curve, this indicates that it is easier to find quadruples for bright stars compared to faint stars.Almost all quadruples have a magnitude of G ≲ 18, while only half of the input sample is fainter than this limit.We conclude that the light curve signalto-noise ratio is too low for sources fainter than G ≈ 18.As can be seen in Masci et al. (2019), the RMS precision is ≈1% for A164, page 6 of 10 sources brighter than 18 magnitude, but the precision degrades rapidly for fainter sources.
To summarise, we expect that future research can find at least three to four times more quadruples by including the following considerations.First, the search should be not limited to the Gaia eclipsing binary catalogue but the entire ZTF database can be searched instead.Secondly, systems with a relative peak height and S/N lower than 10 should also be considered.Thirdly, longer observation time of the ZTF will lead to more data points which has a positive effect on quadruple finding.Lastly, by using an automatic detection method that can replace visual detection, for instance machine learning, all the light curves can be analysed more efficiently.

Conclusions and future work
We have developed a method to systematically search for double eclipsing binaries in sparsely sampled data.Using this method we found 198 new double eclipsing binaries candidates in the ZTF data.The periods of the objects ranged from 0.11 days to 323 days, with a mean of periods A and B of 1.21 days and 10.52 days respectively.By searching the entire ZTF database, increasing the number of data points per light curve, and implementing an automatic detection mechanism that replaces visual inspection, we expect that at least three to four times more quadruples can be found using this method.Other recommendations for future work are; to apply this method to data collected by other telescopes such as Gaia.In addition, we recommend that the multi star systems are included in the target lists of large multiplex spectroscopic surveys such as SDSS, DESI, and WEAVE, in order to obtain phase-resolved spectra of the quadruple candidates to further investigate their properties.Lastly, detection of eclipse-timing variations in the double eclipsing binaries can definitively prove their quadruple nature.

Fig. 1 .
Fig. 1.ZTF light curves of Gaia DR3 source ID 2216420454386370176.Upper panel: normalised flux as a function of time.It shows two sets of eclipses with P A = 1.072780 days and P B = 0.404213 days.Middle panel: phase folded light curve for binary A with P A = 1.072780 days.We also see that not all data points fall nicely on the curve, indicating the presence of another signal.Lower panel: phase folded light curve for binary B, after removal of the signal at period A, with P B = 0.404213 days.

Fig. 2 .
Fig. 2. Periods A and B of 13 test objects and 203 quadruple candidates.Here P short is defined as the shortest period of A and B, and P long as the longest period of A and B. The figure shows that most objects have both periods shorter than a few days.Using a kernel density estimate, this distribution seems to be the same for the candidates and test objects.

Fig. 3 .
Fig. 3. Double eclipsing binaries in HR-diagrams.Left: HR-diagram of all 1 210 001 binary objects in the sample.The double eclipsing binary candidates and test objects, separately plotted as dots, are scattered throughout the diagram and do not seem to have any preferred location.Using the sun as a reference with M G = 5 and G BP − G RP = 1, a larger number of objects seem to be similar to the sun.Right: an HR-diagram of all stars in the Gaia DR3 archive within 200 pc of the Sun.The candidates and test objects are systematically above the single stars in this diagram and are mostly sun-like stars.The arrow with a length of 1.505 magnitude indicates where one star in the quadruple system would lie, assuming all stars in the system are equally bright.

Fig. 4 .
Fig. 4. Orbital periods of double eclipsing binaries of the three largest samples: the sample from OGLE (Zasche et al. 2019), the sample from TESS by Kostov et al. (2022) and Zasche et al. (2022b), and the ZTF sample from this work.Here P short is defined as the shortest period of A and B, and P long as the longest period of A and B.

Fig
Fig. A.1: Right ascension (RA) and declination (DEC) of the double eclipsing binary candidates and test objects in this research.As expected, most candidates lie in the Galactic plane.No objects are seen below a declination of -30 • as this is not visible by the ZTF.

Fig
Fig. A.2: Periodogram of Gaia DR3 source ID 2216420454386370176 showing the lspvals as a function of orbital period.Left: The periodogram power values for period A. Right: The periodogram power values for period B after removing the signal of period A.