© The Authors 2025

1. Introduction

Binary stars have played a key role in stellar astrophysics as they provide a means to accurately measure the physical parameters of stars, including their masses. In the case of eclipsing binaries (EBs), the light curves (LCs) embed information about the component stars’ radii, temperatures, and luminosities. Among the least understood characteristics of EB LCs is the fact that the out-of-eclipse maxima are different. This asymmetry is known as the O’Connell effect (O’Connell 1951). Although it affects a relatively small fraction of EBs (of order 7%, according to Knote et al. 2022), it complicates the measurement of their physical parameters, leading to errors in estimated stellar radii, luminosities, and orbital inclinations. Furthermore, the effect may mimic the presence of a third body, further complicating interpretation.

O’Connell (1951) found, based primarily on visual and photographic data, predominantly brighter first maxima. Davidge & Milone (1984) revisited and expanded his work using UBV photoelectric photometry, arriving at somewhat different conclusions. In particular, they observed a more even distribution of asymmetries and found that the brighter maximum was often redder. Their study suggested that multiple mechanisms, including magnetic activity, circumstellar material, and tidal distortions, can contribute to the effect, thus highlighting the need for further investigation.

Several interpretations have been proposed to explain the O’Connell effect. They include the presence of cool stellar spots due to the magnetic activity of the component(s) (Binnendijk 1960), hot spots as the result of gas stream impact due to mass transfer (Zhai & Fang 1995; Hilditch et al. 1997; Pribulla et al. 2011), and clouds of circumstellar material (dust and gas; Liu & Yang 2003). In particular, cool spots have been demonstrated as a viable explanation by means of Doppler imaging techniques (e.g., Hendry & Mochnacki 2000; Şenavcı et al. 2011; Niu et al. 2022). The other two scenarios have gained popularity recently. The magnetic activity interpretation applies more to late-type systems due to the presence of deep convective zones and rapid rotation. For earlier-type contact systems (A–F spectral types), on the other hand, the LC asymmetries are more probably due to mass transfer from the more massive star: as the primary star transfers material through the first Lagrangian point to its companion, the impact generates a region of enhanced brightness on the receiving star. When this hot spot rotates into view during the orbital cycle, it causes one of the LC maxima to appear brighter than the other, thus producing the O’Connell effect (Pribulla et al. 2011).

Several models have been used for the LC phenomenological modeling of EBs (e.g., Prša et al. 2008; Tkachenko et al. 2015; Mikulášek 2015) and/or to search for contact systems presenting the O’Connell effect (Wilsey & Beaky 2009; Johnston et al. 2019; Papageorgiou et al. 2014, 2017; Knote et al. 2022). One of the primary quantitative indicators of the presence of the O’Connell effect is Δm_max, defined as the difference in magnitude between the two maxima of the LC. Using HIPPARCOS contact binaries presenting the O’Connell effect, Pribulla et al. (2011) showed that, for late-type systems, the signs of Δm_max at quadrature phases are random, while for the earlier systems there is a clear tendency for the first quadrature maximum to be brighter than the second.

Although Δm_max can give us some information about how strong the effect is in different systems, other useful indicators have also been defined. In particular, McCartney (1999) introduced the O’Connell effect ratio (OER) and light curve asymmetry (LCA) metrics, both of which provide further insight into the asymmetries of the LCs. OER is defined as the ratio of the area under the curve for phases ϕ = 0 − 0.5 and ϕ = 0.5 − 1. This means that a LC with an OER = 1 has the same amount of flux in both halves, while an OER ≠ 1 means that one of the two halves has more flux than the other, indicating an asymmetry. LCA quantifies deviations from symmetry in EB LCs by directly assessing the point-wise asymmetry between corresponding phase bins. Unlike the OER, which measures total flux differences, LCA provides a more localized evaluation of asymmetry. In like vein, LCA values greater than 0 indicate that the LC is asymmetric, whereas LCA = 0 means that the LC is highly symmetric.

This work presents a multiband LC analysis of the O’Connell effect in 14 EBs using photometric data from the Zwicky Transient Facility (ZTF; Bellm et al. 2019; Masci et al. 2019) and the Asteroid Terrestrial-impact Last Alert System (ATLAS; Tonry et al. 2018a). We also present a new way to model LCs with strong asymmetries, which is different from the commonly used Fourier series fit, and provide the OER, LCA, and Δm_max parameters for each of these systems in four different passbands.

In Sect. 2 we present the data used in our study. In Sect. 3 we describe our procedure for identifying O’Connell EBs therefrom. Our new phenomenological model for these systems is described in Sect. 4, and the results of its application to our multiband LCs are presented in Sect. 5. Finally, a discussion of our results and conclusions can be found in Sect. 6.

2. Observations

The data used in this project came from two main sources, namely the ZTF and ATLAS surveys.

2.1. ZTF

ZTF is a Samuel Oschin 48-inch Schmidt telescope located at Palomar Observatory with a 47 deg² field of view (FOV; Bellm et al. 2019; Masci et al. 2019). This large FOV ensures that ZTF can scan the entire northern sky every night. It can observe 3760 deg² per hour, reaching a 5σ detection limit of 20.5 mag in r with a 30 s exposure. Its primary science aims are the physics of transient objects, stellar variability, and Solar System science (Graham et al. 2019). ZTF photometry for public surveys is obtained in two observation modes in the g, r, and i bands, namely three-night cadence for the entire visible northern sky from Palomar and one-night cadence for the Galactic plane (Galactic latitude |b|≤7°). In our analysis, the dataset that was used corresponds to ZTF Data Release 15 (DR15). Perfectly clean extractions (catflag = 0, |sharp| < 0.25, and mag < limitmag − 1.0) were used in both g and r bandpasses at every epoch. Since ZTF’s DR15 did not include i-band data for many of our sources, we did not include this band in our analysis. DR15 contains on average ∼ 637 epochs of observations in g and ∼850 in r per star in our sample, with observations covering the time span between March 2018 and November 2022.

2.2. ATLAS

ATLAS is a survey designed to detect small asteroids (10–140 m) on their “final plunge” toward impact with Earth (Tonry et al. 2018a). It observes about 13 000 deg² of sky at least four times per night, which also makes the resulting dataset very useful for the discovery of transient events and the follow-up of variable stars. ATLAS uses Wright Schmidt telescopes with a 4.5 deg² FOV each, and each ATLAS telescope takes four 30 s exposures of 200–250 target fields per night. Although ATLAS is tuned for investigating near-Earth asteroids, it collects the LCs of transients and variable stars down to a magnitude limit of r ∼ 18 in the c (∼ g + r) and o (∼ r + i) bands, for which approximate color transformations linking them to the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) g, r, and i bands (Magnier et al. 2020) were provided by Tonry et al. (2018a).

ATLAS observations of all our objects were taken between July 2015 and October 2020. The ATLAS photometry we incorporated in our analysis was extracted directly from the ASCII files in the ATLAS main storage computer. Finding our objects of interest required going through the ATLAS structured query language (SQL) database, which organizes the exposures based on multiple keywords. A coordinate-based search for the stars previously selected in the ZTF survey was performed in the database to locate the exposures that include the objects of interest. Knowledge of the exposures allowed us to locate the final photometry catalogs of those images (the “.dph” files). We developed scripts to scan these catalogs for detections matching the object coordinates within one ATLAS pixel, and the final photometry of the matching objects was extracted together with the exposure dates.

ZTF J020930.46+513905.4 ZTF J014355.46+354402.7 ZTF J054353.09+014433.7 ZTF J093054.62+745619.7 ZTF J091456.62+394615.5 ZTF J011757.02+305050.8 ZTF J164516.11+321905.2 ZTF J005348.21+465518.6 ZTF J181310.67+761840.2 ZTF J155632.27+654934.9 ZTF J201843.21+540656.4 ZTF J200516.81+520253.7 ZTF J014744.51+450636.8 ZTF J182814.00+551351.2

The photometry in ATLAS .dph files was done with a version of DoPhot (Schechter et al. 1993) that was refurbished by Alonso-García et al. (2012), with further improvements by John Tonry (private communication). The .dph files are the final stage of the ATLAS reduction pipeline for non-moving objects; it builds upon previous stages that provide astrometric and photometric solutions based on several different catalogs (Tonry et al. 2018b) and photometric zero points that vary across the large charge-coupled device (CCD) field and include correction for nonuniform cloud cover. More information on ATLAS photometry is given by Heinze et al. (2018).

3. Candidate selection

Our initial search sample consisted of 640 objects (primarily on the Galactic plane) with ZTF g- and r-band LCs, classified as detached or semidetached EBs by the ALeRCE broker (Förster et al. 2021; Sánchez-Sáez et al. 2021). These objects, retrieved via direct database queries between April and June 2020, were required to have a minimum of 150 combined g + r ZTF detections, as reported by ALeRCE, to ensure well-sampled LCs. All the objects were reclassified as contact binaries in the ZTF Catalog of Periodic Variable Stars (Chen et al. 2020). The O’Connell effect is strongly associated with close binary interactions, particularly in contact systems where mechanisms like mass transfer, gas streams, and tidal distortions dominate. Accordingly, this sample provided a suitable initial basis for identifying candidates exhibiting LC asymmetries indicative of the O’Connell effect.

We then performed an automated search for the O’Connell effect on all 640 candidates using the code provided by Papageorgiou et al. (2014), Papageorgiou et al. (2017), which computes Δm_max values for each star. At this stage, the periods reported by Chen et al. (2020) were used. Candidates were chosen that satisfied the criterion Δm_max > 0.05 mag in ZTF’s r band; visual inspection of the LCs revealed that this cutoff is sufficient to reliably identify the presence of the O’Connell effect in the ZTF light curves. We visually inspected all the remaining ZTF DR15 LCs, ending up with a sample of 16 EBs that show clear evidence of the O’Connell effect. These stars are listed in Table 1. Temperatures were found by cross-matching this sample with Transiting Exoplanet Survey Satellite (TESS) Input Catalog Version 8 (Stassun et al. 2019). Out of these 16 stars, 14 with well-defined photometry available in ATLAS were retained for further analysis.

Finally, in order to refine the periods of the sources, we used the photometry in ZTF and ATLAS to compute a multiband Lomb-Scargle periodogram (VanderPlas & Ivezić 2015). This was accomplished using the gatspy package (Vanderplas et al. 2016) in python. The periods listed in Table 1 are the result of this analysis. In all cases, they are very similar to the period values previously reported by Chen et al. (2020).

Fig. 1. Phenomenological modeling of the O’Connell EB candidate ZTF18aabexej. Top panels: Four components of our model, as described in the text (Sect. 4). Bottom panel: final model fit to ZTF g-band data.

4. New phenomenological model

First, we phase-folded the time axis of the LCs of each object using the usual equation:

$$\begin{array}{c}\hfill \varphi =\frac{T-T_0}{P}-\text{Int}\left(\frac{T-T_0}{P}\right),\end{array}$$(1)

where ϕ is the phase, T is the time, T₀ is an arbitrary reference time, and P is the period. We set T₀ so that it corresponds to the time of a primary minimum, thus ensuring that phase zero corresponds to the deeper eclipse in the phase-folded LCs. After that, we converted the magnitudes of each filter into normalized fluxes (Warner 2016) using the following equation:

$$\begin{array}{c}\hfill I[m(\varphi )]={10}^{-0.4\times [m(\varphi )-m(\text{max})]},\end{array}$$(2)

where m(ϕ) is the magnitude at phase ϕ, and m(max) corresponds to the maximum value of the magnitude in each filter. We then performed a model based on a so-called super Gaussian (sG; Murakami et al. 2022), formulated as follows:

$$\begin{array}{c}\hfill I{(\varphi )}_{\text{sG}}=B_0-B_1\text{exp}\{-{\left[\frac{{(\varphi -\mu )}^2}{2{\sigma }^2}\right]}^p\},\end{array}$$(3)

where B₀ is the mean of the non-varying ground base, B₁ is the height, μ is the mean (location), σ is the width, and p is the squareness of the sG. Our model is a combination of this sG and a truncated third-order Fourier fit, given by

$$\begin{array}{c}\hfill I{(\varphi )}_{\text{Fourier}}=a_0+{\displaystyle \sum _{n=1}^3}[a_{n}cos(2\pi n\varphi )+b_{n}sin(2\pi n\varphi )],\end{array}$$(4)

where a₀ corresponds to the mean, and (a_n, b_n) are the n-th Fourier coefficients. The modeling process is composed of five stages:

1. First, we calculate m(max) by performing a simple second-order polynomial fit around the brighter maximum. Magnitudes are then converted into normalized flux using Eq. (2), whereupon we remove bad-quality data by filtering out all points with errors >  0.08 in normalized flux.

2. Following that, we extend the LC for each filter from phases ϕ = −0.2 to ϕ = 1.2. This is done to facilitate the modeling of the primary eclipses, which take place at phase 0.

3. We then model the two eclipses independently using the sG from Eq. (3). The primary eclipse is modeled over the range ϕ = 0.8 to ϕ = 1.2 in phase (hereafter the sG-1 component), whereas for the modeling of the secondary we use phases from ϕ = 0.2 to ϕ = 0.8 instead (the sG-2 component).

4. Next, we replicate the same sG primary eclipse model at ϕ = 1.0 (i.e., the sG-1 component above) to match the primary eclipse centered on phase ϕ = 0 (the sG-3 component)¹.

5. Finally, when we find the best model for each section of the LC, all of them are combined, and the third-order Fourier fit given by Eq. (4) is added:

   $$\begin{array}{c}\hfill I{(\varphi )}_{\text{fit}}=I{(\varphi )}_{\text{sG}\ast }+I{(\varphi )}_{\text{Fourier}}.\end{array}$$(5)

   Here I(ϕ)_sG* corresponds to the sum of the three sG models described above, and I(ϕ)_Fourier is the truncated third-order Fourier fit, which we use to ensure that the sum of the three sG models properly matches the combined data. Initial parameters for this Fourier fit are estimated using the symfit python library (Roelfs & Kroon 2021). Although modeling is done from phases ϕ = −0.2 to ϕ = 1.2, we use the LC fit from phases 0 to 1 for further analysis. An example of the four components and the final model is presentedin Fig. 1.

Appendix A describes the difference between our model and the most commonly used 12-term Fourier fit (Wilsey & Beaky 2009). The results of the LC fit to our 14 O’Connell EB candidates, which was carried out using lmfit (Newville et al. 2014, 2016) in python, are presented in Appendix C.

Using the model thus obtained, we next calculate ΔI_max using Eq. (2). This is done using two different methods. In the first method, we calculate ΔI_max from the model values at quadrature phases. This value is denoted in Table 2 as ΔI_sG. Uncertainties in this quantity were computed using a bootstrap approach and assuming that the errors of the measurements were normally distributed. In the second method, we again perform a simple second-order polynomial fit around each of the maxima derived from our model. This gives us another value for the difference between maxima, which is denoted ΔI_poly in Table 2. The difference between ΔI_poly and ΔI_sG is further explored inAppendix B.

Following the definition provided by McCartney (1999), we also calculate the OER and the LCA as

$$\begin{array}{c}\hfill \text{OER}=\frac{{\displaystyle {\int }_{0.0}^{0.5}}[I{(\varphi )}_{\text{fit}}-I_0]d\varphi }{{\displaystyle {\int }_{0.5}^{1.0}}[I{(\varphi )}_{\text{fit}}-I_0]d\varphi }\end{array}$$(6)

and

$$\begin{array}{c}\hfill \text{LCA}={\left\{{\displaystyle {\int }_{0.0}^{0.5}}\frac{{[I{(\varphi )}_{\text{fit}}-I{(1.0-\varphi )}_{\text{fit}}]}^2}{I{(\varphi )}_{\text{fit}}^2}d\varphi \right\}}^{1/2},\end{array}$$(7)

respectively, where I(ϕ)_fit is our model given by Eq. (5) and I₀ is our baseline, defined as the minimum value of our model (at ϕ = 0).

A user-friendly version of the code used for the phenomenological modeling described in this section is available on GitHub²

5. Results

The results for all 14 objects, including our measured ΔI, OER, and LCA values, are provided in Table 2. The uncertainties in the derived OER and LCA values were calculated using Monte Carlo integration to solve Eqs. (6) and (7) with a 95% confidence interval.

As expected, a strong anticorrelation between ΔI and OER can be seen in Fig. 2. In Fig. 3, where OER values are shown as a function of LCA, a bifurcated distribution of OER values around OER = 1 can be seen, confirming previous findings (McCartney 1999). Stars in our sample are, by and large, split between a subset whose OER values increase with LCA and another for which the OER decreases with increasing LCA, with few exceptions; in both cases, this indicates increasing asymmetry for larger LCA values. Asymmetry tends to be greater in the bluer passbands, again as expected, as cooler spots on the stellar surface have a stronger contrast against their surroundings in the blue due to the increased sensitivity of shorter wavelengths to temperature variations.

Fig. 2. Relation between the OER and ΔIsG for the 14 EBs in our sample. Different symbols and colors denote different filters, namely r, g, c, and o, following the key at the top right of the plot. The dashed line represents a least-squares fit to the data for all filters combined. The coefficient of determination is r2 = 0.831, and the resulting fit is given by OER = −2.659 ΔIsG + 1.0.

Fig. 3. Same as Fig. 2 but for the relation between the OER and LCA. The dashed line represents the OER = 1 locus.

This is explored further in Fig. 4, where we examine how the difference between the maxima varies as a function of the effective wavelength (λ_eff) for each star in our sample. While for most objects the absolute ΔI_sG value indeed decreases toward the redder filters, there are exceptions – which is perhaps not unexpected, since the filters we used are broad and have significant overlap in wavelength coverage (compare the ZTF and ATLAS transmission curves in Bellm et al. 2019 and Tonry et al. 2018a, respectively). Figure 5 is similar to Fig. 4 but shows instead mean ΔI_sG values for each band over the entire sample. There is again a clear trend, with the exceptions of filters g and c, whose behavior may also be partly explained by their significant overlap in wavelength coverage.

Fig. 4. Absolute values of the difference between the maxima for all the objects in each of the four filters used in this study, plotted as a function of their effective wavelengths. Symbols have the same meaning as in Fig. 2.

Fig. 5. Absolute mean values of the difference between the maxima for all the objects in each of the four filters used in this study, plotted as a function of their effective wavelengths. Symbols have the same meaning as in Fig. 2.

To explore correlations between ΔI_sG and other stellar properties, we used Spearman’s correlation coefficient, ρ (Spearman 1904) and its associated p-value. This p-value helps indicate if the values of ρ are consistent with the null hypothesis that the characteristics are uncorrelated with ΔI_sG.

Unlike Knote et al. (2022), we find that ΔI_sG is significantly correlated with both T_eff and the period. The period–color and period–temperature relations in contact EBs that underpin this result have been known for a long time (e.g., Eggen 1967; Gettel et al. 2006; Qian et al. 2017; Jayasinghe et al. 2020). In Fig. 6 we present the trends between ΔI_sG and temperature (upper panel) and |ΔI_sG| and temperature (bottom panel) alongside Knote’s sample. We find that the anticorrelation between T_eff and ΔI_sG in each filter is stronger than the corresponding anticorrelation involving |ΔI_sG| in the same filter. In particular, we find that a positive O’Connell effect, in which the brighter maximum occurs before the secondary minimum, mainly occurs in relatively cool systems (T_eff ≲ 5000 K), while a negative O’Connell effect occurs mainly in systems with T_eff ≲ 6000 K. This is supported by Knote et al. (2022), who find a similar temperature threshold for negative O’Connell systems. Both of these results agree with the idea that starspots might be the main cause of the O’Connell effect, as their impact is expected to be more pronounced in cooler stars (Strassmeier 2009).

Fig. 6. Effective temperatures of our sample stars plotted as a function of ΔIsG (upper panel) and |ΔIsG| (bottom panel). Symbols have the same meaning as in Fig. 2. Gray corresponds to data from Knote et al. (2022).

Although the anticorrelation between ΔI_sG and temperature is similar for different filters, it becomes slowly less pronounced toward redder effective wavelengths. When comparing positive and negative O’Connell effects (as revealed by the ΔI_sG values), we observe differences in correlations between ΔI_sG and the aforementioned parameters (see Table 3); despite the lower statistical significance due to the reduced sample sizes, in almost all filters the correlation has different signs, which depend on the sign of the O’Connell effect. The exception to this is the o filter, which shows the same anticorrelation for both positive and negative ΔI_sG. Another particular case is c, which shows no correlation at all for negative ΔI_sG. While a much larger sample is needed to put these results on a firm footing, this preliminary analysis supports the idea that the O’Connell effect is caused by different mechanisms, depending on the ΔI_sG sign. Unlike Knote et al. (2022), we find that |ΔI_sG| does not present significant correlations, and also that ΔI_sG has a strong anticorrelation with the temperature. Although statistical uncertainties are much larger in our case, our sample seems to behave differently than theirs. This difference might be explained by the fact that their sample lacks systems with larger ΔI_sG values. Figure 6 shows how the majority of our sample extends beyond both extremes of the data in their study, and how the correlation is much weaker, especially for positive values of ΔI_sG. These results are also consistent when we limit Knote’s sample to between 4500 and 6000 K to match temperatures similar to those of our sample. When this limit is imposed on Knote’s sample, we find similar correlation coefficients: (ρ = 0.269, p = 0.0003) and (ρ = −0.135, p = 0.075) for ΔI_sG and |ΔI_sG|, respectively. This might suggest that systems with higher ΔI_sG values (negative or positive) do not behave the same way as systems with a weaker O’Connell effect. To further investigate this possibility, a large sample of systems with high ΔI_sG values and a much wider range of temperatures than in our study is needed.

The other correlation worth noticing is the one between ΔI_sG or |ΔI_sG| and the period (Fig. 7, upper and lower panels, respectively). Similarly to correlations involving the temperature, we find a stronger anticorrelation with ΔI_sG than |ΔI_sG| across all passbands. The correlations are shown in Table 4, separately for positive and negative |ΔI_sG| systems. We find that generally this anticorrelation is stronger when the O’Connell effect has a negative sign. Although this is not in agreement with the results obtained by Knote et al. (2022) – who find that a correlation between ΔI_sG and the period (the “O’Connell effect size” in their terminology) is slightly positive in all cases, but stronger in negative O’Connell effect systems – they also find that these correlations are noticeable for systems with P ≤ 0.5 d. This is also true for our sample with short periods, P ≤ 0.3 d.

Fig. 7. Correlation between ΔIsG and the period for our sample. Symbols have the same meaning as in Fig. 2. Gray corresponds to data from Knote et al. (2022), limited to the range 0.2–1 days for clarity.

The presence of correlations for ΔI_sG but not |ΔI_sG| indicates that the physical mechanisms driving positive and negative O’Connell effects may differ fundamentally. While starspots have been proposed as the dominant explanation, with larger ΔI_sG values corresponding to greater spot coverage (Kouzuma 2019; Knote et al. 2022), our results reveal some inconsistencies with this interpretation.

The observed anticorrelations with both period and temperature are consistent with the known period–color relation in contact binaries, where shorter-period systems tend to be cooler. However, if spots were the primary cause, we would expect |ΔI_sG| to correlate strongly with both parameters, as demonstrated by Kouzuma (2019) for W-type (W Ursae Majoris, or EW) binaries (see Sect. 6). Instead, we find these correlations to be weak or absent in certain passbands.

The band dependence of these trends, including the complete lack of correlation for negative ΔI_sG in the c filter, further complicates the spot interpretation. While our results remain qualitatively consistent with Kouzuma (2019)’s findings regarding spot properties, the quantitative discrepancies suggest that additional mechanisms may contribute to the O’Connell effect, particularly for systems exhibiting large asymmetries. This conclusion aligns with recent suggestions that close binary interactions may play a significant role in generating the observed LC asymmetries (Guo et al. 2025; Wang et al. 2025).

6. Discussion and conclusions

As shown in Table 2, the percentage difference between maxima is in the range [1–13%], while the ratio of positive-to-negative signs of ΔI_sG is close to 0.5. Cross-matching our sample with TESS Input Catalog Version 8 (Stassun et al. 2019), we found that all of the 14 EBs in our sample are late-type contact binaries with spectral types [K4V-G7V] that have effective temperatures (T_eff) cooler than 5600 K, as shown in Table 1. Combined with the orbital periods (≤ 0.3 d), this suggests that the contact EBs are rather W subtypes. Since both components are fast-rotating cool stars, their surfaces should show deep convective activity. Thus, the most probable mechanism for the LC asymmetries in our sample is related to the presence of magnetic fields. In addition, due to high chromospheric activity, one might expect these systems to also show strong X-ray emissions. Nevertheless, none of these EBs are found in the fourth XMM (X-ray Multi-Mirror Mission)-Newton serendipitous source catalogue (Webb et al. 2020) within a 10″ search radius. The absence of X-ray detections is most likely attributable to the faintness of these particular systems. This is unlike brighter late-type contact binaries, which are frequently identified as X-ray emitters and for which it is often feasible to model their X-ray spectra (e.g., Brickhouse et al. 2001). Despite the X-ray non-detections, data from the Galaxy Evolution Explorer (GALEX) Data Release 6 (Bianchi et al. 2014) indicate that 6 out of these 14 EBs (within a 3″ radius) are likely associated with near-UV emission, which can also be an indicator of chromospheric activity.

Spectroscopic observations around the activity-sensitive Ca II H&K and Balmer (Hα) spectral lines (Wilson 1966; Long et al. 2021) could show equivalent width variations and/or emission due to chromospheric activity. Long-term photometric monitoring could also reveal out-of-eclipse light variations as the result of starspot migration and/or evolution or the presence of magnetic activity cycles (Marsh et al. 2017; Papageorgiou et al. 2018; Hu et al. 2020). Radial velocities for the brighter systems, combined with archival multiband photometric LCs, will provide accurate absolute parameters for the components of the systems. For the total eclipsing contact system candidates (e.g., ZTF18aabexej and ZTF18aabfljv), approximate absolute physical parameters could be derived from multiband photometry only (Papageorgiou et al. 2023). Using their absolute physical parameters, it is possible to evaluate the ability of the stellar engine to produce period and LC variations through the Applegate mechanism (Applegate 1992).

We find that, although none of our systems are particularly hot, systems with a negative O’Connell effect (ΔI_sG < 0) have higher temperatures in general. This could mean that in these systems there is a correlation between the O’Connell effect and spots. Systems with negative ΔI_sG also have temperatures T < 6000 K (Knote et al. 2022). None of our systems have temperatures greater than 6000 K, likely due to the small sample size in our study.

These results also show that ΔI_sG is anticorrelated with temperature and period, and that this anticorrelation is stronger when compared with |ΔI_sG|. These correlations were found in all the different passbands, and they are consistent for most of them, with the exception of r and its corresponding correlation between |ΔI_sG| and the period, which is different from all the anticorrelations from other passbands. This might be due to the small sample size or the fact that positive and negative O’Connell effects follow different correlations, as seen in Table 4, where we see a significant difference when comparing negative and positive ΔI_sG in the r band. The same seems to happen in g, as the corresponding correlation with |ΔI_sG| is nonexistent. This, combined with the fact that c and o follow only weak anticorrelations, suggests that more data are needed to put our results on a firmer footing. These few differences between correlations for positive and negative ΔI_sG seem to indicate that the results from this study are not fully consistent with Kouzuma (2019)’s study on starspots.

In all cases, correlations were much more significant for ΔI_sG than for |ΔI_sG|, even though the correlation with the period seems to be weaker in passbands o and c. Although the difference is not very large, more data in these two passbands would be useful to obtain more conclusive results. This is also not related to the anticorrelation between ΔI_sG and the temperature, where all passbands correlate relatively similarly with no exceptions, albeit these correlations do seem to become weaker toward longer wavelengths. Observations of a larger sample, particularly at short wavelengths, would certainly be of help in establishing these correlations.

Furthermore, while our targets are intrinsically faint, some of the scatter in the OER, ΔI_sG, and LCA parameters may be due to time variability. The evolution and migration of starspots over time can contribute to this scatter, as evidenced by the dispersion seen in our phase-folded LCs. When the data for our stars were split into two or more segments covering smaller time spans to check for long-term patterns, in almost all cases the LC segments were incomplete, making further analysis unreliable. Thus, although our current dataset covers many cycles, a more detailed temporal analysis with higher-cadence data would help us disentangle intrinsic long-term variability patterns from asymmetry effects that are (nearly) invariant in time.

From this, at least for the stars in our sample, we cannot rely solely on the idea of starspots as the only cause of the O’Connell effect, suggesting that multiple causes may be at play. This is supported by the fact that we observe differences in the correlations between positive and negative values of ΔI_sG and the period and temperature. Although the results for different passbands are not fully conclusive, we also find that, for the majority of our sample, the O’Connell effect is most predominant at short wavelengths. Cool starspots have a stronger contrast against the stellar photosphere in bluer wavelengths, which may contribute to the pronounced effect observed in these bands. Still, starspots might not provide a sufficient explanation for all O’Connell effect systems. A larger sample, and an extension of our coverage toward shorter effective wavelengths, would help in further investigating these results.

Starspots are typically the reason given for the O’Connell effect. While this could be true for some systems, we find evidence that, at least for the stars in our sample, starspots may not be its only cause. If so, the O’Connell effect may also be connected to the interaction between the components rather than, or in addition to, starspots. Extending the size of the samples studied using a multiwavelength approach would be of crucial importance in unveiling the root cause(s) of this intriguingphenomenon.

Acknowledgments

We thank an anonymous referee for several useful comments and suggestions. Support for this project is provided by ANID’s FONDECYT Regular grants #1171273 and 1231637; ANID’s Millennium Science Initiative through grants ICN12_009 and AIM23-0001, awarded to the Millennium Institute of Astrophysics (MAS); and ANID’s Basal project FB210003. AP gratefully acknowledges the support provided by the grant co-financed by Greece and the European Union (European Social Fund – ESF) through the Operational Programme “Human Resources Development, Education and Lifelong Learning” in the context of the project “Reinforcement of Postdoctoral Researchers – 2nd Cycle” (MIS-5033021), implemented by the State Scholarships Foundation (IKY). This work has made use of data from the Asteroid Terrestrial-impact Last Alert System (ATLAS) project. ATLAS is primarily funded to search for near earth asteroids through NASA grants NN12AR55G, 80NSSC18K0284, and 80NSSC18K1575; byproducts of the NEO search include images and catalogues from the survey area. The ATLAS science products have been made possible through the contributions of the University of Hawaii Institute for Astronomy, the Queen’s University Belfast, the Space Telescope Science Institute, the South African Astronomical Observatory (SAAO), and the Millennium Institute of Astrophysics (MAS), Chile.

References

Appendix A: Model comparison

Here we provide a brief comparison between our model and the 12-term Fourier fit proposed by Wilsey & Beaky (2009). The latter method can perform equally or better than the one proposed in this study when well-sampled data are available, as in the case of our ZTF LCs (see Fig. A.1). However, in the case of noisier data, as happens with our ATLAS LCs, the 12-term Fourier fit may introduce artifacts, as shown in Fig. A.2, whereas our model remains well behaved.³ In most cases, the primary problem is to get a good fit around ϕ = 0.5. While this may not be relevant when measuring ΔI_max, it is crucial to calculate accurate OER and LCA values. The model presented in this work tries to work around that problem by using different sG fits around the different eclipses (θ = 0, 0.5, 1), which ensures that the minima of each LC will be accurately reproduced. In our case, the Fourier fit is then simply used to correct the sum of all three different models and get accurate maxima.

Although our model is more flexible and can be used as an alternative phenomenological description of the LCs of EBs presenting the O’Connell effect, one of the advantages of using a 12-term Fourier fit is that it allows for an early classification of the system using its cosine coefficients, specifically the second and fourth terms, a₂ and a₄ (Rucinski 1997). This can be of great utility for unclassified systems. We defer finding such correlations with all different parameters that enter our phenomenological model to future work, in which we plan to use an enlarged sample containing previously classified O’Connell effect EBs.

Fig. A.1. Comparison between a 12-term Fourier fit (red) and the model proposed in this work (black), using ZTF g-band data. The reduced χ2 values for the Fourier and new models are 7.9 and 8.2, respectively.

Fig. A.2. Same as Fig. A.1 but for the o band. The Fourier fit is shown in blue. The reduced χ2 values are 19.9 and 13.8, respectively.

Appendix B: Comparison of the maxima

Although our model performs better in some cases when calculating O’Connell LC features like OER and LCA, this is not necessarily true for ΔI_max. For this reason, we carried out a direct comparison between the two methods used in this work, namely the polynomial fit and our model (ΔI_poly and ΔI_sG, respectively, both of which are shown in Table 2). As seen in Fig. B.1, for the ZTF filters the values are sufficiently similar that we cannot distinguish between their results. However, when ATLAS data are used, we do see some important discrepancies in a number of cases; an example of this is shown in Fig. B.2, where we can see a noticeable difference in the inferred secondary maxima. Through visual inspection, we conclude that the maxima obtained using the sG model (red crosses in Fig. B.2) are more reliable in the majority of cases since.

Fig. B.1. Comparison of |ΔImax|. Top panel: ZTF g and r filters. Bottom panel: ATLAS c and o filters. In both cases the dashed line corresponds to the identity line.

Fig. B.2. Comparison between the polynomial fit and our proposed model. Top panel: Complete LC of ZTF17aaawfrc in ATLAS c. Bottom-left panel: Primary maxima of the LC. Bottom-right panel: Secondary maxima of the LC. The solid black line corresponds to the sG model, and the dashed line to the polynomial fit. Blue and red crosses indicate the maxima obtained using the polynomial fit and the sG model, respectively.

Appendix C: Modeled light curves

Here we present all the LCs used in this work with each corresponding model described in Sect. 4 for all the different passbands.

Fig. C.1. LCs of the O’Connell EB candidate ZTF17aaadztm in the four filters used in this study: ZTF r, ZTF g, ATLAS c, and ATLAS o. The gray line corresponds to the best fit using our model (Eq. 5). The period of the star is also shown at the top of the figure, along with its ZTF ID.

Fig. C.2. Same as Fig. C.1 but for ZTF17aaaesno.

Fig. C.3. Same as Fig. C.1 but for ZTF17aabvxoo.

Fig. C.4. Same as Fig. C.1 but for ZTF18aaacuaj.

Fig. C.5. Same as Fig. C.1 but for ZTF18aaaclzz.

Fig. C.6. Same as Fig. C.1 but for ZTF18aaakvwc.

Fig. C.7. Same as Fig. C.1 but for ZTF18aabexej.

Fig. C.8. Same as Fig. C.1 but for ZTF18aabfljv.

Fig. C.9. Same as Fig. C.1 but for ZTF18abaaway.

Fig. C.10. Same as Fig. C.1 but for ZTF18aagrehv.

Fig. C.11. Same as Fig. C.1 but for ZTF17aaawfrc.

Fig. C.12. Same as Fig. C.1 but for ZTF17aabumbt.

Fig. C.13. Same as Fig. C.1 but for ZTF18aabfkyu.

Fig. C.14. Same as Fig. C.1 but for ZTF18aapmolr.

All Tables

All Figures

	Fig. 1. Phenomenological modeling of the O’Connell EB candidate ZTF18aabexej. Top panels: Four components of our model, as described in the text (Sect. 4). Bottom panel: final model fit to ZTF g-band data.
In the text

	Fig. 2. Relation between the OER and ΔIsG for the 14 EBs in our sample. Different symbols and colors denote different filters, namely r, g, c, and o, following the key at the top right of the plot. The dashed line represents a least-squares fit to the data for all filters combined. The coefficient of determination is r2 = 0.831, and the resulting fit is given by OER = −2.659 ΔIsG + 1.0.
In the text

	Fig. 3. Same as Fig. 2 but for the relation between the OER and LCA. The dashed line represents the OER = 1 locus.
In the text

	Fig. 4. Absolute values of the difference between the maxima for all the objects in each of the four filters used in this study, plotted as a function of their effective wavelengths. Symbols have the same meaning as in Fig. 2.
In the text

	Fig. 5. Absolute mean values of the difference between the maxima for all the objects in each of the four filters used in this study, plotted as a function of their effective wavelengths. Symbols have the same meaning as in Fig. 2.
In the text

	Fig. 6. Effective temperatures of our sample stars plotted as a function of ΔIsG (upper panel) and |ΔIsG| (bottom panel). Symbols have the same meaning as in Fig. 2. Gray corresponds to data from Knote et al. (2022).
In the text

	Fig. 7. Correlation between ΔIsG and the period for our sample. Symbols have the same meaning as in Fig. 2. Gray corresponds to data from Knote et al. (2022), limited to the range 0.2–1 days for clarity.
In the text

	Fig. A.1. Comparison between a 12-term Fourier fit (red) and the model proposed in this work (black), using ZTF g-band data. The reduced χ2 values for the Fourier and new models are 7.9 and 8.2, respectively.
In the text

	Fig. A.2. Same as Fig. A.1 but for the o band. The Fourier fit is shown in blue. The reduced χ2 values are 19.9 and 13.8, respectively.
In the text

	Fig. B.1. Comparison of |ΔImax|. Top panel: ZTF g and r filters. Bottom panel: ATLAS c and o filters. In both cases the dashed line corresponds to the identity line.
In the text

Fig. B.2. Comparison between the polynomial fit and our proposed model. Top panel: Complete LC of ZTF17aaawfrc in ATLAS c. Bottom-left panel: Primary maxima of the LC. Bottom-right panel: Secondary maxima of the LC. The solid black line corresponds to the sG model, and the dashed line to the polynomial fit. Blue and red crosses indicate the maxima obtained using the polynomial fit and the sG model, respectively.

In the text

	Fig. C.1. LCs of the O’Connell EB candidate ZTF17aaadztm in the four filters used in this study: ZTF r, ZTF g, ATLAS c, and ATLAS o. The gray line corresponds to the best fit using our model (Eq. 5). The period of the star is also shown at the top of the figure, along with its ZTF ID.
In the text

	Fig. C.2. Same as Fig. C.1 but for ZTF17aaaesno.
In the text

	Fig. C.3. Same as Fig. C.1 but for ZTF17aabvxoo.
In the text

	Fig. C.4. Same as Fig. C.1 but for ZTF18aaacuaj.
In the text

	Fig. C.5. Same as Fig. C.1 but for ZTF18aaaclzz.
In the text

	Fig. C.6. Same as Fig. C.1 but for ZTF18aaakvwc.
In the text

	Fig. C.7. Same as Fig. C.1 but for ZTF18aabexej.
In the text

	Fig. C.8. Same as Fig. C.1 but for ZTF18aabfljv.
In the text

	Fig. C.9. Same as Fig. C.1 but for ZTF18abaaway.
In the text

	Fig. C.10. Same as Fig. C.1 but for ZTF18aagrehv.
In the text

	Fig. C.11. Same as Fig. C.1 but for ZTF17aaawfrc.
In the text

	Fig. C.12. Same as Fig. C.1 but for ZTF17aabumbt.
In the text

	Fig. C.13. Same as Fig. C.1 but for ZTF18aabfkyu.
In the text

	Fig. C.14. Same as Fig. C.1 but for ZTF18aapmolr.
In the text