BU CMi as a quadruple doubly eclipsing system

We found that the known spectroscopic binary and variable BU~CMi = HD65241 ($V$=6.4-6.7 mag, Sp~=~A0~V) is a quadruple doubly eclipsing 2+2 system. Both eclipsing binaries are detached systems moving in an eccentric orbits: pair"A"with the period $P_A$~=~$2^{d}.94$($e$=0.20) and pair"B"with the period $P_B$~=~$3^{d}.26$ ($e$=0.22). All four components have nearly equal sizes, temperatures and masses in the range $M$~=~3.1--3.4 M$_\odot$ and $A0$ spectra. We derived the mutual orbit of both pairs around the system barycenter with a period of 6.54 years and eccentricity $e$ = 0.71. We detected in pairs"A"and"B"the fast apsidal motion with the periods $U_A$~=~25.0 years and $U_B$~=~25.2 years, respectively. The orbit of each pair shows small nutation-like oscillations in periastron longitude. The age of the system estimated as 200 mln. years. The photometric parallax calculated from the found parameters coincides perfectly with the $GAIA~DR2$ $\pi$=$0.00407"\pm0.00006"$.


INTRODUCTION
The system BU CMi is designated in GCVS as an eclipsing variable EA with the period P = 2.93 d and strongly displaced secondary minimum MinII -MinI = 0.37P, Samus et al. [1], which means the large eccentricity of the orbit. So, it was included in our program of investigation of the inner structure of the stars Volkov and Volkova [2]. We started the observations in 2012 at Stará Lesná observatory, Slovakia. We managed to fix only one exit from minimum strongly apart of the GCVS ephemeris. Our further observations in the same year did not bring any result, the depth of the only observed eclipse turned out to be insignificant and the system was classified as unpromising for the inner structure investigation. But when the reviews of the bright stars such as MASCARA Burggraaff et al. [3] and TESS Stassun et al. [4] became available we returned to analysis of the system. The ASAS data Pojmanski [5] for this bright star were useless as their precision was poor, may be due to overexposure. The first glance to the Light Curve (LC) built with MASCARA data revealed that there are two groups of eclipses with different but stable periods and displaced secondary minima (see Fig. 1 with more illustrative TESS observations). Using the corrected ephemeris obtained from MASCARA observations we have resumed intensive observations of the star at the Simeiz observatory of INASAN with the 60-cm reflector and U BV -photometer constructed by I.M.Volkov, Volkov and Volkova [6]. Here we present the corrected results firstly published in Volkov et al. [7]. We managed to correct those results by disentangling far removed from the current epoch Hipparcos data, Perryman et al. [8].

OBSERVATIONS AND DATA REDUCTION
The log of our observations is presented in Table I. During observations in 2020 all four minima -two primary and two secondary ones, were recorded for both eclipsing systems observed simultaneously as a single star. HD64963 (V =8.23, A2) at a distance 22' of the variable served as a standard star for observations with arXiv:2107.09116v1 [astro-ph.SR] 19 Jul 2021 U BV -photometer equipped by photomultiplier. The observations were carried out according to the standard star -variable star -standard star scheme. The expositions in each of the photometric passbands were from 20 to 30 sec, sometimes the background near the standard and variable stars was measured. We used aperture 27.5". The maximum value of the signal from the star did not exceed 200000sec −1 in B. The signal was weaker in other photometric passbands. Such rate of the signal needs a correction due to the nonlinearity of the equipment not more than 0 m .005, which was taken into account according to the well-known nonlinearity formula for the Poisson distribution of impulses in the photon-counting flux. The dead time of our registration channel: photomultiplier + amplifier + pulse counter averaged 32.7 nanoseconds for the entire period of observations. This value has been carefully controlled, it remained stable within 5 percent. All observations were corrected for atmospheric extinction. More details of our observational method are described in Volkov et al. [9], the instrumental system is presented in Volkov and Volkova [6]. Since the number of nights suitable for observations with a singlechannel photometer is limited by weather conditions, we attempted to use CCD observations with standard and variable stars measured simultaneously. To attenuate the peak signal from such a bright star as BU CMi, images were defocused and exposures were shortened. The result was unsatisfactory. We failed to achieve the required accuracy better than 0 m .01 and the shape of the minima was distorted by systematic errors. The CCD observations in Stará Lesná observatory, Slovakia were much better. We used a Maksutov with a 15-cm aperture, so it was not necessary to attenuate the signal artificially. In addition, a larger field of view ensured the presence of a bright comparison star in one frame with a variable star. In the spring 2021 we obtained at the Simeiz observatory of INASAN only the incomplete observations in four minima due to unsatisfactory weather conditions. But these observations turned out to be very important for determination of the parameters of the external orbit of A+B pairs. BU CMi was also measured six times during observations of the SAI catalogue of bright stars ? ]. We recalculated the ultraviolet W values given in this catalogue into the standard U system, and the instrumental values B, V were corrected a little to correspond precisely to Johnson system. Two of these measurements were obtained in minima, which enable us to fix the moment of the minimum farthest from the current epoch, and four measurements turned out to be located between the minima and were averaged. Data of the brightness of BU CMi on the plateau and comparison star are given in Table II.
When analyzing the observations of Hipparcos Perryman et al. [8], we noticed that some of the observations of the object fall into minima, but, as a rule, not more than one or two. These observations are slightly less distant from the current epoch than the SAI catalogue observations and we were not able to use them prior to building the model of the system published in Volkov et al. [7]. Now we are in a position to correctly classify each of the brightness attenuation and supplement our analysis with these data. We found that the observed brightness attenuations correspond to four secondary and one main minima of the "B" component.
Also, at our request, T. Pribulla and R. Komžik carried out spectroscopic observations of the system using the high-resolution echelle spectrograph mounted on 1.35m reflector at the Skalnaté Pleso observatory, Slovakia Döhring et al. [10]. Obtained spectra contain moving lines of four different stars, with the spectral types close to A0 and approximately of the same brightness. We used the spectra to construct the radial-velocity curves (RVCs) for all four components of the system BU CMi.

INTERSTELLAR EXTINCTION AND TEMPERATURES OF THE COMPONENTS
The star under investigation is at a distance of 250 parsecs from us and is located far enough from the Galactic equator, b = 18 • , so we cannot expect large interstellar absorption. This is confirmed by the position of the star on the two-colour U −B, B −V diagram. We do not show it here, but a similar one can be seen in Fig. 3 of our previous work Volkov et al. [9]. The colour indices of the standard star HD64963 do not show any significant interstellar reddening, too. According to Green et al. [11], zero interstellar absorption in the direction of BU CMi follows up to a distance of several kiloparsecs. The maps of interstellar reddening Schlegel et al. [12], Schlafly and Finkbeiner [13] suggest for BU CMi interstellar extinction E(B − V )=0 m .007 ± 0.0006. Since it is difficult to determine the interstellar extinction with the same accuracy from our two-colour U − B, B − V diagram, we accepted for further analysis this value of extinction and corrected the star colour index from Table II. Despite the fact that reddening is close to zero, its value in the range of the colour indices corresponds to a temperature difference of 140K according to calibration used Popper [14], Flower [15].

THE PHOTOMETRIC ANALYSIS OF THE SYSTEM
In our studies, we solve the complex problem of determining the entire set of interconnected characteristics of a multiple system, and the problem is complicated by the fact that there are two eclipsing binaries and the brightness of all four stars is measured simultaneously. Thus, the number of free parameters is at least double. A number of natural conditions limit the range of possible parameters of the model, which facilitates its construction. These conditions are as follows: First, the combined flux of the system must match the colour indices given in Table II. Second, Kepler's laws must be fulfilled for each of the two eclipsing pairs separately, and for the external orbit for both pairs together. There may be some deviations from these laws, because each of the systems cannot be considered closed, the stars themselves cannot be considered as material points, as evidenced by the apsidal motion and small brightness variations between minima, see Fig. 1. But as a first approximation, the laws should work. Third, the photometrically determined distance to all four stars should be the same. Fourth, the total brightness of all four stars is taken as unity. The last condition immediately implies that we will solve the LCs with a significant contribution of the third light. Since the depths of all four minima are approximately the same and close to 0 m .2, and the power of all lines is also approximately the same, it can be expected that all four stars have a similar surface brightness or, equivalently, temperature. We also assume that the empirical massluminosity law is fulfilled.
The applied technique is described in detail in a number of our previous articles, such as Volkov et al. [9], Volkov et al. [19], Bagaev et al. [20], Kravtsova et al. [21] and Volkov and Kravtsova [22]. Here, we just note that the stellar disk is represented by an ideal circle, the brightness distribution over it is taken into account by concentric circles of different brightness, Khaliullina and Khaliullin [23], in the model of the linear law of darkening to the edge, and the final solution is obtained by the differential corrections method. Edge-darkening coefficients were fixed from the theoretical models Wade and Rucinski [24]. Within the framework of this model, there are no brightness variations between the minima. The analysis revealed that in very accurate TESS observations, see Fig. 1, such changes at the level of fractions of a percent exist, but within the errors of the determined parameters, they do not affect the solutions at the minima.
There are three independent photometric sets in our disposal. The first set is the MASCARA unfiltered observations Burggraaff et al. [3]. We solved them under the assumption that the average wavelength of the observations corresponds to the Johnson R band. The second set comprise our own U BV observations and the third set the TESS observations Stassun et al. [4], to which the Cousins Ic band was attributed. Each of the available observational sets was analyzed separately and the results were compared. We got approximately the same independent parameters, and the light of "A" and "B" components was equally divided in all solutions. This is an important point. At the beginning, we received the same errors for our V observations and TESS points, at the 0 m .006 level. But after a closer look at the residual deviations in TESS observations, it was found that by adding the small corrections to the time of observation at each date, the accuracy of the LC solution can be improved tenfold, from 0 m .006 to 0 m .0006. Later, having plotted on a large scale the graphs of the course of the O − C residuals at the minima timings, we explained the introduction of these corrections by nutation, see the section VI. The real weight of each individual observation increased by a factor of 100. After correcting the TESS observations for this effect, we obtained the most accurate solution, which was taken as the final one. The approximation of our observations and MASCARA observations by this model does not worsen the scatter of points on the corresponding LCs. The results for all LCs are shown in Fig. 2. All these curves are approximated by the same parameters, except for the longitude of the periastron, which changes with time due to the rapid apsidal rotation. The parameters obtained from the solutions are presented in Table III. The masses of the stars were estimated from photometric data using an indirect method, see Volkov et al. [25]. The way in which we estimate the mass errors of the indirect method is described in Volkov et al. [9]. The parameters found from the photometry served as the basis to disentangle the spectroscopic observations. Table III lists the masses already refined from the RVCs.
The obtained relative luminosities and temperatures, in accordance with the calibrations Flower [15], correspond to the non-reddened colour index B − V = −0 m .031. This value is larger by 0 m .004 than that measured in the catalogue Kornilov et al. [26], see Table II. We believe that such a discrepancy is an objective assessment of both the accuracy of the calibrations Flower [15] and the uncertainty of our solution, i.e., the match is satisfactory.

THE SPECTRAL OBSERVATIONS AND THEIR ANALYSIS
Disentangling the spectra in order to measure the radial velocities of the components turned out to be difficult task, because the line profile of each component was almost always in blend with the line profiles of other components. The line profiles of all four components were approximated by Gaussians of the same amplitude and width, and the best fit was sought between the resulting general approximation function and the observed spectrum. We used for spectral type A0 the most powerful optically thin lines MgII (4481Å), see Fig. 3. The results of the analysis are presented in Table V and in Fig. 4, 5. The significant scatter is due to the constant superposition of lines in the spectra. When solving the radial-velocity curves, the apsidal rotation was taken into account. The masses of the components are shown in Table III. These results are not final; work on the spectra is ongoing.

THE APSIDAL ROTATION AND THE LIGHT TIME EFFECT IN MINIMA TIMINGS
To study both effects, it is necessary to know the exact times of the minima. We have determined all possible timings from the MASCARA, TESS and our observations. According to the conditions of the star's visibility above the horizon at the Simeiz observatory of INASAN, we were unable to observe any full eclipse. From the observed individual branches of the minima the average minima were constructed for our observations in 2020. Thus, we have obtained four average eclipse timings, one for each of the stars in the system. Another individual moment of minimum was obtained from our observations at Stará Lesná in 2012. The bibliographic search revealed only one published minimum timing Ogloza et al. [27]. One more timing we derived from the observations of parts of the minima in 2021 using the established geometric model of the star, the four most recent minima timings were obtained. For each of the available timings it was determined which component of the quadruple system is eclipsed, the blending minima were discarded. All obtained data are presented in Tables VI, VII, VIII, IX. The largest contribution to the deviation of the course of the minima timings from the linear formula is made by the apsidal rotation. The amplitude of the effect reaches 0.2 days. The next in importance is the Rømer light-time effect (LITE) with the amplitude of 0.022 days due to the change in the distance to the stars "A" and "B" when they move in a common orbit. Then there are slightly smaller nutation fluctuations with an amplitude of about 0.015 days. The problem was solved using the programs specially developed for this case taking into account our experience in searching for invisible satellites in eclipsing systems Volkov et al. [25], Volkov et al. [28], Volkov and Volkova [29], Volkova and Volkov [30], Volkov [31]. The search for a solution was carried out by the method of successive approximations. Both the parameters of the apsidal rotation and the orbit of the third body were searched simultaneously. In this case, both external orbits for components "A" and "B" had to be identical and differ only in the longitude of the periastron (by 180 degrees) and the amplitudes of the LITE -due to the difference in the masses of both systems. It was the observations of the "A" component that made it possible to fix the parameters of the outer orbit most reliably, since for it, the number of observed minima is maximal and they are distributed over a significant time interval of 34 years. For the apsidal rotation of both eclipsing systems, (asini) A = 3.47±0.08 AU (asini) B = 3.34±0.08 AU e = 0.71 ±0.03, The errors are underestimated, they correspond to the specific model, but there could exist several models which correspond to their own local minima. So the parameters presented will not necessarily be final. The period of the outer orbit is not fully covered by observations yet. Our study of the entire area of feasible solutions showed that there is a set of parameters which differ from each other beyond the errors presented. The real period of the outer orbit surely falls in 5.9-7.7 years interval. We can also estimate that the ratio of the masses of the system "A" to system "B" is approximately 1.0. This result is in good agreement with photometric solutions, in which the luminosities of components "A" and "B" are equally divided. Application of Kepler's third law to the obtained parameters of external orbit gives the total value of the mass of the entire system as M (sin i) 3 = 6.6 M . The solutions to the radial-velocity curves assume M = 13.1 M (the inclination angles of the orbits are known with high accuracy from the photometric solutions). Hence, we conclude that sin i ≈ 0.8, i.e., the angle of inclination of the outer orbit is close to 50 degrees to the line of site.
For the final conclusions about the external orbit of the system further observations are required. We cannot rule out that the above set of parameters of the outer orbit corresponds to one of the local minima of the space of solutions.
A large-scale O − C plots for individual TESS timings is shown in Fig. 10. It can be seen that the primary and the secondary minima for each of the "A" and "B" components slant in antiphase. Such a picture can be obtained if we assume that the elliptical orbits of each of the "A" and "B" components wobble within one degree in the longitude of the periastron. We estimated the periods of such wobbling at 60 days, and for both eclipsing stars they are approximately the same. It turns out that in addition to the usual effects observed in double elliptical systems, such as apsidal motion and the LITE due to the influence of the third body, nutation can also be present.
We do not give here the calculations of the theoretically expected apsidal rotation. For such orbits and masses of the components, it is obviously smaller than the measured values. The reasons for the rapid apsidal rotation must be related to the mutual influence of both systems on each other. Synchronization is clearly observed in the system in the orbital, nutation and apsidal motions of the components "A" and "B": the periods of the eclipsing pairs are P A /P B =0.9, the periods of nutation are approximately equal, the apsidal periods are also the same and U/P 3 ≈ 4. The orbital planes of both eclipsing systems are collinear. It can be concluded that this multiple system has already reached an equilibrium state and is stable.

Conclusions
We have revealed the nature of the unusual quadruple eclipsing system BU CMi, consisting of two eclipsing binaries. The periods of both binaries are tied up by a resonance of 9/10. A fast apsidal rotation is observed in each binary. We built a preliminary orbit in which both binaries rotate around a common center of mass. We found the physical characteristics of all four stars: temperatures, sizes, masses. All stars appear to be very young, aged up to 200 million years. A preliminary comparison of the obtained parameters with theoretical models indicates that the chemical composition of the stars does not correspond to the solar one, and the secondary components are younger than the primary ones. This fact is not surprising in the case of young stars. Under the condition of common origin, less massive stars should set on the ZAMS later than more massive ones. Further intensive photometric and spectral observations are needed to refine the parameters of the mutual orbit and to study nutation. We estimate the angular distance between the "A" and "B" components at the current moment to be about 0.01". So, it is possible to directly measure the components of a binary star using speckle interferometry. BU CMi system is very interesting from the point of view of stellar formation and evolution theory, as well as for studying the dynamics of close multiple systems.
The observations were fulfilled with the 60-cm reflector of the Simeiz observatory of INASAN. This work has used the SIMBAD database of the Strasbourg Center for Astronomical Data (France) and the bibliographic reference service of the ADS (NASA, USA). This document includes data collected by the TESS mission. Funding for the TESS mission is provided by NASA's Science Missions Office. We are grateful to the staff members of the Astronomical Institute of the Slovak Academy of Sciences T. Pribulla and R. Komžik for the spectroscopic observations carried out at our request.