More than half of the stars in the universe are part of binary or multiple stellar systems. However, only a fraction of these systems exhibit eclipses observable from the Earth (Guinan, 2004). Eclipsing binaries serve as fundamental building blocks in stellar astrophysics, allowing the determination of the masses and radii of the component stars through photometric and spectroscopic observations (Torres et al., 2010). Over recent decades, tens of thousands of eclipsing binaries have been identified using these methods. Notably, Slawson et al. (2011) unveiled the Second Data Release of Kepler Eclipsing Binary Stars, encapsulating 2,165 eclipsing binaries with remarkably precise light curves. Additionally, Prša et al. (2022) cataloged 4,584 eclipsing binaries identified in the initial 2-year sky survey conducted by the Transiting Exoplanet Survey Satellite (TESS) (Ricker et al., 2015). Complementarily, Qian et al. (2017), Qian et al. (2018) compiled data on 7,938 EW-type and 3,196 EA-type eclipsing binaries, as observed by the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST, also known as the Guo Shoujing Telescope), employing a low-resolution spectrograph (LRS) (Cui et al., 2012; Zhao et al., 2012). In a more recent advancement, Mowlavi et al. (2023) introduced the first Gaia catalog of eclipsing binary candidates in Gaia Data Release 3 (Gaia DR3), encompassing over 2,184,477 candidates across a brightnesses range in the Gaia G-band from a few to 20 magnitudes covering the entire sky. These efforts signify a promising advancement in the research of eclipsing binaries, demonstrating a growing interest and expanding knowledge base in this field.

Recent intensive studies have shown that many eclipsing binaries are part of more complex systems, often including three or more stars (Duchêne and Kraus, 2013). In an in-depth survey of 165 solar-type spectroscopic binaries, Tokovinin et al. (2006) examined the impact of higher-order multiplicity on orbital period distributions. Their findings highlighted notable differences between the orbital period distributions of spectroscopic binaries with and without tertiary companions, suggesting that tertiary companions might influence the primary binary’s orbital period through angular momentum exchange via the Kozai mechanism (Kozai, 1962). Furthermore, Zhou et al. (2019) suggested that dynamic interactions within these systems, particularly between the eclipsing binary and any additional companions, could alter the binary’s orbital period. Delving into specific cases offers further insights into these complex relationships. Welsh et al. (2012) disclosed that the systems Kepler-34 and Kepler-35 both harbor circumbinary planets, providing compelling examples of the diverse architectures within stellar systems. Meanwhile, Wang and Zhu (2021) raised the intriguing possibility of an undetected black hole in orbit around the eclipsing binary RR Dra, underscoring the vast range of celestial phenomena yet to be fully understood. These investigations shed light on the intricate dynamics of multiple-star systems.

The study of GZ And, a member of the visual multiple system ADS 1693, spans nearly a century. Initially, as shown in Figure 1, ADS 1693 was determined to comprise four components: ADS 1693A, B, C, and D, with GZ And (ADS 1693A) being identified as the brightest. Espin (1908) made early visual observations of the system’s variability but could not definitively attribute the brightness changes to either ADS 1693A or B. Subsequently, Walker (1973) analyzed the UBV band light curves of ADS 1693A, identifying it as a W UMa-type binary with a spectral type of G7 or K0, and confirmed that the other components—ADS 1693B, C, and D—were not variable, classifying their spectral types as G5, K4, and G7 based on their color indices. This led to the conclusion that ADS 1693 is a quintuplet system. GZ And was then designated as ADS 1693A (Kukarkin et al., 1975). Further light curve analyses by Liu et al. (1987) reconfirmed its EW-type binary classification. Walker’s (1991, 1996) measurements indicated that ADS 1693B, C, and D were positioned 8, 13, and 24 arcseconds away from A, respectively, suggesting that GZ And is part of a close triple system, with the wider pair orbiting approximately every 5.3 years. Although Liu et al. (1987) made these claims, comprehensive photometric solutions and power spectrum analyses are not available now. Later, Lu and Rucinski (1999) unveiled radial velocity curves of GZ And, indicating it as a W-subtype contact binary with a mass ratio of q = 0.514 (08), but did not detect a third component spectroscopically. While both Yang and Liu (2003) and Baran et al. (2004) provided photometric solutions for GZ And, this intriguing visual multiple system still demands further in-depth research to unravel its complexities fully.

This study unveils new findings on the W UMa-type contact binary GZ And, situated within a visual multiple system. Section 2 details the photometric data employed in our subsequent analyses. In Section 3, we analyze the O− C curve of the light minima timings to ascertain the orbital period changes of GZ And. Section 4 is dedicated to deriving the photometric solutions for the binary system by fitting its observed light curves. Finally, in Section 5, we discuss our findings and present the conclusions drawn from our analyses.

Orbital period variations are a common phenomenon in close binary systems, typically resulting from mass transfer between components, the gravitational influence of distant stellar or substellar companions, and other factors (Zhou et al., 2016; Zhou and Soonthornthum, 2019). The O-C method is commonly used to analyze orbital period variations, in which O refers to observed times of light minima and C refers to calculated times of light minima while a designated reference epoch is given. The observed data covering light minima, as shown in Figures 2 and 3, are fitted by a parabola with the least squares method to determine the observed times of light minima. Then, we compiled all previously published times of light minima and included them in the Supplementary Material, alongside the minima determined from our observations and TESS data, to analyze the orbital period changes of GZ And.

A linear ephemeris, M i n . I B J D = 2459900.121896 30 + 0 d . 3050173 × E , (1) serves as the basis for calculating the O−C (observed minus calculated) values. We used Eq. 1, where the time of the primary light minimum observed by the Xinglong 85-cm telescope is designated as the reference epoch. An initial orbital period of p = 0.3050173 days, sourced from the O-C gateway, ³ was utilized. The computed O−C values are graphically depicted in Figure 4: black dots illustrate data from the existing literature, while red and green dots correspond to findings from TESS and the Xinglong 85-cm telescope, respectively. Observing a downward parabolic trend in the O-C data, we applied a quadratic polynomial to fit the O-C curve. Due to some uncertainties not accounted for in Supplementary Table S1, we assigned equal weight to each data point during the fitting process. The resulting equation (Eqs 2, 3) is O − C = − 1.8 ± 0.8 × 1 0 − 4 − 1.78 ± 0.01 × 1 0 − 6 × E − 3.16 ± 0.03 × 1 0 − 11 × E 2 . (2) Subsequently, the ephemeris was refined to M i n . I B J D = 2459900.1217 + 0 d . 3050155 × E . (3) The orbital period of GZ And shortens at a rate of dP/dt = −7.58 (7) × 10⁻⁸ day ⋅ year ⁻¹, indicating a continuous decrease in the system’s orbital period. The parabolic fit, depicted with black solid lines, is also presented in Figure 4 for comparison. Since only a few times of light minima are obtained before BJD 2446016, and the same weight was given to all O-C values while fitting the O-C curve, the O-C values prior to −40000 are not fitted well, as shown in Figure 4.

The Wilson–Devinney (W-D) code (Wilson and Devinney, 1971; Wilson and Van Hamme, 2014) was employed to model the light curves collected using the Xinglong 85-cm telescope and TESS. The light curve morphology of GZ And is classified as EW type, which corresponds to contact binaries. Consequently, we applied mode 3, designated for overcontact binaries, in our analysis. For the star undergoing primary eclipses (commonly referred to as the primary star or star 1), we set its surface effective temperature at T ₁ = 5348 K, a value sourced from Gaia DR3 (Gaia et al., 2016; Gaia et al., 2023). Although this temperature cannot be verified at which phase of the GZ And orbit, previous studies from other researchers (Walker, 1973) and Gaia data indicate that GZ And is a solar-type contact binary. Even though the actual temperature of the primary star may be several hundred degrees higher or lower, changes to the determined parameters such as orbital inclination, mass ratio, and the temperature ratio between two components will be minimal, as they mainly rely on the variations of the light curve (Jiang et al., 2015). Given that GZ And is akin to a solar-type system, we presume the presence of convective envelopes for both stars. Accordingly, we assigned bolometric albedos of A ₁ = 0.5 and A ₂ = 0.5 (Ruciński, 1969) for the respective components. The gravity-darkening coefficients are established at g ₁ = 0.32 and g ₂ = 0.32 (Lucy, 1967), reflecting the behavior of convective stellar envelopes. For limb-darkening coefficients, we selected the monochromatic, passband-specific, and bolometric coefficients as per the values listed in van Hamme’s (1993) tables.

With the radial velocity curves for GZ And already provided by Lu and Rucinski (1999), there was no necessity to use the q-search method to determine the mass ratio (q) of the binary system. In our analysis using the W-D program, we concurrently fit Lu and Rucinski’s (1999) radial velocity curve along with the BVR _c I _c light curves from the Xinglong 85-cm optical telescope and the TESS light curve. We synchronized our study with the TESS light curve observed on 16 November 2022, matching the observation date of the Xinglong telescope’s BVR _c I _c light curves. The primary adjustable parameters in our model include orbital semi-major axis in solar radii (A), systemic radial velocity (V _γ ), star 2 average surface temperature (T ₂), the mass ratio (q) of the binary system, orbital inclination (i) observed from Earth, the modified dimensionless surface potential of star 1 (Ω₁), and passband-specific luminosities (L ₁). Notably, as illustrated in Figures 6 and 7, the observed light curves exhibit asymmetry, a characteristic often associated with the O’Connell effect (O’Connell, 1951), prevalent in contact binaries of solar or later spectral types, typically attributed to stellar magnetic activity (Kouzuma, 2019). To accommodate this asymmetry, spot parameters are integrated as variable components in our model, encompassing the spot’s latitude, longitude, angular radius, and temperature factor. Given the correlation between angular radius and temperature factor, only one is varied during the fitting process. We also explored the inclusion of the third light (l ₃) as an adjustable parameter, yet the findings confirmed the absence of tertiary luminosity contribution in GZ And. The solutions are tabulated in Table 1, with the observed and synthesized radial velocity and light curves showcased in Figures 5–7. The mass ratio listed in Table 1 is nearly the reciprocal of the value given by Lu and Rucinski (1999) since we used the opposite designation of primary and secondary stars. The geometric structure of the primary and secondary stars at phase 0 and phase 0.5, depicted in Figure 8, clearly demonstrates GZ And’s total eclipse nature, with the secondary star’s spot illustrated in dark black.

The photometric analyses presented in Table 1 characterize GZ And as a W-subtype shallow contact binary system with a filling factor of f = 14 (2) %. The system demonstrates a notably high orbital inclination of i = 89.99°, enabling the primary star’s total eclipse at orbital phase 0, vividly depicted in Figure 8. The credibility of the solutions in Table 1 is reinforced by their derivation from comprehensive radial velocity curves and fully eclipsed light curves, ensuring their dependability. Using the mass function (M ₁ + M ₂)sin ³ i = 1.763 (44)M _⊙ (Lu and Rucinski, 1999), along with the given orbital inclination of i = 89.99° and mass ratio M ₂/M ₁ = 2.07 (1), we determined the masses of the component stars: M ₁ = 0.57 (4)M _⊙ for the primary and M ₂ = 1.19 (9)M _⊙ for the secondary. By applying Kepler’s third law, we can articulate the following equation: a 3 P 2 = G 4 π 2 × M 1 + M 2 . (4) Consequently, the orbital semi-major axis of GZ And is determined to be a = 2.30 (6)R _⊙. The calculated radii and luminosities of the two-component stars are as follows: R ₁ = 0.75 (2)R _⊙, R ₂ = 1.04 (3)R _⊙, L ₁ = 0.42 (2)L _⊙, and L ₂ = 0.63 (3)L _⊙. Typically, in contact binary systems, the components, enveloped within a common atmosphere, exhibit nearly identical temperatures. Yet, the noted temperature disparity of 310K between the two stars, coupled with the modest filling factor of 14%, suggests that GZ And could represent an incipient contact binary system that has not yet attained thermal equilibrium (Zhou et al., 2015; Wang et al., 2022).

The analysis of the orbital period reveals a decreasing trend for GZ And, quantified at dP/dt = −7.58 (7) × 10⁻⁸ day ⋅ year ⁻¹. This decline is likely due to mass transfer from the more massive secondary star to the less massive primary star. Under the assumption of conservative mass transfer, the transfer rate can be calculated using the following equation: d M 2 d t = − M 1 M 2 3 P M 1 − M 2 × d P / d t . (5) Based on this, the mass transfer rate is determined to be d M 2 d t = − 9.06 ( 8 ) × 1 0 − 8 M ⊙ / y e a r . However, it is important to note that the actual mass transfer rate might significantly differ when the effects of angular momentum loss are taken into account. Furthermore, the observed parabolic variations in the orbital period, depicted in Figure 4, could potentially be attributed to the light-travel time effect, which arises due to the presence of close stellar companions within the ADS 1693 visual multiple system in which GZ And resides.

While multiple stellar systems are prevalent across our skies, the subset featuring eclipsing binaries within such systems is relatively rare (Chambliss, 1992). These configurations offer unique insights into the dynamic interactions and evolutionary processes within multiple stellar systems. Walker (1991, 1996) posited that GZ And resides within the quintuple system ADS 1693, although detailed results of this claim are pending. Gaia mission observations have provided astrometric data for this system, as documented in Gaia DR3 and summarized in Table 2. Analysis indicates that ADS 1693A (GZ And) and ADS 1693B share close proximities in distance to Earth, suggesting a potential gravitational binding between the two.

Considering that the parabolic fit in Figure 4 could be a part of a periodic variation caused by a tertiary, the distance between ADS 1693A and ADS 1693B is estimated to be approximately 1340.7AU, based on their separation of 8 arcseconds (Walker, 1991; 1996), and a distance of 167.6 pc from Earth. Using Kepler’s third law and assuming the mass of ADS 1693B is 1 M _⊙, the orbital period of ADS 1693B around ADS 1693A is estimated to be 2.95 × 10⁴ years. The period is significantly longer than the time span covered in Figure 4, reinforcing the possibility that ADS 1693A (GZ And) and ADS 1693B are gravitationally bound.

Conversely, ADS 1693C and ADS 1693D appear significantly more distant from GZ And, likely positioning them as background stars along the same line of sight rather than integral system components. An intriguing aspect is the similarity in surface gravity and metallicity between ADS 1693A and ADS 1693B, hinting at a common origin from the same molecular cloud and nearly simultaneous formation timing. Consequently, ADS 1693 emerges not as a quintuple but maybe as a physical triple system, with ADS 1693A (GZ And) as a W-subtype contact binary at its core and ADS 1693B acting as a distant companion orbiting the central binary.