On the long cycle variability of the Algol OGLE-LMC-DPV-065 and its stellar, orbital and disk parameters

We included OGLE-II (Szymanski, 2005)¹¹1http://ogledb.astrouw.edu.pl/$\sim$ogle/photdb/ and OGLE-III/IV data²²2 OGLE-III/IV data kindly provided by the OGLE team.. The OGLE-IV project is described by Udalski, Szymański, & Szymański (2015). Poleski et al. (2010) published the OGLE-II and OGLE-III $I$-band³³3ftp://ftp.astrouw.edu.pl/ogle/ogle3/OIII-CVS/lmc/dpv/phot/I/OGLE-LMC-DPV-065.dat and $V$-band⁴⁴4ftp://ftp.astrouw.edu.pl/ogle/ogle3/OIII-CVS/lmc/dpv/phot/V/OGLE-LMC-DPV-065.dat data of this star. We also considered 460 $B$ magnitudes from the Digitalized Harvard plates (DASCH project)⁵⁵5http://dasch.rc.fas.harvard.edu/project.php covering 59.4 years, since August 1893 to January 1953. In addition, we obtained new photometry with the CTIO 1.3m telescope operated by the SMARTS consortium in service mode between November 2014 and October 2015, with the ANDICAM camera and filters $V$ and $I$. Another data set was collected by Ian Porritt in Turitea Observatory, New Zealand, with the 0.36-meter Meade telescope and a yellow filter. These new data were reduced in the usual way, removing bias and performing flat field corrections in the images and calculating differential magnitudes with respect to constant comparison stars. Finally, 664 $V$-band magnitudes were included from the ASAS-SN catalogue. The photometric observations analyzed in this paper amount to 3099 data points, cover 124 years and are summarized in Table 1.

We separated the light curve into long- and a short-period components. For that we used an algorithm especially designed to disentangle multi-periodic light curves through the analysis of their Fourier component amplitudes. The method is described in Mennickent et al. (2012) and a short summary is given here. The main frequency $f_1$ is found with a period searching algorithm, this is usually the orbital frequency. A least-square fit is then applied to the light curve with a fitting function consisting of a sum of sine functions representing the main frequency and their additional significant harmonics. Afterwards the residuals are inspected for a new periodicity $f_2$. This new periodicity (the long cycle in the case of DPVs) and their harmonics are included in the new fitting procedure. Finally, we obtain the light curve represented by a sum of Fourier components of frequency $f_1$ and $f_2$ and their harmonics. Data residuals with respect to the second and first theoretical light curves are the photometric series representing the orbital and long cycle, respectively.

We were granted 25 hours for spectroscopic observations of the target with the ESO Ultraviolet and Visual Echelle Spectrograph UVES at the Kueyen telescope in the Paranal Observatory in service mode. This is a cross-dispersed echelle spectrograph designed to operate with high efficiency from the atmospheric cut-off at 300 nm to the long- wavelength limit of the CCD detectors (about 1100 nm). The light beam is split into two arms, UV-Blue and Visual-Red, within the instrument. The two arms can be operated separately or in parallel with a dichroic beam splitter. The instrument provides accurate calibration of the wavelength scale down to an accuracy of at least 50 m/s.

With the aim of covering both the orbital as well as the long cycle 27 spectra were secured between October 1, 2013 and February 1, 2015 with the dichroic#2 mode in the ranges 3760$-$4985, 5700$-$7520 and 7665$-$9445 Å. The slit widths of 0$\text{aas}@@fstack\prime \prime$9 at the blue channel and 0$\text{aas}@@fstack\prime \prime$8 at the red channels provided a resolving power of 50 000 and 55 000, respectively. The object was observed at typical airmass 1.4 and with 3000 s exposure time per single exposure. A typical S/N of 65 was achieved at 480 nm. A summary of the spectroscopic observations is given in Table 7.

Zero points can be an issue when combining photometry obtained at different sites with different detectors, filters and sky conditions. We have shifted OGLE IV, ASAS-SN, Turitea, DASCH and CTIO magnitudes to the OGLE-II and OGLE-III averages before performing the analysis described in this paper. In Table 1 we provide only the original magnitude average for OGLE-II and OGLE-III data and also for the DASCH $B$-band magnitude. The light curve in the $I$-band is shown in Fig. 1.

We conducted a search for the orbital period using standard methodologies: eclipse times were measured interactively in the light curve with the computer cursor and a straight line fit was performed with the measured (epoch, time) pairs; the resulting slope gave the orbital period and their error. We also used the Period04 program, that calculates errors based on a Monte Carlo technique (Lenz & Breger 2005). The $I$-band residuals were obtained after removing the long-term cycle (see next section). The periodicities found in different datasets are given in Table 2. We can see that the data are consistent with a constant orbital period; we find the following ephemerides for the main eclipse:

Maxima of the long cycle were measured directly from the light curves and compared with the ephemerides for a 240 day test period, as reported in Table 3. The observed minus calculated ($O-C$) diagram, constructed with the observed times of maxima ($HJ{D}_O$) and the predicted ones ($HJ{D}_C$), shows that the long cycle decreased at the beginning of the observations then remained more or less constant during about 14 cycles (Fig. 2). We notice that considering the MACHO data analyzed by Mennickent et al. (2005, HJD: 2448900-2451500, not included in this paper), which is previous to the OGLE data reported here, the long period has decreased from about 350 to 218 days continuously during about 13 years, before entering in a phase of almost constant period, that lasted for slightly more than 10 yr.

At every epoch we defined a local long cycle period $P_l$, subtracting the observed maximum timing from the previous one and dividing by the number of elapsed times, as given in Table 3. After inspection of Fig. 2 we choose six data intervals characterized by a more or less constant long cycle and large number of observations (Table 4). This procedure allowed to apply the disentangling to every data block considering the variability of the long cycle. The resulting disentangled light curves are shown in Fig. 3, they reveal that the long cycle is double-humped and that it shape remains relatively constant. In addition, the orbital light curve shows a small but significant variability (Figs. 4 and 5): (i) on the 5th interval between $HJD$ 2 455 804 and 2 456 405 the system is brighter at quadratures, and produces larger scatter in the long cycle light curve, (ii) on the first interval the main eclipse seems to be shallower, (iii) significant variability is observed during secondary eclipse; the secondary eclipse seems to occur earlier in interval 1 than in interval 2, and (iv) the shape of the eclipses vary minimally during the maximum, the minimum and the secondary maximum of the long cycle, perhaps the egress of the main eclipse around phases 0.1-0.2 is shallower during the low stage. The changes in timing of minima during the secondary eclipse might indicate changes in the photo-center of the eclipsed or eclipsing source, or changes in circumstellar matter or the donor hemisphere facing the gainer. An unseen/undetected body that dynamically affects the photo-center is another possibility.

We did the same exercise with the $V$-band but we had to use a smaller number of intervals due to the smaller number of observations in this band. The intervals are documented in Table 5. The long cycle usually has a smaller amplitude than in the $I$-band and the orbital light curve shows subtle variability. These changes are better visualized in the combined light curve (Fig. 5).

In order to obtain the radial velocities and orbital parameters we used the KOREL code (Hadrava, 1995, 1997) based on the method of Fourier disentangling, yielding directly the orbital parameters together with the decomposed spectra of the multiple stellar system under study. In addition, we also used the code FOTEL (Hadrava, 1990) to estimate the errors of the orbital parameters.

We notice that the system can be classified as SB2, i.e. both stellar components are detected in the spectrum, in particular in helium and hydrogen lines. The detected components correspond to an early B-type (primary or gainer), and an early A-type (secondary or donor). The method of spectra disentangling does not use any template or another information about the laboratory wavelengths of the spectral lines, therefore the systemic velocity is set to zero. For this reason we adopted an average of systemic velocities calculated by Gaussian adjustments for different spectral lines of each component. For the gainer we obtain ${\gamma }_{pri}=275.1\pm 2.3$ km s${}^{-1}$ and for the donor ${\gamma }_{sec}=279.8\pm 2.8$ km s${}^{-1}$. The lines used in this calculation are shown in Table 6.

We notice that S II/III measurements systematically differ from other lines, suggesting a different formation place. For this reason they are not included in the above calculation. The disentangled spectra are shown in Figs. 6 and 7.

We performed the calculation of radial velocity in seven regions of every spectrum. These regions were chosen because they include several narrow, unblended and well identified metallic lines. All our spectra were prepared with a routine written in IRAF6 and the sampling auxiliary code PREKOR (Hadrava, 2004) was used. To diminish the numerical errors of the disentangling we sampled each spectral region in the maximum number of bins allowed by the code, viz. 4096. This results in the average resolution 0.726 km s${}^{-1}$ per bin, which is higher than the original resolution on the spectrograph detector.

Radial velocities obtained with KOREL for the cases of circular orbit are given in Table 7 and their best fit is shown in Fig. 8. The operation of the KOREL code is described in Hadrava (2004).

The radial velocity for each components is given by

where the sum is realized on the orbits that influence the movement of the star. The true anomaly $\upsilon$ is calculated according to

where $E$ is obtained from the solution of Kepler’s equation.

The orbital parameters obtained by disentangling of the seven spectral regions are summarized in Table 8. The Solution I, which we accept for our modeling of photometry, has been obtained using an independent disentangling of each region separately and then calculating mean solutions and standard deviations of each parameter. Solution II is the simultaneous (”multi-region”) disentangling of all the regions together. The errors of the parameters were obtained using the Bayesian estimate, i.e. from the moments of the Bayesian probability distribution (Hadrava, 2016) We have also solved the radial-velocity curve using the FOTEL code with the input radial velocities obtained from the disentangling. The resulting values of parameters were within the error-bars of the Solution I, but their errors were for about one order underestimated, so we skipped this solution. Finally, the multi-region Solution III is to verify that a possible eccentricity of the orbit can be neglected.

Once disentangled the donor spectrum, we compared it with a grid of solar-metallicity synthetic models constructed with SPECTRUM⁶⁶6http://www.appstate.edu/ grayro/spectrum/spectrum.html and search for the synthetic spectrum minimizing residuals. We find the best fit with a stellar spectrum of $T_{eff,2}$ = 9825 $\pm$ 75 K, $v_{rot,2}sini$ = 53 $\pm$ 3 km s${}^{-1}$ and log g = 3.2 $\pm$ 0.2. Comparisons of the donor disentangled spectrum with the best fit model are shown in Fig. 6. Similarly, from the region 4120$-$4199 Å  we obtained a model with T${}_{eff,1}$ = 22000 K and $v_{rot,1}sini$ = 70.6 km s${}^{-1}$ for the gainer.

Part of the phenomenology of DPVs has been associated with the presence of an optically thick disk around the gainer, probably feed by a Roche-lobe filling donor star (e.g. Garcés L. et al., 2018). Consistently, we model the orbital light curve of OGLE-LMC-DPV-065 considering the stellar fluxes of the two stars, the contribution of an accretion disk around the primary and eventually the light contribution of hot/bright spots located in the outer disk rim. The basic elements of the binary system model with a plane-parallel disk and the corresponding light curve synthesis procedure are described by Djurašević (1992a, b, 1996). The code has been successfully applied to several close binaries including the well-studied binary system $\beta$ Lyrae (e.g. Djurašević et al., 2010, 2012; Garrido et al., 2013; Mennickent et al., 2012; Mennickent & Djurašević, 2013).

We assume that the disk is optically and geometrically thick and that its outer edge is approximated by a cylindrical surface. The vertical thickness of the disk can change linearly with radial distance, allowing different disk’s conical shapes: plane- parallel, concave or convex. The geometrical parameters of the disk are its radius ($R_d$), its vertical thickness at the outer edge (d${}_e$) and the vertical thickness at the inner boundary (d${}_c$). The cylindrical edge of the disk is characterized by its temperature, T${}_d$, and the conical surface of the disk by a radial temperature profile inspired in the temperature distribution proposed by Zola (1991):

We further assume that the disk is in physical contact and thermal equilibrium with the gainer, so its inner radius and corresponding temperature are equal to the radius and temperature of the star ($R_h$, T${}_h$). The temperature of the disk at the edge (T${}_d$) and the temperature exponent (a${}_T$), as well as the radii of the gainer ($R_h$) and of the disk ($R_d$) are free parameters, determined by solving the inverse problem (see Section 4.2).

Motivated by previous research on DPVs (Mennickent et al. 2016), our model includes active regions on the edge of the disk. These active regions are usually revealed in Doppler tomography maps of disks in Algol binaries (Richards, 2004). These regions have higher local temperatures than the disk, and produce a non-uniform distribution of radiation. We consider two active regions: a hot spot (h${}_s$) and a bright spot (b${}_s$), characterized by their temperatures T${}_{hs,bs}$, angular dimensions (radius) ${\theta }_{hs,bs}$ and longitudes ${\lambda }_{hs,bs}$ . The longitude $\lambda$ is measured from the line joining the gainer and donor centers in the direction opposite to the orbital motion in the orbital plane. These parameters are also determined by solving the inverse problem. We also consider a possible departure of symmetry of light emerging from the hot spot due, for instance, to the impact of the gas stream in the disk. This deviation is described by the angle ${\theta }_{rad}$ between the line perpendicular to the local disk edge surface and the direction of the hot spot maximum radiation in the orbital plane. The second spot in the model (here named bright spot), simulates the spiral structure of an accretion disk, observed in hydrodynamical calculations (Heemskerk, 1994). The tidal forces exerted by the donor star produce a spiral shock, observed as one or two extended spiral arms in the outer disk regions. This bright spot can also be interpreted as a region where the disk significantly deviates from the circular shape.

Two potential limitations of the code need to be briefly mentioned: the lack of treatment of the donor irradiation by the hot spot, and the lack of inclusion of a possible not eclipsed additional third light, considering that the long-cycle light was already removed with the process of disentangling. However, the very good fit obtained (based on ${\chi }^2$ minimization) suggests that these additional light sources, if present, are much fainter than the stars and the disk/spots. In addition, while the donor irradiation by the hot spot is not included, the much larger effect of the donor irradiation by the gainer is implemented in our code.

The fit to the orbital light curve was performed using the inverse-problem solving method based on the simplex algorithm, and the model for the binary system described in the previous section. The inverse-problem method is the process of finding the set of parameters that will optimally fit the synthetic light curve to the observed data. We used the Nelder-Mead simplex algorithm (Press et al., 1992) with the optimization described by Dennis & Torczon (1991). For details see Djurašević (1992b).

Based on results of the previous sections we fixed the spectroscopic mass ratio to $q$ = 0.203 and the donor temperature to T${}_2$ = 9825 K. In addition, we set the gravity darkening coefficient and the albedo of the gainer and the donor to ${\beta }_{h,c}$ = 0.25 and A${}_{h,c}$ = 1.0, following von Zeipel’s law for radiative shells and complete re-radiation (von Zeipel, 1924). The limb darkening for the components was calculated as described by Djurašević et al. (2010).

We assume that the donor is rotating synchronously, i.e. the non-synchronous rotation coefficient, defined as the ratio between the actual and the Keplerian angular velocity is f${}_c$ = 1.0. This is justified since it is assumed that the donor has filled its Roche lobe (i.e. the filling factor of the donor was set to F${}_c$ = 1.0), then the synchronization is expected as consequence of the tidal forces.

The case for the gainer is different, since the accreted material from the disk is expected to transfer angular momentum increasing the rotational speed of the gainer up to the critical velocity as soon as even a small fraction of the mass has been transferred (Packet, 1981; de Mink, Pols, & Glebbeek, 2007; Deschamps et al., 2013). For this reason we assumed critical rotation for the gainer, and estimated a non-synchronous rotation factor f${}_h$ = 8.9 in the critical rotation regime.

The best fit along with the $O-C$ residuals, individual donor, disk and gainer flux contributions and the view of the optimal model at orbital phases 0.25, 0.50 and 0.75, are shown in Fig. 9. We note that the residuals show no dependence on orbital or long-cycle phases, except a larger random scatter around main eclipse. Parameters are given in Table 9 and the sensitivity of ${\chi }^2$ with some parameters is illustrated in Fig. 10. We find that at quadrature and $I$-band, the gainer contributes 27% more flux than the donor and the disk only 48% of the donor to the total flux.

We find that the system contains a 13.8 M${}_{\odot }$  star and a 2.81 M${}_{\odot }$  star with absolute magnitude $M_{bol}$ = -6.4 and -3.0 respectively, separated by 49.9 R${}_{\odot }$. The stellar temperatures are T${}_h$ = 25460 K and T${}_c$ (fixed) = 9825 K. The best-fitting model contains an optically and geometrically thick disk around the hotter, more massive gainer star. With a radius of $R_d$ $\approx$ 25 R${}_{\odot }$, the disk is 2.8 times larger than the central star ($R_h$ $\approx$ 8.8 R${}_{\odot }$). The disk has a convex shape, with central thickness d${}_c$ $\approx$ 1.6 R${}_{\odot }$ and edge thickness d${}_e$ $\approx$ 3.5 R${}_{\odot }$. The temperature of the disk decreases from T${}_d$ = 9380 K at its edge, to T${}_h$ = 25460 K at the inner radius.

We notice that the hot spot has 20$\text{aas}@@fstack\circ$5 angular radius and covers 12% of the disk outer rim, and it is situated at longitude ${\lambda }_{hs}$ = 312$\text{aas}@@fstack\circ$4, roughly between the components of the system, at the place where the gas stream falls on to the disk (Lubow & Shu, 1975). The temperature of the hot spot is approximately 18 per cent higher than the disk edge temperature, i.e. T${}_{hs}$ = 11068 K. Although including the hot spot region into the model improves the fit, it cannot explain the light curve asymmetry completely. By introducing one additional bright spot, larger than the hot spot and located on the disk edge at ${\lambda }_{bs}$ = 114$\text{aas}@@fstack\circ$9, the fit becomes much better. This bright spot is 26% hotter than the disk at its edge, i.e. T${}_{bs}$ = 11819 K and has an angular radius of 27$\text{aas}@@fstack\circ$9, covering 16% of the disk outer rim.