A star under multiple influences Magnetic activity in V815 Her, a compact 2+2 hierarchical system ⋆

Context. Close binaries with magnetically active components are astrophysical laboratories for studying the e ff ects of binarity on activity. Of particular interest are binary and multiple star systems that contain a solar-type active component with an internal structure similar to the Sun, allowing us to study how the dynamo of a solar-type star would work under di ff erent conditions. Aims. We are conducting a comprehensive investigation of V815Her using photometric and spectroscopic data to understand the origin of the activity and what influences it in the short and long term. Methods. Using space photometry we performed light curve modeling in order to derive astrophysical and orbital parameters for the eclipsing binary subsystem V815Her B. Using archival photometric data covering a century we carried out a time frequency analysis.Spectral synthesis was applied to determine the basic astrophysical parameters of the rapidly rotating primary using high-resolution STELLA spectra recorded in 2018. Results. Photometric analysis of archived data revealed multiple cycles on timescales between ∼ 6.5 and ∼ 26 years, some of which may be harmonic. From TESS photometry we obtained orbital solution for the V815Her B subsystem. By placing the primary component on the Hertzsprung − Russell-diagram, we can deduce an age of ≈ 30Myr, in line with the high Li-6707 abundance. The STELLA spectra covering the 200 day-long observing season enabled to create 19 time-series Doppler images, which revealed a constantly changing spotted surface on a time scale of a few weeks. From the consecutive image pairs we built up the average cross-correlation function map to measure the surface di ff erential rotation of the spotted star, from which we derive a weak solar-type surface shear. Conclusions. We found evidence that the V815Her B component previously apostrophized as a ’third body’ is actually an eclipsing close binary subsystem of two M dwarfs with a period of 0.5d, i.e., V815Her is a 2 + 2 hierarchical quadruple system. The system is apparently young, only a few times ten million years old, consistent with the spotted primary V815Her Aa being a zero-age main sequence star. Spot activity on the primary was found to be vivid. Fast starspot decay suggests that convective-turbulent erosion plays a more significant role in such a rapidly rotating star. The weak surface shear of V815Her Aa due to di ff erential rotation is presumably confined by tidal forces of the close companion V815Her Ab. The slowly increasing photometric cycle of about 6.5 years on average is interpreted as a spot cycle of V815Her Aa, which is probably modulated by the eccentric wide orbit.


Introduction
Close binaries that contain magnetically active components may be regarded as astrophysical laboratories for studying the effect of binarity on activity.Tidal forces in such systems are supposed to modify the operation of the magnetic dynamo, by altering the flux emergence patterns at the surface (Holzwarth & Schüssler 2003a,b).Amplified dynamo operation is expected, for instance, ⋆ Based on data obtained with the STELLA robotic observatory in Tenerife, an AIP facility jointly operated by AIP and the Instituto de Astrofisica de Canarias.
when fast rotation is sustained by weakening magnetic braking through angular momentum transport in a close binary system (Song et al. 2018).Although, other binarity-related effects like early-evolution accretion processes can also result in faster rotation (Meibom et al. 2007).In all probability, the gravitational influence of a close companion can suppress the differential rotation of the convection zone (Scharlemann 1982;Collier Cameron 2007) or even force the formation of active longitudes at certain phases fixed to the orbit (Oláh 2006;Kővári et al. 2021).Especially interesting are those binaries and multiple star systems which contain active G-type main-sequence (MS) stars with inner structure (convective zone) similar to that of the Sun, making it possible to investigate how the solar-type dynamo would work under different circumstances.
In accordance with the above, this time we choose the active G dwarf component of the single-lined but in fact a quaternary star system V815 Her (=HD 166181) for our study.In the beginning, it was believed that the G star with its enhanced magnetic activity inferred from strong Ca ii H&K emission is the primary member of a single-lined binary (SB1) system with an orbital period of ≈1.81 days (Nadal et al. 1974).Later using new spectroscopic data taken in 1992 at Kitt Peak National Observatory Dempsey et al. (1996) found a significant radial velocity difference of 11 kms −1 compared to the data of Nadal et al. (1974), which was attributed to a possible third body with an orbital period of many years.This idea was warmed up by Fekel (2004) who, using new observations from 2002, calculated long-period orbital elements for the outer companion in the system for the first time with a preliminary period of 6.3 yr.With additional radial velocity measurements from 2003, this value was soon refined to 2092.2±5.8 d (Fekel et al. 2005).From their mass estimates for the system components of the long-period orbit Fekel et al. (2005) also suggested that the unseen component might also be a close binary.An important step forward in the story was that the outer component has been resolved by direct imaging within the frame of Gemini Deep Planet Survey (GDPS) by Lafrenière et al. (2007) who confirmed the astrometric solution of Fekel et al. (2005) for the third body.Going even further, in our paper we present space photometry from NASA's Transiting Exoplanet Survey Satellite (TESS) to demonstrate that the suspected third body itself is indeed an eclipsing close binary system of two unseen but very likely M dwarfs, i.e.V815 Her is in fact a four-star system consisting of two close binaries orbiting each other.
According to photometric observations, it has long been known that the G (most probably G5-G6) active component features cool spots on its surface (Mekkadan et al. 1980).Strassmeier et al. (1989) reported automatic photoelectric telescope (APT) data from 1984-1986 in U, B and V colours and found that the photometric period associated to the rotation of the G star was 1.8195 d and 1.8356 d for 1985 and 1986, respectively, i.e. ≈1% longer than the orbital period from Nadal et al. (1974).From an extended study using long-term APT data taken in BV between 1984and 1998Jetsu et al. (2000) found seasonal photometric period changes attributable to surface differential rotation of the spotted star.In addition, their time-series analysis revealed that one dominant active longitude was present for the full 14 year-long observing period while another weaker one emerged during some of the seasons.This APT data were re-analyzed in Savanov (2009) to reconstruct temperature inhomogeneities of the stellar surface.Again, two active longitudes were found, separated by about 0.5 rotation phase (cf.Jetsu et al. 2000), showing the so-called flip-flop phenomenon, i.e. when dominance switches quasi-periodically from one active longitude to the other.
Like the photosphere, the chromosphere and corona of the G component show the signatures of strong magnetic activity over a wide range of the electromagnetic spectrum.In addition to the strong and narrow Ca ii H&K emission (Nadal et al. 1974;Bopp 1984) UV observations from the IUE satellite (Fekel et al. 1986) also indicate an active chromosphere.The absolute luminosity in the Mg ii k line supports this enhanced activity (Cardini et al. 2003), however, the lack of hot component in the UV spectra suggests that the companions without contribution should be M-dwarfs or fainter.Moreover, the G star performs extreme-ultraviolet (EUV) radiation measured during the ROSAT Wide Field Camera All-Sky Survey (Pounds et al. 1993) and by the Extreme Ultraviolet Explorer satellite (Malina et al. 1994).Due to the ROSAT mission, V815 Her is known for its bright X-ray corona (Dempsey et al. 1993a,b).Indeed, its coronal luminosity of L X =3.183 10 30 erg (Makarov 2003) classifies V815 Her as one of the 100 brightest stellar X-ray sources within 50 pc.Accordingly, with such an active corona, the star is known as a radio source as well (Drake et al. 1992).
V815 Her is not associated with any particular moving group (but see Wichmann et al. 2003, who cited the star as Pleiadeslike).The strong Li line in the optical spectrum (Fekel et al. 2005) indicates that the G component is indeed young, near the zero-age main sequence (ZAMS).This observation is consistent with the star being listed in the Two Micron All Sky Survey Point Source Catalog (2MASS/PSC, Haakonsen & Rutledge 2009) due to its infrared excess, most probably related to primordial protoplanetary material.
In our paper, we perform a Doppler imaging study for the G star for the first time with reconstructing 19 time-series Doppler images using high-resolution optical spectra taken during a 9month long observing run in 2018.The paper is organized as follows.In Sect. 2 we present the analysis of photometric data.Using TESS observations we analyze the newly discovered eclipses and give an orbital solution for the unseen, but most likely M+M dwarf binary companion V815 Her B in the quadruple system.In Sect. 3 we derive precise astrophysical parameters of the active G component for the Doppler imaging study presented in Sect. 4. Results are discussed in Sect.5, and summarized in Sect.6.To prepare data for the computation of the orbital solution of this new close binary we fitted the spotted light curve of the G star with a time-series three-spot model (see Fig. 1), which could properly follow the ever-changing light variation instead of a constant rotational period and its harmonics.For the fit we used an upgraded version of the SpotModeL code originally presented in Ribárik et al. (2003).For the light curve modeling the fixed parameters used were only crudely approximated (temperature, limb darkening, etc.) since the main goal was only to fit the light curve as smoothly as possible.The data were also cleaned from obvious outliers and a possible flare.The resulting model was subtracted from the original data and binned to 0.005 days (7.2 minutes) which was then used for binary modeling.For the eclipsing binary (EB) light curve modeling the residual curves of the three TESS sectors were converted into the flux domain, and normalized to unity.Then they were phase-folded with the previously determined eclipsing period of 0 d .5206-dayand binned into 500 evenly spaced phase cells.The fluxes in each cell were simply averaged, and these 500 mean flux values were rendered to the mid-phase moments of each cell.In such a manner we obtained evenly spaced, folded, binned, mean light curves separately for sectors 26, 40 and 53, which are plotted in Fig. 2.

Photometric data and analysis
The EB light curve analyses of these prepared data sets were carried out with software package Lightcurvefactory (see Borkovits et al. 2019, 2020, andfurther references therein).We analysed the data of the three sectors separately.During our MCMC-based parameter search process seven light curve parameters were adjusted, as follows: the total duration of an eclipse (which is a direct observable, and closely relates to the sum of the fractional radii of the two stars); the R Bb /R Ba ratio of the stellar radii; the T Ba effective temperature of the primary; the logarithm of the ratio of the effective temperatures log(T Bb /T Ba ); the "third" light ℓ, which allows us to take into account the light contamination caused by the much brighter binary V815 Her A; the EB's inclination i B and, finally, the T 0 time of the reference epoch. 2 Regarding the other orbital elements, the folded mean light curve in Fig. 2 reveals that the secondary eclipses occur just at (or, at least, very close to) phase 0.5, i.e. in the mid-times between the consecutive primary eclipses and, moreover, their durations look similar than that of the primary ones.Thus, we assume a circular EB orbit and, hence, keep the eccentricity fixed to e B = 0.The original orbital solutions without assuming surface inhomogeneities are drawn with grey lines, while the models drawn with red lines fit also the light curve variations arising from stellar spots.The residuals of both fits are plotted on the bottom parts of all panels.Note, that the eclipses in sectors 40 and 53 (middle and bottom panels) shifted from phase 0.0 and 0.5 indicating timing variations which is discussed in the text in detail.
The adjustment of the effective temperature of the primary T Ba needs some further explanation.In general, in the case of a single photometric band EB light curve, the temperature ratio of the two stars is well determined, while the effective temperatures of each star, themselves, play only a minor role in forming the light curve and hence, they are ill-determined.So, the usual treatment is to take one of the temperatures from some outer sources (e.g.spectroscopy or SED analysis), and to keep it fixed.However, in the present case, T Ba was treated as a free parameter in order to check the physical reliability and consis-tency of our solution.This was done as follows.We know that the brighter binary should be a very young object.We reasonably assume that the two binaries are coeval, and thus, the two stars of binary B locate also on (or, very near to) the zero-age main sequence.Hence, we expect that the empirical mass-luminosity and mass-radius relations of Tout et al. (1996) are applicable for these stars.Therefore, with the combined and inverted use of these relations, we can compute the stellar masses and the radii, as well, from the effective temperatures.Then we can compare the total mass of binary B, obtained from our analysis to the minimum mass inferred from the RV solution of Fekel et al. (2005).Moreover, we can compare the stellar radii, calculated directly through the M − R relations of Tout et al. (1996) to that calculated from the combination of the fractional radii (which are direct outputs of the light curve analysis) and the masses, inferred again from the Tout et al. (1996) relations.
Finally, we also note, that a first inspection of the folded, mean light curve of binary B reveals that despite the averaging over ∼25 orbital cycles, the out-of-eclipse light variations did not average out perfectly.One possible origin of these out-of-eclipse light variations might be the surface brightness inhomogeneities (e.g.stellar spots), which, for the (almost) synchronized stellar rotations may have 'survived' the averaging for the orbital period.We model these variations together with the EB light curve modeling simultaneously in a purely mathematical way, i.e. independent of their real origin, as follows.In each trial step, after the removal of the model EB light curve from the input light curve data, the residual is modelled with harmonic functions of two fixed frequencies (namely, the twice of the orbital frequency and, the orbital frequency itself), of which the four (plus one) coefficients are obtained via matrix inversion.Then, this mathematical model of the residual light curve is added to the EB model light curve and the actual χ 2 value is calculated for this mixed model curve.
In Table 1 we tabulate the median values and the 1σ statistical uncertainties of the parameters (including several derived, physical ones) that were obtained from our analysis for the sector 26 data.We also display the synthetic model fits (both with and without the mathematically modelled out-of-eclipse light variations) derived from the best-fit solution in Fig. 2. We note that in case of the analysis of sector 40 and 53 light curves we obtained mostly similar parameters within the tabulated 1-σ uncertainties and, hence, we do not tabulate these results separately.The two exceptions are the ℓ extra flux and the T 0 epoch of the primary minimum.Regarding the former one, in sectors 40 and 53 we obtained smaller ℓ parameters with ∼5−10%, which we explain with a smaller level of contaminated fluxes in the TESS pixels.What physically more interesting is, the shift of the phase of the primary eclipses, which can also be seen in the consecutive panels of Fig. 2. In order to understand the nature and origin of this effect, we carried out a detailed analysis of the times of eclipsing minima.

Eclipse timing variation analysis of V815 Her B
We determined mid-eclipse times for all the observed primary and secondary eclipses of V815 Her B for the available TESS cleaned time-series (see in Sect.2.1).The individual times of minima are tabulated in Table A.1 of Appendix A. We derived an eclipse timing variation (ETV) curve using the linear ephemeris where the zero epoch is given in barycentric Julian Days (BJD).We plot the obtained ETV curve in Fig. 3.As one can see, the ETV curve displays clearly a non-linear behaviour, which reveals that the observed eclipsing period is varying.
As a next step, we demonstrate that this variation can be explained by the light-travel time effect (LTTE) caused by the revolution of the eclipsing binary V815 Her B around its brighter counterpart V815 Her A, and the parameters of this LTTE orbit is in accordance with the former RV solution of Fekel et al. (2005).This indicates, that this currently discovered EB is identical with the third component suspected formerly through RV measurements.In order to demonstrate this fact, we derive LTTE solution from the outer RV orbit of Fekel et al. (2005, their Table 2).While doing so, one should be cautious about a few things.First, the argument of pericenter of the outer orbit (ω A ) given in the RV solution refers to the orbit of component A, while in the LTTE solution of component B, the argument of pericenter of component B (ω B = ω A + 180 o ) is present.Moreover, the projected semi-major axis a A sin i out of the orbit of component A around the center of mass of the whole quadruple system, which was deduced from the amplitude of the RV curve, should also be transformed to the projected semi-major axis of component B (a B sin i out = m A /m B × a A sin i out ).For this step, we must know the mass ratio of the two binaries.For the formerly known, brighter component we accept m A = m Aa + m Ab = 1.27M⊙ (Fekel et al. 2005), while our solution gives the mass of component B to be m B = m Ba + m Bb = 0.63 M⊙, which results in an outer mass ratio of m A /m B = 2.02.Taking into account these constraints, and plotting the theoretical ETV curve against the measured ones, one can see a nice agreement (see Fig. Table 2. Comparison of the periods and geometry of the wide (AB) and the two close orbits (A and B) in the 2+2 quadruple system V815 Her.
.17 in our readings demonstrate nicely that the spectroscopic multiplicity of V815 Her is confirmed independently through TESS photometry, and the third spectroscopic component itself is a close, eclipsing binary, forming such a way a close quadruple system with a 2+2 hierarchy.
For comparison, the periods and geometrical parameters of the wide (V815 Her AB) and the two close orbits (V815 Her A and B) in the 2+2 hierarchical system are given in Table 2.

Long-term photometric behavior using archival data
V815 Her has data from scanned photographic plates in the Digital Access to a Sky Century @ Harvard (DASCH, Grindlay et al. 2009) database for about hundred years between 1890-1990 supplemented with photoelectric data collected in Jetsu et al. (2000) with an overlap of 5 years.We stitched the two data sets together by shifting the photoelectric B band light curve to match the median of the DASCH light curve in the 5 year overlap.These data cannot depict the rotational modulation of the G star, but can be used to search for long-term variability and trends.We combine the photographic and photoelectric data sets in Fig. 4.
Looking for possible activity cycles we used our Time Frequency Analyzer package TiFrAn (Kolláth & Oláh 2009) with Choi-Williams distribution kernel which gives good resolution in the frequency domain and somewhat less in the time-domain.
The data were averaged in bins of 365 days.Although there is a huge difference between the observational noise of the photographic and photoelectric data, the magnitude of the noise does not alter the frequency pattern.The different kernels of TiFrAn and their features are discussed in detail in Kolláth & Oláh (2009).
The data and results are plotted in Fig. 4. We reveal a cycle with a slowly increasing period, which has a steady amplitude at ∼6.5 years.Two other cycles of ∼9.1 and ∼13 years with similarly strong amplitudes are present throughout the observations.Between them, we also find an ∼11.5-year cycle which is still clearly visible despite its weaker amplitude.This fits quite well to twice the period of the outer orbit, i.e. 11.45 years (=4184 d); see the lower dashed line in Fig. 4. The longest persistent feature (ignoring the one corresponding to the total data length) points towards a timescale around ∼24-26 years with weakening amplitude.The cycle lengths of ∼13 and ∼26 years are perhaps harmonics of ∼6.5 years, although in contrast to the slow growth of the 6.5-year cycle over time, the 13-year cycle seems to decrease towards the end of the time series.At this point, it is difficult to interpret the possible physical background of each cycle that appears in the brightness change.In any case, we will attempt to list some ideas in Sect. 5.

STELLA-SES spectra
The spectroscopic observations were carried out with the 1.2-m STELLA-II telescope of the STELLA robotic observatory (Strassmeier et al. 2010) located at Izaña Observatory in Tenerife, Spain.STELLA-II is equipped with the fibre-fed, fixedformat STELLA Echelle Spectrograph (SES) providing an average spectral resolution of R = 55, 000 along the covered 3900-8800 Å wavelength range.Further details on the performance of the system and the data-reduction procedure can be found in Weber et al. (2008Weber et al. ( , 2012)).Altogether 621 high-resolution spectra were recorded between 07 March and 11 November, 2018.With a default exposure time of 3000 s an average signal-to-noise ratio (S/N) of ≈200:1 could be reached and only 18 spectra had to be discarded due to their poor quality.This data set were used for both the spectral synthesis detailed in the next section and the Dopplerinversions (see Sect. 4).The observing log is summarized in Table B.1.Assuming synchronized rotation for the G star in the close binary subsystem V815 Her A, the rotational phases are calculated using the orbital period from Fekel et al. (2005) according to the following equation: (2) The zero phase is the time of maximum radial velocity, i.e., the primary component is closest to us in ϕ=0.75 phase, accordingly, the λ=270 • longitude of the primary is facing the secondary.

Spectral synthesis
For the determination of some of the basic astrophysical parameters, we used the spectral synthesis code SME (Piskunov & Valenti 2017).Ten good quality (S/N∼200) STELLA-SES echelle spectra were randomly selected to carry out spectral synthesis independently on each.During the computations, local thermodynamical equilibrium (LTE) was assumed and MARCS atmospheric models (Gustafsson et al. 2008) were used.Atomic line data were taken from the VALD line database (Kupka et al. 1999).A global macroturbulence relationship was adopted from Valenti & Fischer (2005, see their Eq.1).The parameters, T eff effective temperature, log g surface gravity, [Fe/H] metallicity, and v mic microturbulence were determined iteratively; for details on the fitting method see Kriskovics et al. (2019).For the fits the equatorial rotational velocity v sin i was kept at 30 kms −1 since the Doppler-inversion (see Sect. 4 in this paper) proved to be artifact-free at this value, while using e.g.27 kms −1 and 33 kms −1 caused severely artifact-ridden results having unrealistic equatorial brightening with polar darkening or polar brightening with equatorial darkening, respectively.We note that if v sin i was treated as a free parameter as well, the process yielded only slightly different results compared to the value of 30 kms −1 .
Lithium 6707Å abundance was determined also by SME using LTE synthesis together with 3D-NLTE departure coefficients taken from Harutyunyan et al. (2018).We selected the same ten spectra of high signal-to-noise ratio as above, for which we determined the lithium abundance separately.Averaging these yielded logarithmic abundances of A(Li) LTE = 1.66 ± 0.04 and A(Li) NLTE = 1.75 ± 0.04.Here, logarithmic abundance is defined as A(Li) = log(N Li /N H ) + 12 where N Li and N H are the number densities of lithium and hydrogen, respectively.Fig. 5 shows  32.087 ± 0.127 an example fit to the observed Li i 6707 Å spectral region.In a solar-like star like the G component in the V815 Her system, Li-6707 is gradually depleting during the pre-MS evolution over a few tens of Myr (e.g.Hillenbrand 2009).According to Somers & Pinsonneault (2015), in the Pleiades, A(Li) typically decreases from the initial 3.3 to below ∼1.0 by the age of ∼100 Myr, however, highly dependent on the evolutionary models applied.Nevertheless, the correspondingly estimated age of a few tens of Myr is also consistent with the result of our spectral energy distribution analysis; see Sect.3.3.

Analysis with VOSA
We performed a spectral energy distribution (SED) synthesis using the virtual observatory (VO) SED analyzer tool VOSA (Bayo et al. 2008) to build the SED of V815 Her from the available VO catalogues and surveys (GALEX, Tycho, SLOAN/SDSS, Pan-Starrs, Gaia, 2MASS, AKARI/IRC, WISE, and IRAS).For fitting the resulting SED in Fig. 6 ATLAS9 Kurucz ODFNEW/NOVER models (Castelli & Kurucz 2003) were used with the nearest grid values of T eff =5500 K effective temperature, log g=4.5 surface gravity, and [Fe/H]=0.2metallicity.
It is noteworthy that the SED (see Fig. 6) shows an excess towards infrared (IR) wavelengths.This could be the characteristics of a very young age since pre-main sequence stars of ≲10 Myr (Williams & Cieza 2011) usually possess primordial protoplanetary disc, which causes excess in near-IR (cf.Sect.3.2).
In Fig. 7 we plotted the star on the Hertzsprung−Russelldiagram (HRD) along with theoretical isochrones and evolutionary tracks for the most suitable Geneva model (Haemmerlé et al. 2019).According to the plot, the age of the G star with M/M ⊙ ≲ 1.0 mass ratio is ≈26 million years.We note, however, that errors of age determination increase towards the zeroage main sequence (ZAMS) as the isochrones line up more and more densely.Nevertheless, this age determination is in agreement with the measured Li-6707 abundance in Fig. 5. Fig. 6.Spectral energy distribution for V815 Her generated by the VOSA SED analyzer tool (Bayo et al. 2008).The synthetic spectrum (blue line) is fitted to the archival photometry (red dots) from the available VO catalogues and surveys using ATLAS9 Kurucz ODFNEW/NOVER models with T eff =5500 K effective temperature, log g=4.5 surface gravity, and [Fe/H]=0.2metallicity.

Luminosity and size
The Gaia EDR3 parallax of π = 31.1649±0.1225mas (Gaia Collaboration et al. 2016, 2021) yields a distance of d = 32.087± 0.127 pc for V815 Her.Assuming V br ≈ 7 m .56 for the brightest V magnitude with B − V of 0 m .71 (cf.Fekel et al. 2005), and neglecting interstellar extinction yields an absolute visual magnitude of M V = 5 m .03 ± 0 m .03. Taking the BC = −0.132bolometric correction from Flower (1996) gives a bolometric magnitude of M bol = 4 m .90 ± 0 m .03.This gives L/L ⊙ = 0.87 ± 0.03 when using a value of M bol,⊙ = 4 m .74 for the Sun.
Assuming an inclination of i=75±5 • (cf.Fekel et al. 2005), the projected equatorial rotational velocity v sin i of 30 kms −1 with ±1.5 kms −1 error (see Sect. 3.2) would yield a stellar radius of R/R ⊙ = 1.1 ± 0.1.We note, that at a distance of ≈32 pc the estimated limb-darkened angular diameter of the G star of 0.346(±0.009)mas (Mid-infrared stellar Diameters and Fluxes compilation Catalogue V10, Cruzalebes et al. 2019) would yield a similar radius of ≈1.2R ⊙ .Although the ≳1.1R ⊙ size is a bit larger than predicted for a non-magnetic G5 ZAMS star, still, such inflation is indeed to be expected for a ZAMS star with enhanced surface magnetic activity (cf.Jackson et al. 2016, and  In Table 3 we summarize the astrophysical parameters of the G star determined in our study, along with the others adopted from the literature.

Imaging code iMAP
To reconstruct the surface temperature maps we use the state-ofthe-art Doppler imaging code iMAP (Carroll et al. 2012).This code carries out multi-line inversion on a list of photospheric lines.We selected 38 non-blended absorption lines from the 5049-6750 Å region with suitable line-depth, temperature sensitivity and well-defined continuum.The stellar surface is modeled on a 5 • ×5 • spherical grid.Each local line profile is computed with a full radiative solver (Carroll et al. 2008).The local line profiles are disk integrated, and the individually modeled diskintegrated lines are averaged.Atomic line data are taken from the Vienna Atomic Line Database (VALD, Kupka et al. 1999).Model atmospheres from Castelli & Kurucz (2004) are interpolated for the necessary temperature, gravity, or metallicity values.When solving the radiative transfer, local thermodynamical equilibrium (LTE) is assumed.Additional input parameters are micro-and macroturbulence, and the projected equatorial velocity (see Sect. 3).For the surface temperature reconstructions iMAP uses an iterative regularization based on a Landweber algorithm (Carroll et al. 2012).The iterative regularization has been proven to always converge on the same image solution,

Time-series Doppler images
From the STELLA-SES spectra described in Sect.3.1 we formed consecutive and independent data subsets in such a way that each of them covers the entire rotation phase in sufficient density.As a result, we created enough subsets for 19 time-series Doppler images.For the temporal distribution of the Doppler images see Table 4.The 19 time-series Doppler images of V815 Her Aa are plotted in Fig. 8, while the corresponding line-profile fits are given in Figs.C.1-C.3 in Appendix C. In order to follow the change in the spot distribution even more spectacularly, in Appendix D, we also present our Doppler images in the Mercator projection.The surface temperature maps indicate a constantly changing surface structure on a time scale of a few weeks.Even the coverage of the visible pole with a spot is not continuous: in images S05-S06 and S08-S10 and in S18 we see strong polar spots, which, for example, are almost absent in images S03, in S07 and S19.Spot temperature ranges are not constant either, the highest contrast (see e.g., S08) reaches 1800 K, while the smallest temperature contrast (S07) is only a few hundred degrees.The lifespan of the longest-lived spots can typically be estimated at 2-3 weeks.However, with few exceptions, it is difficult to clearly point out the corresponding spots in successive images.The obvious main reason for this is rapid spot evolution, in which surface differential rotation is most likely involved.What is certain is that during the observed period we do not see a permanent polar spot, nor is there any sign of the formation of an active longitude.For further discussion of surface evolution see Sect.5.4.

Surface differential rotation
The latitude dependent surface rotation can usually be inferred from the cross-correlation of two successive Doppler images (Donati & Collier Cameron 1997).However, random changes in spot distribution can easily mask the expected cross-correlation pattern of differential rotation.Therefore, for the detection of the surface shear operating in the background, we use ACCORD (e.g., Kővári et al. 2012Kővári et al. , 2015)), a technique which enables the averaging of latitudinal cross-correlations in the case of a sufficient number of pairs of Doppler images not too distant in time, this way suppressing the effect of randomness and amplifying the correlation pattern attributed to surface differential rotation.
The overall correlation pattern in the resulting average crosscorrelation map (Fig. 9) is fitted by a rotational law in the form where Ω(β) is the latitude (β) dependent angular velocity, while Ω eq is the angular velocity at the equator, and ∆Ω = Ω eq − Ω pole gives the difference between the equatorial and polar angular velocities.The dimensionless surface shear parameter α is defined as α DR = ∆Ω/Ω eq .The resulting fit indicates a weak solartype differential rotation, with Ω eq = 199.06± 0.40 • /day and ∆Ω = 1.99 ± 0.60 • /day, yielding α DR = 0.010 ± 0.003 shear parameter.In comparison, the solar pole-to-equator angular velocity difference is ∆Ω ⊙ = 4.1 • /day, i.e. twice that of V815 Her.
In the average correlation map, the stronger correlations are observed at latitudes higher than ≈30 • , this is probably due to the fact that at lower latitudes the spots have less contrast, and perhaps in connection with this, their lifetime is also shorter.According to the rotational law that best fits the average crosscorrelation pattern, it is most likely that the equatorial regions of the G star are (most) synchronized to the orbit.We note that the errors of the fitted rotation law plotted in Fig. 9 were estimated based on the error bars of the locations with the best correlations at the different latitude bands.However, this is the error of the fitting process, which does not necessarily reflect the true errors resulting from the method.

Multiplicity and dynamical evolution
The multiple nature of the V815 Her system makes it very special, as the size of the wide (i.e.AB, in other words 'outer' or 'long') orbit of the 2+2 hierarchical quadruple is relatively small, only a AB ≈1.6AU (cf.Table 2).The investigation of similar systems is also exciting because their secular gravitational evolution can have a significant impact on the future fate of the components.Long-term gravitational effects in such quadruple systems can have a serious impact on the evolution of the orbital eccentricities, since both binaries can act as a distant perturber on the Fig. 9.The average cross-correlation function map derived from 18 individual cross-correlations indicates how much longitude shift occurs at a given latitude due to surface shear during ∼6P rot (which is the average time difference between the mean-HJDs of the consecutive maps; see Table 4).
other pair.Unfortunately, the evolution of such a 2+2 hierarchy is still an under-researched topic (but see Vokrouhlický 2016).
It is certain, however, that this interaction is not only dynamic (e.g.Kozai oscillations, also known as von Zeipel-Lidov-Kozai mechanism or Lidov-Kozai cycles; for a recent review see Ito & Ohtsuka 2019), but due to the small size and large eccentricity of the wide orbit, as well as the compactness of the subsystems, tidal effects must also be taken into account (cf.Fabrycky & Tremaine 2007).
Due to the Kozai mechanism, in hierarchical triple systems a distant third body (as well as a distant close binary in a 2+2 system) may even decircularize the inner orbit (e.g.Hamers 2021, see also their references).The time scale of these Kozai cycles is in the order of P 2 out /P, i.e. few times 10 3 years.As a matter of fact, this can be seen in action on HD 74438 (Merle et al. 2022), a 2+2 hierarchical system in a highly eccentric wide orbit, where the eccentricity of one of the inner orbits was probably pumped by gravitational effects.According to this, the circularized close binaries of V815 Her will not necessarily remain so either.However, Kozai oscillations can be suppressed drastically when tidal friction enters the picture, shrinking and circularizing the inner orbit(s) on the tidal dissipation timescale (Fabrycky & Tremaine 2007).For both inner binaries, this timescale should be on the order of 10 5 years (Van Eylen et al. 2016), i.e. much shorter than ∼10 7 -year age of the Aa component and so the assumed age of the four-star system (in case the components are coeval).For this reason, we do not expect a large-amplitude Kozai oscillation in the case of V815 Her.
We note that in Merle et al. (2022) HD 74438 is reported to be one of the shortest period 2+2 quadruple systems with P ≈5.7yr for the wide orbit, coincidentally almost exactly the same as the AB orbital period of V815 Her.Indeed, along with a few compact and coplanar (or even doubly eclipsing) 2+2 four-star system (Tokovinin 2018;Borkovits et al. 2021;Kostov et al. 2023;Zasche et al. 2023), V815 Her is among the shortest period 2+2 hierarhical systems.However, considering the size of the wide orbit, the V815 Her system is much more compact compared to HD 74438, which can be a significant difference in terms of the strength of secondary (i.e.tidal) effects; see also Sect. 5.3 2020), and V815 Her Aa, Ba and Bb (present paper).Non-magnetic, solar metallicity isochrones are plotted from Baraffe et al. (2015).
The components of the low-mass eclipsing binary V815 Her B of the quadruple system are plotted on the mass-radius diagram on Fig. 10 together with other eclipsing binary components from the literature.Note that the oversize of the primary of V405 And were determined by two different methods using independent data (cf.caption of Fig. 10).In the figure the non-eclipsing component V815 Her Aa is also given for reference to the age determination.Together with V815 Her Ba+Bb three other systems are plotted in color (see the legends in Fig. 10).While the masses of both the primaries and secondaries of these four systems are very similar to each other, laying on the two sides of the full-convection limit of ≈0.35 M ⊙ (Chabrier & Baraffe 1997), the primaries show about twofold difference in radius.The fully convective secondaries are quite similar to each other in size, except for Thor-42 (Murphy et al. 2020).This detailed study of Thor-42 was not able to simultaneously model the mass, radius, temperature and luminosity of the components.
The relative oversize of low-mass magnetically active stars compared to inactive stars is a known phenomenon (cf.Mullan & MacDonald 2001).Kraus et al. (2011) suggested rotation as a third parameter of the mass-radius relation.Fast rotation, maintained by the binarity influence the magnetic field and through this, the strength of the magnetic activity of the late-type dwarfs.MacDonald & Mullan (2021) modelled Thor-42 using both magnetoconvection and starspots and through this they were able to reconcile all stellar parameters of the system, but at a different age that Murphy et al. (2020) found.This result underlines the importance of magnetic field in estimating fundamental parameters of active stars.We may say that the apparent discrepancy between the ages of the components of the V815 Her quadruple system shown in Fig. 10 could be largely due to the magnetic activity of the components which, in the same time, could also be different star-by-star.

Brightness change along the wide orbit
In the top panel of Fig. 11 we plot the V observations of V815 Her between 1984-98 from Jetsu et al. (2000).In these data we found a long-term brightness modulation, which showed an exciting coincidence with the approximate period of 5.73 yr of the outer (AB) orbit.Therefore, in the bottom panel of Fig. 11 we replot the V data along with the folded light curve phased with the period of the wide orbit taken from Fekel et al. (2005).It is interesting that near zero phase (which is the periastron of the wide orbit), the brightness decreases, while away from it the brightness gradually increases.The question arises as to whether this is just a coincidence or really a change related to the wide orbit, and if it is the latter, then what is causing it.
On the other hand, the time-frequency analysis of the photometric data on a much longer term plotted in Fig. 4 does not show a strong cycle corresponding exactly to the period of the outer orbit, instead we see a strong signal of a 6.5-year cycle with gradually increasing period, and strangely, a signal that appears almost exactly at twice the orbital period.Nevertheless, it cannot be ruled out that the slowly changing, roughly 6.5-year photometric period is the changing activity cycle of the Aa component, which is perhaps triggered along the wide orbit.According to Fekel et al. (2005) the eccentricity of the wide orbit is quite large (≈0.77), therefore gravitational effects exerted on the G star by V815 Her B change significantly along the outer orbit.Tidal interactions have indeed a potential impact on the internal structure and the turbulent/mixing processes (e.g.Toledano et al. 2007;Koenigsberger et al. 2021), ultimately, therefore, for the operation of the magnetic dynamo.This may also result in a change of the strength of spot activity of the G star, which ultimately causes the long-term change in brightness.We note here that such a possible binary-induced, i.e. tidally triggered activity was reported in Strassmeier et al. (2011).
An alternative explanation could be that the brightness change along the wide orbit is somehow related to the dust-gas material component of any kind (circumstellar or circumbinary) present in the quadruple star system.Assuming that there is a circumbinary disk around the Ba+Bb binary (consistently with the infrared surplus indicated by Fig. 6), such a disc if tilted to the plane of the wide orbit, can reflect the radiation of the G star differently, depending on their mutual position, i.e. the orbital phase.In our case, the reflected light can only be estimated cautiously in the case of a disc of unknown composition, size and geometry.However, if we estimate the contribution of the reflection from the disc to ≲1% (which is a rough upper estimate, cf.Bromley et al. 2021), then we can only explain a few percent of the long-term change in the mean brightness with reflection (say, ≲0.02 mag).Therefore, we conclude that changing spot activity results in the modulated brightness along the outer orbit, even though light reflected from any disc contributes little, if at all.

Spot decay and turbulent-driven magnetic diffusion
According to observations, the lifetime of sunspots depends on their area, based on which linear and nonlinear spot decay theories were born (e.g., Bumba 1963;Krause & Ruediger 1975;Petrovay & van Driel-Gesztelyi 1997;Hathaway & Choudhary 2008;Litvinenko & Wheatland 2015;Muraközy 2021;Forgács-Dajka et al. 2021, etc.).The underlying theories assume different mechanisms that also depend on the morphology of the sunspots/sunspot groups.When starspot lifetimes and starspot decay are discussed based on the solar analogy, the question arises as to whether it is a matter of scaling or other considerations are necessary (see e.g., Namekata et al. 2019, and their references).The average spot area on V815 Her is about a factor of ∼100 larger compared to long-lived big sunspots (see Fig. 8).
We note however, that the spatial resolution of our Doppler images does not reach the size range of sunspots, but in principle they are able to resolve surface structures corresponding to the largest sunspot groups.In this context, it is perhaps more reasonable to compare the starspots of our Doppler images with groups of sunspots instead of large individual sunspots.
As for the spot lifetime, based on our time-series Doppler images, the spotted surface of V815 Her is completely rearranged in a few steps, typically within ∼25-50 days (cf.Fig. 8 and Table 4).This is therefore the approximate upper limit of spot lifetime.At the same time, the lifespan of individual spots can vary widely.The shortest-lived spots can only be observed in one or at most two consecutive images, so their lifetime is ∼10-15 days.
We emphasize here, that it is difficult to discuss on the starspot decay in the context of sunspot theory without knowledge of the true morphology of starspots, i.e., do we see monolithic spots on Doppler images or clusters of smaller spots?It is even possible that one does not exclude the other, nevertheless, usually the limited spatial resolution of observations does not allow us to make a clear distinction (but see Järvinen et al. 2018).Moreover, spot decay may strongly be influenced by surface and subsurface flows.Işik et al. (2007) modeled the temporal decay of the surface magnetic flux when assuming different spot configurations while large-scale flows (differential rotation, meridional circulation) are switched on and off.Their results indicate that spot conglomerates decay faster compared to large monolithic spots, and flows also play an important role in speeding up the decay.We note that the ∆Ω ≈ 2 • /d surface shear on V815 Her is much weaker than on the Sun (see Sect. 4.3).If we take into account the t l lap time, i.e. the time taken for the equator to lap the poles, then in the case of V815 Her we get 360 • /∆Ω=180 d (i.e. about twice the solar lap time of ∼88 d).This is much longer than the observed spot lifetimes on V815 Her, therefore, differential rotation is likely to play a minor role (if at all) in spot decay.How-ever, even with flows, the simulations in Işik et al. (2007) generally predict longer spot lifetimes than what we see on V815 Her.Bradshaw & Hartigan (2014) argued that, assuming the same physical background and starting from the Gnevyshev−Waldmeier (hereafter GW) rule, the decay of sunspots and starspots can be scaled by the value of magnetic diffusivity through the turbulent scale length (which is connected to the supergranulation on the Sun).They suggested an anomalous turbulent-driven magnetic diffusion, which can account for the decay of either sunspots or starspots by scaling the characteristic supergranule/diffusion lengths.
By following the empirical GW rule, i.e. assuming a linear relationship between the maximum A s spot area and its t s lifetime, based on our time-series Doppler images, taking t s ∼20 d and R s ∼10 • , equivalent of 1.3 10 5 km, as average spot lifetime and corresponding characteristic maximum size, respectively, we can estimate a characteristic value of the spot decay rate ∆A s /∆t s of about −4.0 10 4 km 2 s −1 .In Appendix F, we also present a more detailed study to measure the spot decay rate.Although the method is not mature enough yet, it is certainly more quantitative than the above estimate.We note that this detailed study yields even faster decay of about −7.0 10 4 km 2 s −1 .Now, using t s =20 d and R s =1.3 10 5 km (see above), we get from t s ∼ R 2 s /η T (Bradshaw & Hartigan 2014, their Eq.3) a value of η T ≈ 10 4 km 2 s −1 for the turbulent diffusivity.According to Chae et al. (2008, their Eq. 19), from this value would yield a turbulent scale length of l=4.2 10 5 km or ≈0.55 R * , i.e., in the order of the stellar radius.Indeed, if we assumed an even faster spot decay from Appendix F, we would get about twice this amount.This large turbulent scale length is probably an overestimate, which may result from an overestimation of the maximum spot size.However, if, for example, several smaller spots were assumed instead of a monolithic spot, the turbulent scale length could also be reduced to a more plausible value (cf.Bradshaw & Hartigan 2014).Finally, we should also note that our measurement of the spot decay rate can be significantly distorted if, in addition to the spot shrinking, a rapid flux emergence occurs simultaneously in the same region.This would limit the validity of the interpretation that the observed spot decay rate and turbulent diffusion are directly related.

Latitudinal and longitudinal distribution of spots
In order to make it easier to assess the distribution of spots on the time-series Doppler maps, for each temperature reconstruction we calculated the sums of the f filling factor values along the λ longitude and β latitude (we have 72 longitude and 36 latitude pixels on a 5 • ×5 • spherical grid).When calculating f λ and f β , we took into account the geometric distortion on the spherical surface according to the following formula: where T i and β i are the temperature and the latitude of the ith pixel, while the limit of the summation is the image pixel number over either the longitude (n λ =72), or the latitude (n β =33, since we do not see the star below β=−75 • ).The resulting spot filling factors for the nineteen Doppler images are shown in Appendix E. Based on the latitudinal distribution of the spots, we see that especially in the case of low spot coverage, the spots tend to appear at mid-latitudes.In these cases (e.g.S01-S02-S03-S04, S07, S12-S13, S18) a somewhat equatorially symmetric distribution is seen.On the other hand, at other times spots tend to appear at high latitudes, closer to the visible pole.We note however, that, since we cannot see the other pole and its immediate surroundings, we cannot say whether the appearance of the polar spots shows any symmetry to the equator.Granzer et al. (2000) studied the dynamics of magnetic flux tubes in young stars by numerical simulations and found that at a rotation rate of V815 Her Aa the spot emergence latitude is most likely between 25 • -50 • (see their ZAMS star model of 1 solar mass with ∼10Ω/Ω ⊙ rotation).However, similarly to V815 Her Aa, the coexistence of low-to-mid latitude spots and high latitude/polar spots on G-K type main sequence stars has already been supported by numerous Doppler imaging studies (e.g., Strassmeier 2009, and references therein).Moreover, EK Dra, the rapidly rotating young Sun, showed a very similar starspot distribution in the recent Doppler imaging study by Şenavcı et al. (2021), in which the polar and mid-latitude (>20 • ) spots are also consistent with flux emergence and transport simulations.According to Granzer et al. (2000), after flux emergence at mid latitudes a poleward flux transport by meridional circulation and/or magnetic stresses could explain the formation of polar spots (for the formation of polar spots on rapid rotators see also Schrijver & Title 2001;Işık et al. 2018).Maybe a weak sign of this can be seen in some of our subsequent Doppler image pairs (S04-S05, S13-S14, S16-S17), where the formation of the polar spot is preceded by a nearby lower latitude spot, although the available time resolution is not enough for a strict conclusion.Here, we note that the internal structure of V815 Her Aa is indeed affected by its close companion Ab, and this may be important also from the point of view of magnetic flux tube dynamics (e.g., Schuessler et al. 1996;Holzwarth & Schüssler 2003a,b).
The longitudinal spot distribution, in contrast, is more random, we see no definite indication of the existence of active longitudes.Although the changes can still mostly be followed in successive images, however, no systematic pattern can be seen in the evolution of the spotted surface over the entire period.

Surface differential rotation
From our cross-correlation study, we infer that the active G star has a weak solar-type differential rotation.Actually, this finding is in a fair agreement with the recent results obtained for two similar close binary systems V471 Tau (Kővári et al. 2021) and EI Eri (Kriskovics et al. 2023) in which the active component preserves its rapid rotation as a result of synchronization, but at the same time the degree of surface shear due to differential rotation is presumably confined by tidal interaction.These observations support the idea that differential rotation is suppressed in low-mass stars tidally locked in close binary systems (see also Collier Cameron 2007, and their references).From a theoretical point of view, this kind of behavior is not unexpected, since a marked initial differential rotation structure would converge to rigid body-like rotation due to tidal shear energy dissipation (Koenigsberger et al. 2021).Our Fig. 9 also suggests that the equatorial region of V815 Her Aa is synchronized the most with the orbital period of V815 Her A. We note that this is consistent with tidal perturbations being expected to be strongest around equatorial latitudes.

Summary
The V815 Her system undoubtedly evolves dynamically and most likely magnetically, while the two subsystems, V815 Her A and B, interact gravitationally.Our main conclusions are summarized as follows: -The wide component in V815 Her previously apostrophized as a 'third body' is actually an eclipsing close binary subsystem of two M dwarfs, i.e., V815 Her is a 2+2 hierarchical quadruple system.-V815 Her is presumably a young system, supported by the infrared surplus indicated by the spectral energy distribution in Fig. 6.The young age is also supported by the lithium abundance that has not yet been depleted.In accordance with this, the evolutionary state of V815 Her Aa is a zero-age main sequence star, with an age of ∼30M years.-A rapid spot evolution can be observed on the surface of the active G component V815 Her Aa, the spot distribution can change completely on a time range of few times ten days.-We measured a weak solar-type surface differential rotation on the G star, which is probably suppressed by the tidal force of the close companion.-Based on the long-term photometric data, we found a slowly increasing cycle length of about 6.5 years on average, which we interpret as a spot cycle.The length of the spot cycle is not necessarily closely related to the 5.73-year period of the wide orbit, but we would not rule out the possibility of this at this point.We suspect that the reason for the slow change may be the modulation due to the eccentricity of the wide orbit.
The case of V815 Her is a good example of the particular importance of tidal effects for the dynamo mechanism.For this reason, our quartet is certainly interesting for further study of how gravitational evolution of the system affects activity.

Fig. 1 .
Fig. 1.TESS sector 40 light curve of V815 Her (blue dots) fitted with a time-series three-spot model (red line).

Fig. 2 .
Fig.2.Cleaned and folded TESS Sector 26, 40 and 53 light curves (blue dots) for the wide binary companion V815 Her B. The original orbital solutions without assuming surface inhomogeneities are drawn with grey lines, while the models drawn with red lines fit also the light curve variations arising from stellar spots.The residuals of both fits are plotted on the bottom parts of all panels.Note, that the eclipses in sectors 40 and 53 (middle and bottom panels) shifted from phase 0.0 and 0.5 indicating timing variations which is discussed in the text in detail.

Fig. 3 .
Fig.3.Eclipse timing variations of V815 Her B, calculated from the mid-minima times of primary and secondary eclipses (red dots and blue squares, respectively), observed with TESS.The black line represents the theoretical LTTE curve derived from the former third-body RV solution ofFekel et al. (2005).See text for details.

Fig. 4 .
Fig. 4. Long-term photometry and time-frequency analysis of V815 Her.Top: DASCH light curve of ≈100 years for V815 Her from scanned photographic plates (black dots) supplemented with photoelectric data (blue dots) spanning more than a century between 1890-1998 with 5 years' overlap.The spline smoothed data are indicated by a red line.Bottom: Time-frequency analysis for the available photo-graphic+photoelectric data of V815 Her using the TiFrAn code (Kolláth & Oláh 2009).The plot indicates dominant cycles of different lengths and amplitudes on timescales corresponding to ∼6.5, ∼9.1, ∼11.5, ∼13 and ∼26 years, the latter two are probably multiples of the 6.5-year period.Dashed lines indicate the 2092 d period of the wide orbit and its double at 4184 d (≈11.5yr).See the text for more.

Fig. 5 .
Fig. 5.An example plot of the observed Li i 6707 Å spectral region (black) fitted with a synthetic spectrum (blue).

Fig. 7 .
Fig. 7. Position of V815 Her Aa (black dot) on the HRD plotted by using the VOSA SED analyzer tool (Bayo et al. 2008).Evolutionary tracks for the Geneva model (Haemmerlé et al. 2019) indicate an age of 30 Myr for a M≲ 1M ⊙ ZAMS star.

Fig. 8 .
Fig. 8. Nineteen time-series Doppler images of V815 Her for the 2018 observing season.From top to bottom, in chronological order, Doppler images S01-S10 are on the left, and Doppler images S11-S19 are on the right.

Fig. 11 .
Fig. 11.Top: V light curve of V815 Her fromJetsu et al. (2000) between 1984-98.Bottom: Folded light curve using the 5.73 yr (=2092.2d) period of the wide orbit.As zero phase the time of the periastron passage was set.The continuous color gradient helps to distinguish successive phases.

Table 1 .
Adjusted and derived astrophysical and orbital parameters for V815 Her B.

Table 3 .
Adopted/calculated astrophysical parameters for V815 Her Aa

Table 4 .
Temporal distribution of the subsequent data sets for each individual Doppler image.
below in this paper.

Table F .
1. Identified spots with their longitude and latitude coordinates and their calculated area.In the last two columns, the equivalent spot radius and average temperature values are provided for information.