Simulation of fluctuating interference pattern for light scattered on erythrocytes in suspension

Orbital period modulation of the eclipsing binary system ER Vulpeculae: Real or not? Alexandru Pop a,⇑ , Ca˘lin Vamos ß b a Romanian Academy, Astronomical Institute, Astronomical Observatory Cluj-Napoca, Str. Cires ßilor 19, RO-400487 Cluj-Napoca, Romania b Romanian Academy, ‘‘T. Popoviciu’’ Institute of Numerical Analysis, P.O. Box 68, 400110 Cluj-Napoca, Romania highlights " We detected a periodicity in the O–C curve of ER Vul and proved its detectability. " We investigated the ‘‘standard’’ mechanisms for orbital period modulation. " We emphasised a statistically significant increase of the orbital period of ER Vul. " We proposed a methodology for estimating the detectability of periodic signals. article info Article history: Received 28 November 2012 Received in revised form 7 February 2013 Accepted 16 February 2013 Available online 26 February 2013 Communicated by E.P.J van den Heuvel Keywords: Methods:data analysis Methods:statistical Binaries:Eclipsing Stars:individual:ER Vul abstract The reanalysis of the timing data of the binary system ER Vulpeculae revealed the superposition of a rel- atively slow increasing linear trend in the orbital period, having a relative period change rate of 2:73  10 8 yr 1 , and a low amplitude (3.5 min) periodic modulation with a period of about 17.78 yr. This periodicity is significantly shorter than that previously published, i.e., 30.6 yr. We propose a meth- odology relying on Monte Carlo simulations which is intended to evaluate the detectability of the inferred signal. Both frequentist and Bayesian approaches have been considered. The actual amplitude of the orbi- tal period modulation might be lower than that obtained via nonlinear least-squares modelling of the timing data, namely about 2.0 min. Consequently, the amplitude of the underlying mechanism will be lower too. In the case of ER Vul, within the hypothesis of a third companion in the system, it could be a low mass star of about 113.8 M Jupiter , but also a brown dwarf of about 64.7 M Jupiter . Further theoretical and observational studies are needed in order to discriminate between this mechanism and the cyclic magnetic activity occurring in both component stars, as possible causes of the observed period variation. Ó 2013 Elsevier B.V. All rights reserved. 1. Introduction The eclipsing binary ER Vulpeculae (HIP 103833, HD 20039, BD+27°3952) is a detached, short-period, partially eclipsing system (P  0.6981 d) consisting of two dwarf stars with spectral types G0 V + G5 V, according to Northcott and Bakos (1967), or G1– 2 V + G3 V, according to Duemmler et al. (2003). It has been classi- fied by Hall (1976) as a short-period RS CVn system (Ekmekçi et al., 2002; Çakirli et al., 2003; Duemmler et al., 2003; Kjurchieva et al., 2003; Shkolnik et al., 2005). Duemmler et al. (2003) argued that the system is still detached, although its primary component al- most fills its Roche lobe. The spectroscopic studies performed by different authors evidenced the presence of circumstellar matter possibly resulted from mass loss phenomena from the stellar com- ponents (e.g., Kjurkchieva et al., 2003), or that of a broad, high- velocity stream of hot gas flowing onto the secondary component of the system (Shkolnik et al., 2005; see also Harmanec et al., 2004). Shkolnik et al. (2005) mentioned that the components of this system ‘‘are the most chromospherically active solar-type stars known’’ (Gunn and Doyle, 1997; Çakirli et al., 2003; Kjurchieva et al., 2003; see also references therein). Duemmler et al. (2003) confirmed the previous result of Gunn and Doyle (1997), that the secondary star is the most active component in the system. The light curve of ER Vul displays distortions and also variability phe- nomena at short-time scales (Mennella, 1990; _ Ibanog˘lu et al., 1993; Oláh et al., 1994; Kjurchieva et al., 2003; Harmanec et al., 2004). Dryomova et al. (2005) included ER Vul in their catalogue of DW-type eclipsing variable stars, which are pre-contact binary systems of short-periodic RS CVn-type evolving towards the W UMa type contact systems. In the SIMBAD database ER Vul is clas- sified as an eclipsing binary of W UMa type. Both the observational features of the DW-type binary systems (e.g., Dryomova et al., 1384-1076/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.newast.2013.02.003 ⇑ Corresponding author. Tel.: +40 264594592; fax: +40 264596889. E-mail addresses: andi_pop@yahoo.com (A. Pop), cvamos@ictp.acad.ro (C. Vamosß). New Astronomy 23-24 (2013) 27–35 Contents lists available at SciVerse ScienceDirect New Astronomy journal homepage: www.elsevier.com/locate/newast

2005), and intricate photometric and spectroscopic behaviour of ER Vul related to the intense chromospheric activity occurring in its component stars, make this system an interesting target for the investigation of the time behaviour of its orbital period. The results of the available studies on the orbital period varia- tion of ER Vul are relatively scarce and rather contradictory. Thus, Milano et al. (1986) and Mennella (1990) established parabolic ephemerides involving a period decrease and a period increase, respectively. Later, Qian et al. (1998) applied the method proposed by Kalimeris et al. (1994) with a fifth order polynomial model for the O–C curve. They inferred that the orbital period of ER Vul dis- play a small amplitude modulation (0.32  10 5 d  0.28 s) with a periodicity of 30.6 yr superposed on a decreasing secular trend. Qian et al. (1998) and Qian (2001) assumed that the periodic mod- ulation may be caused either by the cyclic magnetic activity in one or both components of the system, or by the presence of a very low-mass companion. They explained the secular period decrease of the orbital period of ER Vul through angular momentum loss via magnetic braking. Duemmler et al. (2003) argued that a para- bolic model, corresponding to a period decrease, does not improve significantly the run of the O–C residuals with respect to a linear ephemeris. Moreover, they concluded that the orbital period of ER Vul did not display evidences for a change ‘‘over at least 60 yr’’ (see also Wilson and Raichur, 2011). Harmanec et al. (2004) established a new, improved linear ephemeris. The hypothesis of the constancy of the orbital period was adopted by them after per- forming some numerical tests. However, they did not exclude the possible presence of some ‘‘mild variations of the orbital period’’. Lastly, Liao and Qian (2010) considered that the available data on this system are inadequate for judgement of the character of its orbital period variation. Nanouris et al. (2011) suggested that in case of ER Vul the orbital period evolution could be mainly driven by the mass loss through stellar wind. The main purpose of the present approach is to explore the presence of a low amplitude orbital period modulation in the bin- ary system ER Vul, possibly embedded in the observational noise. Such a possibility is particularly interesting, having in view the spectral types of the two components and their intense chromo- spheric activity (Hall, 1989, 1990). Moreover, as it is well-known, Applegate’s (1992) mechanism explains the orbital period modula- tion phenomena in active binary systems in which only the sec- ondary component is chromospherically active (see also Lanza et al., 1998). Therefore, in case of ER Vul one could expect the coex- istence of two orbital period modulation phenomena. On the other hand, according to the regression relationships obtained by us (Pop and Vamos ß, 2012) reconsidering the correlation between the length of the activity cycle and the angular velocity emphasised by Lanza and Rodonò (1999) for active stars, we found for ER Vul a modulator period either of 27.9 yr or 26.6 yr. These values are obviously in good agreement with the periodicity of 30.6 yr found by Qian et al. (1998), although this periodicity has been proposed without a statistical substantiation, and its estimated amplitude is very low. Because of the above mentioned contradictory results concern- ing the behaviour of the orbital period of ER Vul, and having in view the interest in clarifying the presence of a periodic modula- tion, we proposed in this study a methodology relying both on periodicity detection and O–C curve modelling, and on the evalua- tion of the detectability of the inferred periodic signal. The latter problem has been approached through decision-making, from both frequentist and Bayesian viewpoints (e.g., Appourchaux, 2011). In this context, we introduced the notion of signal detection probabil- ity, by extending the notion of trend detection probability previ- ously defined by Pop et al. (2010). The final goal is to compare the effects of the null hypotheses H 0 (noise) and of the alternative hypothesis H 1 (periodic signal + noise) and to establish if their sig- natures would be perceived in the synthetic data through ampli- tude spectrum analysis. Finally, the astrophysical implications of the obtained results are discussed. 2. Observational data The present study of the orbital period variation of the eclipsing binary system ER Vulpeculae relies only on photoelectric and CCD times of minimum light. We collected 92 primary minima and 78 secondary minima from the Lichtenknecker-Database of the BAV (http://www.bav-astro.de/LkDB) and ‘‘O–C Gateway. Database of times of minima (E) and maxima (RR)’’ of Variable Star Section of Czech Astronomical Society (http://var.astro.cz/ocgate/index.php), as well as the paper of Harmanec et al. (2004), and that of Samolyk (2012). The above data sets have been obtained after all timing data have been checked according to their original sources. During preliminary stages of data analysis we decided to reject several times of minimum light which were found to be outliers according to a kr criterion, where k > 3. Typical values of k used by us were: 5.0, 4.0, 3.5. This test has been applied in an iterative way. The Grubbs test was also used (see Pop and Ba˘rbosu, 2006). The follow- ing primary minima have been rejected: 2435050.552 (Northcott and Bakos, 1956), 2444836.4751, 2444836.4758 (Pohl et al., 1982), 2445993.1749 (Srivastava et al., 1991), 2446265.46002 (Mennella, 1990), 2449563.24449 (Zeinali et al., 1995). We also re- jected the U, B and V times of secondary minimum light at HJD 2449564 (Zeinali et al., 1995). The observed minima corresponding to the same event (eclipse) have been averaged. Thus, we got 48 primary minima and 44 secondary minima, which led us to a final data set consisting of 91 times of minimum light corresponding to distinct events (eclipses). They cover a time base of about 54.9 yr. The preliminary linear ephemeris has been obtained using the method proposed by Pop (1999) (see also Pop et al., 2003), fol- lowed by an improvement through linear regression t n ¼ 2445221:10924ð34Þþ 0:698094903ð39Þ n. 3. Data analysis methods 3.1. O–C curve modelling In the present study on the orbital period variability of ER Vul, we applied the same methodology used in the previous papers of Pop et al. (2011) and Pop and Vamos ß (2012) (see also the refer- ences therein). In the case of ER Vul, the timing data modelling through the classical O–C method relied on the following ephem- eris (e.g., Pop et al., 2003) t n ¼ t 0 þ s 1 n þ s 2 n 2 þ s 11 sinð2pf 01 n þ U 11 Þ; j ¼ 1; 2; ... ; N d ; ð3:1Þ in which the dimensionless frequency f 01 = P p /P s1 (P p  s 1 being the orbital period, and P s1 the period of a possible periodic modulation of the orbital period) has been previously detected, while N d is the number of analysed timing data. The periodicity detection was per- formed on the O–C residuals resulted after removing a parabolic secular trend, using the Monte Carlo type method previously re- viewed by Pop et al. (2010) (see also Pop, 2005, 2007; Pop and Vamos ß, 2007). Note that for real data, the use of the Monte Carlo method for periodicity detection relying on s max requires the preliminary data standardization. In fact, this procedure was already applied in case of Lomb–Scargle periodogram, in which the input data are trans- formed by removing their mean value, and then normalising by their standard deviation (e.g., Press and Rybicki, 1989). We also applied – as a complementary method for time series analysis – the self-correlation method proposed by Percy et al. 28 A. Pop, C. Vamos ß / New Astronomy 23-24 (2013) 27–35

(1981) (see also Percy et al., 1993; Percy et al., 2003). This method, which is conceptually different from that of amplitude spectrum analysis, is particularly useful because the influence of the magnetic activity on the orbital motion may be featured by the lack of strict periodicity (e.g., Applegate, 1992; Lanza et al., 1998). 3.2. Evaluation of the orbital period modulation detectability It is well-known that in case of unevenly sampled data contain- ing a periodic signal only, its amplitude spectrum will display a peak at the signal frequency. Because of spectral leakage (e.g., Deeming, 1975; Scargle, 1982) other peaks will occur at linear combinations between the signal frequency and the frequencies of the peaks appearing in the spectral window of the data. It is also well-known that in case of noisy data (e.g., consisting in pure Gaussian noise) the amplitude spectrum is also noisy (Scargle, 1982). Moreover, numerical simulations relying on Gaussian noise and the observed uneven data sampling revealed that the highest peak in the amplitude spectrum may occur at any frequency in the considered range, but with different probabilities [see Fig. 5 in Pop (2007), and Figs. 7 and 8 in Pop et al. (2010)]. In real situa- tions, especially when we are dealing with a low amplitude peri- odic signal embedded in an additive noise with a relatively high level, our ability to detect such a signal through amplitude power spectrum analysis is determined by the interplay between the sig- nal, the noise (its nature and level), and the spectral window of the available data. Sometimes, we also have to take into account an additional contribution, namely the unequal precision of individual observations. Let us consider that we have detected the presence of a periodic signal (orbital period modulation phenomenon) in the analysed data, and that we have estimated the values of its parameters (fre- quency, amplitude, and phase) together with their corresponding standard errors: f 01 ± r(f 01 ), s 11 ± r(s 11 ), and U 11 ± r(U 11 ). Let r res be the standard deviation of the residuals obtained after removing the periodic component (together with a possible secular trend). In order to estimate the detectability of the above periodic signal, we tested the null hypothesis H 0 (there is no signal in the data), and the alternate hypothesis H 1 (the data contain a well-defined peri- odic signal) by Monte Carlo simulations. We generated for each case N ex artificial time series (O–C curves), with the time sampling of the original data, i.e., the cycle numbers of the observed timing data: (i) N ex data sets consisting of Gaussian noise with zero mean and r res standard deviation (H 0 ) ðO  CÞðn j ; H 0 Þ¼ r res Nð0; 1Þ; j ¼ 1; 2; ... ; N d ; ð3:2Þ (i) N ex data sets consisting of a superposition of the above Gaussian noise over the inferred periodic signal (see Eq. (3.1))(H 1 ) ðO  CÞðn j ; H 1 Þ¼ðO  CÞðn j ; H 0 Þþ k s ~ s 11 sinð2p ~ f 01 n j þ ~ U 11 Þ; j ¼ 1; 2; ... ; N d ; ð3:3Þ in which the parameters of the periodic component are given by ~ n ¼ n þ rðnÞNð0; 1Þ, where n stands for s 11 , f 01 , and U 11 , while k s is a coefficient which allow us to control the signal’s amplitude. Its role is to take into account the effect mentioned by Welsh (1997): superposing a Gaussian noise with different standard deviations over a sinusoidal signal with given amplitude and frequency, the amplitude estimated via least-squares method, ncreases as the amount of noise increases [see also the clean gain used by Roberts et al. (1987, p. 972) in the CLEAN Algorithm]. As it is expected, we found the same result for amplitudes estimated through amplitude spectrum computation. Our approach to evaluate a given periodic signal in a well- defined observational context relies on the investigation, through Monte Carlo simulations, of two types of random processes occur- ring in the amplitude spectrum of synthetic data, in a given fre- quency range: (i) the occurrence of a particular value of a parameter which features the amplitude spectrum, and (ii) the fre- quency at which occurs the peak of a previously chosen rank in the amplitude spectrum. Hence, for each synthetic data set generated under the H 0 and H 1 hypotheses we computed the amplitude spec- trum in the frequency range of interest. Finally, the amplitude spectra for H 0 and H 1 have been separately averaged. For each spectrum we estimated the values of the amplitude of the highest peak (s max ), and the frequency of this peak in the amplitude spectrum. Finally, we got the following information under both the H 0 and H 1 hypotheses: the distribution of s max , and the number of occur- rences of the highest peak at a given frequency within the consid- ered frequency range (N pj , with j = 0 for H 0 , and j = 1 for H 1 ). If f  01 is the frequency at which the highest peak in the amplitude spectrum occurred most frequently, then let us call it the dominant frequency. 3.2.1. The frequentist approach 3.2.1.1. The decision-making relying on the s max distribution. Being given the distribution of s max established under the H 0 and H 1 hypotheses, and its observed value s max;obs , let us denote: N a – the number of s max values under H 0 for which s max > s max;obs , and N b – the number of s max values under H 1 for which s max < s max;obs . The probabilities of committing a type I error (false alarm, false positive) and a type II error (false negative) are respectively a ¼ N a =N ex ; b ¼ N b =N ex : ð3:4Þ Obviously, the ideal situation occurs when a; b ! 0. 3.2.1.2. The decision-making relying on the dominant frequency distribution. In this section we extended the notion of trend detec- tion probability previously introduced by Pop et al. (2010) to the case of periodic signals. Within the frame of the ‘‘classical’’ ap- proach, we are interested in the distribution of the amplitudes of the first order (highest) peak in the amplitude spectra of the artifi- cial data sets generated under the H 0 and H 1 hypotheses. The idea of our approach is to estimate how often the first order peak in the amplitude spectrum (i.e., the signature of the periodic signal) oc- curs in the close proximity to the frequency of interest both under the H 0 and H 1 hypotheses. Our approach is justified in a natural way: the condition to be fulfilled in order to claim the detection of a periodic signal is the occurrence of the first order peak in the amplitude spectrum of the respective data approximately at the signal frequency. Thus, unlike the ‘‘classical’’ strategy which fo- cuses on the amplitude information only, we take into account the amplitude information only by considering the peak order (the first order peak), but we also take into account the information con- cerning the position of the respective peak quantified by the corresponding frequency. From this view point the distribution value of the signal amplitude acquires a secondary importance for estimation of the signal detectability. Let us assume that the highest peak in each of the N ex amplitude spectra occurs most frequently at a frequency value f  01 (the domi- nant frequency) situated in the immediate neighbourhood of the ‘‘observed’’/detected frequency f 01 . Let us also consider a frequency interval Df 0 centred on f  01 :[f  01  Df 0 =2; f  01 þ Df 0 =2]. We define the signal detection probability under the hypothesis H j (j ¼ 0; 1), by means of the following formulae P SD ðf  01 ; Df 0 ; H 0 Þ¼ P SD ðH 0 Þ¼ N p0 =N ex ; P SD ðf  01 ; Df 0 ; H 1 Þ¼ P SD ðH 1 Þ¼ N p1 =N ex : ð3:5Þ A. Pop, C. Vamos ß / New Astronomy 23-24 (2013) 27–35 29

Concerning the choice of the frequency interval Df 0 , we may con- sider (as a maximum value) the frequency resolution corresponding to the available time base and time sampling, i.e., the width of the main peak (at zero dimensionless frequency) of the spectral win- dow function (Df 0w ) (e.g., Deeming, 1975; Roberts et al., 1987). According to the classification scheme given by Zwintz et al. (2000) and adapted by Pop et al. (2010) to the case of a linear trend, we may consider the following thresholds in order to give an interpretation of the results corresponding to the inferred peri- odic model and to a well defined observational context: P SD (H 1 ) > 0.999 – detectable signal, 0:850 < P SD ðH 1 Þ 6 0:999 – possible detectable signal, P SD ðH 0 Þ < P SD ðH 1 Þ 6 0:850 – signal difficult to be detected, P SD ðH 1 Þ 6 P SD ðH 0 Þ – undetectable signal. In case of the null hypothesis H 0 we may distinguish the following thresholds: P SD (H 0 ) > 0.999 – spurious peak, false detection, 0:850 < P SD ðH 0 Þ 6 0:999 – possible spurious peak, possible false detection, P SD ðH 0 Þ 6 0:850 – noisy peak. 3.2.2. The Bayesian approach 3.2.2.1. The decision-making relying on the s max distribution. For s max values which satisfy the condition s max 6 s max;obs , we adopt the decision D 0 , that the hypothesis H 0 is true, while in case of s max val- ues satisfying the condition s max > s max;obs , we adopt the decision D 1 , that the hypothesis H 1 is true. In these conditions, using the Bayes theorem we can define the following four posterior probabil- ities for the two considered hypotheses (e.g., Liteanu and Rîca˘, 1985; Axelsson, 1999; Fawcett, 2006) PðH i jD j Þ¼ PðH i ÞPðD j jH i Þ PðH 0 ÞPðD j jH 0 Þþ PðH 1 ÞPðD j jH 1 Þ ; i; j ¼ 0; 1; ð3:6Þ where P(D j |H i ) are the likelihoods of making the decision D j , under the hypothesis H i PðD 0 jH 0 Þ¼ 1  a; PðD 0 jH 1 Þ¼ b; PðD 1 jH 0 Þ¼ a; PðD 1 jH 1 Þ¼ 1  b; ð3:7Þ and P(H 0 )= P(H 1 ) = 1/2 are the prior probabilities of the two consid- ered hypotheses, about which we have no prior information about their relative importance (e.g., Sturrock and Scargle, 2009). Having in view our goal to evaluate the detectability of the inferred signal, we derived the following posterior probabilities, which can be ex- pressed, using Eqs. (3.6) and (3.7), in terms of the probabilities to commit type I and type II errors, a and b, respectively PðH 1 jD 1 Þ¼ PPV ¼ 1  b 1 þ a  b ; PðH 0 jD 0 Þ¼ NPV ¼ 1  a 1  a þ b ; ð3:8Þ where PPV denotes the positive predictive value (also called the Bayesian detection rate), which is the probability to obtain a s max value, under the hypothesis H 1 , which allow us to make the decision D 1 , and NPV denotes the negative predictive value, which is the prob- ability to obtain a s max value, under the hypothesis H 0 , which allow us to make the decision D 0 . Obviously, a given parameter will have a higher performance in detectability evaluation, if both PPV and NPV values will be higher, tending ideally toward 100%. 3.2.2.2. The decision-making relying on the dominant frequency distribution. We also used the Bayes theorem in the form given by Eq. (3.6), in which the involved direct probabilities are ex- pressed by the signal detection probabilities (see Eq. (3.5)) specific to the frequentist approach PðD 0 jH 0 Þ¼ 1  P SD ðH 0 Þ; PðD 0 jH 1 Þ¼ 1  P SD ðH 1 Þ; PðD 1 jH 0 Þ¼ P SD ðH 0 Þ; PðD 1 jH 1 Þ¼ P SD ðH 1 Þ: ð3:9Þ In these formulae, the data D j (j ¼ 0; 1) refers to the occurrences of the first order peak in the following spectral ranges: D 0 – the fre- quencies outside the considered interval [f  01  Df 0 =2; f  01 þ Df 0 =2], D 1 – the frequencies inside the same interval. Because one has no prior information concerning the relative importance of H 0 and H 1 , we considered P(H 0 )= P(H 1 ) = 1/2. As in the previous case (Sec- tion 3.2.2.1) we will evaluate the detectability of the inferred signal using the previously defined posterior probabilities (3.8), taking into account the above formulae for the direct probabilities (Eq. (3.9)) PðH 1 jD 1 Þ¼ PPV ¼ P SD ðH 1 Þ P SD ðH 0 Þþ P SD ðH 1 Þ ; PðH 0 jD 0 Þ¼ NPV ¼ 1  P SD ðH 0 Þ 2  P SD ðH 0 Þ P SD ðH 1 Þ : ð3:10Þ The above formulae allow us to find which of the two considered hypotheses is more likely to be true in the well-defined circum- stances of our study. At the same time, we get quantitative infor- mation concerning our chance both to detect the inferred periodic signal when it is really present in the observed data, and to get a non-detection result when the inferred signal lacks from these data. Finally, we emphasise that both posterior probabilities P(H 1 | D 1 )= PPV defined by Eqs. (3.8) and (3.10), quantify our chance to detect the inferred signal, being given the sampling pathology of the observed time series (the O–C curve in our case), and the esti- mated noise level, assuming it to be of Gaussian character. 4. Results 4.1. O–C curve analysis and modelling The O–C residuals of ER Vul proved to display a parabolic secu- lar trend, the corresponding ephemeris being t n ¼ 2445221:10773ð46Þþ 0:698094938ð36Þn þ 1:95ð44Þ  10 11 n 2 : ð4:1Þ The relative period change rate is (1/P p )(dP p /dt) = 2.93(65)  10 8 yr 1 , while the residuals standard deviation is r res = 0.00291 d. The analysis of the O–C residuals obtained after removing the parabolic trend reveal the possible presence of a low amplitude (3.5 min) orbital period modulation with a periodicity of 17.78 yr (see Figs. 1 and 2, and Table 1). The statistical significance of the corresponding peak has been estimated using the Monte Carlo type method presented in Section 3.1 performing 100,000 numerical experiments. It supplied the following confidence levels for rejection null hypotheses built up using Gaussian noise, boot- strap resampling, and random permutation of the original data (keeping the original time sampling): 96.117%, 96.390%, and 96.361%, respectively. Such values suggest a moderately good argument to claim the detection of a signal, and thus they may be interpreted as a strong suspicion that we are dealing with a real periodicity buried in observational noise. The self-correlation analysis of the above O–C residuals indi- cated the presence of a (quasi)periodicity of about 18.6 yr (Fig. 3). On the other hand, having in view the error bars of the first two deepest minima appearing in the self-correlation diagram (SCD), which occur at 18.5 yr and 37.5 yr (close to 2  18.5 yr), we may consider that, within the limits of the observational errors, 30 A. Pop, C. Vamos ß / New Astronomy 23-24 (2013) 27–35

this pattern belongs to a periodic process. Due to the noisy charac- ter of SCD, we estimated the semi-amplitude of this periodicity by taking into account the D(O–C) values belonging to the interval within the minimum value and the first quartile (for estimating the minimum level), and those situated between the third quartile and the maximum value (for estimating the maximum level). Com- puting the corresponding weighted means, we found the semi- amplitude s SCD ¼ 0:00140  0:00011 d. The results supplied by the two complementary methods are in good agreement (see Sec- tion 4.2 below), and consequently, this fact sustains the hypothesis concerning the reality of the inferred periodicity. The final residuals obtained after removing the model given by Eq. (3.1) with the parameters listed in Table 1 are practically cen- tred on zero (their mean value is 2:8  10 12 d). They have a Gaussian distribution, according to the Kolmogorov–Smirnov test, while the application of the runs test proved their random character. 4.2. The evaluation of the detectability of the detected orbital period modulation We applied the methodology described in Section 3.2 using the values of the parameters of the inferred periodic model given in Ta- ble 1. The number of numerical experiments was 100,000. We per- formed Monte Carlo simulations for amplitude of the periodic signal of the form k s s 11 , with k s < 1.0. The amplitude value esti- mated from the amplitude spectrum of the O–C residuals after removing the parabolic trend (see Fig. 1, upper panel) is s max;obs ¼ 0:00167 d. An equal value of the highest peak in the syn- thetic averaged spectrum is obtained for an ‘‘input’’ amplitude (see Table 1) attenuated with a factor k s = 0.577, i.e., k s s 11 ¼ 0:00142 d. Remark the excellent agreement between the above estimated val- ues of s SCD and k s s 11 . This value of k s has to be considered from a statistical perspective: the estimated k s value corresponds to the most probable value of the Fourier amplitude at the frequency f  01 which is equal to s max;obs (see Fig. 4). But for lower probabilities, for other k s values, one may find good agreements between s max;obs and the highest peak of amplitude spectra of other artificial time series. On the other hand, s max;obs also represents a quantity associated to a realisation of a stochastic process: the amplitude spectrum computed on the basis of the observed timing data, which include the contributions of different noise sources related both to the observed binary system, and the observational tech- nique as well as the data processing details. Obviously, we have no information about the actual amplitude of the involved periodic signal. Consequently, it is natural that our estimates should contain some amount of incertitude. However, this approach is able to sup- ply a more realistic perspective on the interplay of these effects. The ‘‘observed’’ spectrum and the averaged spectra generated under the two hypotheses H 0 and H 1 (Eqs. (3.2) and (3.3)) are dis- played in Fig. 4. The spectrum obtained under H 0 is obviously white. The comparison between the ‘‘observed’’ spectrum and the averaged one reveals the persistence of some spectral features at well-defined frequencies, the most persistent being the peak sit- uated at f  01 . The details of signal detectability evaluation relying on the number of occurrences of the highest (first order) peak in the amplitude spectra of artificial data (see Sections 3.2.1.1 and 3.2.2.2) are presented in Fig. 5. We considered the frequency inter- val ½0:000074; 0:000138 of width Df 0w = 0.000064 (see Sec- tion 3.2.1.2) centred on the dominant frequency f  01 ¼ 0:000106. Fig. 1. Amplitude spectra of the O–C curve of ER Vul obtained after removing the parabolic trend (upper panel), and the parabolic trend + the inferred periodic signal (lower panel). Fig. 2. The parabolic + periodic model of the O–C curve of ER Vul. Table 1 Parameters of the parabolic + periodic ephemeris of the timing data on ER Vul. t 0 ¼ HJD2445221:10744 0:00046 s 1  Pp ¼ 0:698094884 dc 1 0:000000037 s 2 ¼ 1:82  10 11 dc 2 0:52  10 11 ð1=PpÞðdPp=dtÞ¼ 2:73  10 8 yr 1 0:79  10 8 f 01 ¼ 0:0001075 0:0000032 P S1 ¼ 17:78 yr 0:53 s 11 ¼ 0:00246 d 0:00048 U 11 ¼ 6:17 rad 0:19 rres ¼ 0:00257 d A. Pop, C. Vamos ß / New Astronomy 23-24 (2013) 27–35 31

Fig. 5 displays the histograms of the number of occurrences of the highest peak under H 0 and H 1 hypotheses, for k s = 0.577. The corre- sponding values of the signal detection probability and of the pos- terior probabilities are also mentioned. The peak appearing at f  01 under H 1 is from far the most prominent peak, its occurrence being near two magnitude orders higher than that which appeared under H 0 at the same frequency, and about 21 times higher than the high- est peak which appeared under H 0 . The results of the above analyses are summarized in Table 2, in which are listed the values of different parameters featuring the detectability of the inferred periodic signal, corresponding to its observed signature: s max;obs ¼ 0:00167 d (see Fig. 1, upper panel). The results corresponding to the Bayesian approach re- veals two interesting conclusions: (i) there is an excellent agree- ment between the values of the posterior probabilities P(H 1 |D 1 ) evaluated in case of s max (about 98.9%, and 98.1%), for the two values of k s , and (ii) the posterior probabilities P(H 1 |D 1 ) evalu- ated in case of the analysis of the dominant frequency distribu- tion (about 95.6%, and 93.9%) also display a good agreement for the two considered values of k s . The overall conclusion is that the hypothesis H 1 which involves the presence of a periodic modulation of the orbital period of ER Vul may be taken into ac- count as a working hypothesis at least. The low amplitude signal with a period of 17.78 yr proved to be possible to be detected relying on the available timing data. 4.3. The secular trend in the orbital period The available timing data enabled us to emphasise the presence of a relatively slow increasing linear trend in the orbital period of ER Vul, which yields a concave shape of the O–C residuals (Fig. 2). The coefficient of the second order term proved to be sta- tistically significant at a significance level of 0.00078. One remarks the good agreement between the values of the second order coef- ficient, and those of the relative period change rate corresponding to the adopted ephemerides (see Eq. (4.1), and Eq. (3.1) with Ta- ble 1). This result is in agreement with the recent inference of Nanouris et al. (2011) according to which the orbital period evolu- tion of ER Vul is expected to be mainly driven by mass loss through stellar wind (see also Kjurkchieva et al., 2003). As we already men- tioned in Section 1, this binary system is detached, and conse- quently this secular orbital period variation cannot be explained through the mass transfer in the system. 4.4. The third body hypothesis Having in view the above results on periodicity analysis of the O–C residuals of ER Vul, we approached the interpretation of the parameters of the detected orbital period modulation phenomenon in terms of a Keplerian motion related to the presence of an unseen Fig. 3. The self-correlation diagram of the O–C residuals of ER Vul obtained after removing the parabolic trend. Fig. 4. Averaged artificial amplitude spectra for k s = 0.577 together with the ‘‘observed’’ amplitude spectrum (identical with that plotted in Fig. 1, upper panel). Fig. 5. Number of occurrences of the highest peak in the amplitude spectrum in the considered frequency range, for k s = 0.577. Table 2 Values of the different probabilities quantifying the detectability of the inferred orbital period modulation in the O–C curve of ER Vul. Frequentist approach Bayesian approach k s = 1.0 smax : a ¼ 0:01025 b ¼ 0:03879 smax : PðH 0 jD 0 Þ¼ 96:229% PðH 1 jD 1 Þ¼ 98:945% P SD ðf  01 ; Df 0 ; H 0 Þ¼ 4:472% P SD ðf  01 ; Df 0 ; H 1 Þ¼ 97:284% PðH 0 jD 0 Þ¼ 97:235% PðH 1 jD 1 Þ¼ 95:605% k s = 0.577 smax : a ¼ 0:01025 b ¼ 0:46318 smax : PðH 0 jD 0 Þ¼ 68:121% PðH 1 jD 1 Þ¼ 98:126% P SD ðf  01 ; Df 0 ; H 0 Þ¼ 4:472% P SD ðf  01 ; Df 0 ; H 1 Þ¼ 68:368% PðH 0 jD 0 Þ¼ 75:124% PðH 1 jD 1 Þ¼ 93:861% 32 A. Pop, C. Vamos ß / New Astronomy 23-24 (2013) 27–35

companion. Moreover, taking into account the monoperiodic char- acter of the inferred modulation, we assumed that the hypothetical companion revolves on a circular orbit around the barycentre of the system (close binary + companion). For the detected periodic- ity (P s = 17.78 yr), we estimated the values of the corresponding parameters: a 12 sin i – the projection of the semimajor axis of the close binary absolute orbit on the line of sight, i – inclination of the normal to the orbit plane on the line of sight, f(M) – the mass function, M 3 – the mass of the hypothetical companion, and K – the semi-amplitude of the radial velocity curve. Relying on the results of Wilson and Raichur (2011), we estimated the mass of the binary system M 12 ¼ 1:939  0:021 M Sun , with M 1 ¼ 0:9946 0:0108 M Sun , and M 2 ¼ 0:9447  0:0113 M Sun . Thus, for the esti- mated value of s 11 (see Table 1), we obtained: a 12 sin i ¼ 0:426 AU, f ðMÞ¼ 2:44  10 4 M Sun ¼ 0:256 M Jupiter , and K ¼ 0:71kms 1 , while for the corrected value of s 11 , i.e., k s s 11 ¼ 0:00142 d, we got significantly different values of the above parameters: a 12 sin i ¼ 0:246 AU, f ðMÞ¼ 4:7  10 5 M Sun ¼ 0:049 M Jupiter , and K ¼ 0:41kms 1 . We estimated the values of the mass (M 3 ), and the orbit radius (a 3 ) of the hypothetic companion of ER Vul as functions of inclination value, for the above two con- sidered situations (see Fig. 6). Taking into account the inclination estimated by Wilson and Raichur (2011), i = 67.95°, we obtained in the case of s 11 : M 3 ¼ 0:109 M Sun ¼ 113:8 M Jupiter , and a 3 ¼ 8:19 AU, while for k s s 11 : M 3 ¼ 0:062 M Sun ¼ 64:7 M Jupiter , and a 3 ¼ 8:32 AU. It is interesting that in case of the corrected (semi- )amplitude value of the periodic component of the O–C curve (k s s 11 ), and for inclination values higher than about 53.3°, the third body in the system should be a brown dwarf (M 3 < 75 M Jupiter ). According to the distance estimates given by Bilir et al. (2008), and Wilson and Raichur (2011), i.e., d ¼ 48:48  52:0 pc, we found for the minimum apparent separation of the wide system 0.062 arcsec – 0.067 arcsec, while for the maximum apparent sep- aration we found 0.165 arcsec – 0.178 arcsec. We investigated the stability of the two hypothetical triple sys- tems in the frame of the two-body approximation using the simple stability criterion of Szebehely and Zare (1977). The coplanarity of the orbits of the inner close binary system (ER Vul) and wide sys- tem (close binary + companion) is obviously assumed. It resulted that in both cases (s 11 , and k s s 11 ) the respective triple systems would be stable either for direct or retrograde motion, but ‘‘more stable’’ in case of direct motion. 4.5. The magnetic activity cycles hypothesis As we already mentioned in Section 1, both component stars of ER Vul display intense stellar activity and therefore a modulation of the orbital period induced by the cyclic magnetic activity is to be expected to occur (e.g., Hall, 1989, 1990). The mechanism used to explain such a phenomenon, and also to investigate its implica- tions, i.e., to estimate the values of some specific physical parame- ters, is that proposed by Applegate (1992), which was further investigated by Lanza et al. (1998), Lanza and Rodonò (1999), and Lanza (2005). From the beginning we have to remind that this approach relies on a well-defined model for the considered binary system: it consists of a magnetically active star and a companion star which has no magnetic activity, and which can be regarded as a point mass. Consequently, having in view that both compo- nents of ER Vul have late spectral types G1-2 V + G3 V, and are magnetically active (Duemmler et al., 2003), Applegate’s (1992) model is not suitable to be applied in this case. However, for the two amplitudes of the periodic term (s 11 and k s s 11 ), we computed some overall parameters (e.g., Applegate, 1992; Lanza and Rodonò, 1999) for the component stars of the sys- tem: A O–C – the semi-amplitude of the O–C curve, DP p /P p – the rel- ative semi-amplitude of the orbital period modulation, DQ – the variation of the quadrupole moment of the active star which would be related to the observed DP p /P p value. The following values have been obtained for these parameters (the orbital period is P S = 17.78 yr, see Table 1):  A O–C = s 11 = 0.00246 d, DP p /P p = 2.38  10 6 , DQ 1 = 4.31  10 42 kg m 2 , DQ 2 = 4.10  10 42 kg m 2 ,  A OC ¼ 0:577 s 11 ¼ 0:00142 d, DP p /P p = 1.37  10 6 , DQ 1 = 2.49  10 42 kg m 2 , DQ 2 = 2.36  10 42 kg m 2 . The DP p /P p values obtained for the two considered O–C curve (semi-) amplitudes are typical for W UMa systems. This result is in agreement with the evolutionary stage of ER Vul (see Section 1). The DQ values for the two component stars are of the same order of magnitude with those computed for other W UMa systems, e.g., AB And (Borkovits et al., 2005), and CK Boo (Pop and Vamos ß, 2012), or for cataclysmic variables (e.g., Lanza and Rodonò, 1999). We also remark that the periodicity found by us (17.78 yr) sig- nificantly differs from the prediction relying on the correlation Fig. 6. The mass (lower panel) and the semimajor axis of the absolute orbit (upper panel) of the hypothetical companion for different inclination values, and for k s = 1.0, and k s = 0.577. The mass at which occurs the brown dwarf/main sequence cut is marked with a dashed line. A. Pop, C. Vamos ß / New Astronomy 23-24 (2013) 27–35 33

suggested by Lanza and Rodonò (1999) (27.9 or 26.6 yr; see Sec- tion 1). It does not necessarily mean that the orbital modulation mechanism could not involve the cyclic magnetic activity of the component stars of ER Vul, but it is another reason to deepen the study of this correlation. An interesting approach to explain the orbital period modula- tion in magnetically active close binaries has been proposed by Yuan and Qian (2007). In their model, the authors take into account the magnetic connexion of the binary system components, when their separation is very small, i.e., a/(R 1 + R 2 )  1–3, and conse- quently, the stars will magnetise each other. Making some simpli- fying assumptions, Yuan and Qian (2007) succeed to derive analytic formulae for the magnetic field inside the two component stars, and for the mean magnetic field that we can observe. Given the above mentioned characteristics of ER Vul (both components are magnetically active, and a/(R 1 + R 2 ) = 1.9), this model can be ap- plied for this binary system. We computed the following estimates for the above magnetic fields:  A O–C = s 11 : B 0 = 75.1 T, B s ¼ 6:0  9:0 T,  A OC ¼ 0:577 s 11 : B 0 = 57.0 T, B s ¼ 4:6  6:8 T. These values are in relatively good agreement with the estimates given by Yuan and Qian (2007) in their Table 1. It is interesting to remark that the above B s values are also in agreement with the values of mean magnetic field involved in the model of Lanza et al. (1998) (1–6 T). 5. Concluding remarks The orbital period variability of the close binary system ER Vul has been reanalyzed relying on 91 times of primary and secondary minimum light which cover a time base of about 54.9 yr. We found that the run of the O–C residuals can be described by a statistically significant secular parabolic trend corresponding to a relatively slow increasing trend in the orbital period featured by a relative period change rate of 2:73  10 8 yr 1 , and a 17.78 yr periodic component, with a low amplitude of about 3.5 min. Although the confidence level for the rejection of the null hypothesis (a little over 96%) was lower than the threshold usually taken into account (99.9%), we decided to claim the detection of this periodic signal because of its confirmation through the complementary method of self-correlation analysis. However, having in view the rather subliminal character of this modulation, we approached the prob- lem of evaluation the detectability of such signal. We have estimated our chance to detect a signal having the same parameters as those estimated, in data sets consisting in the same amount of data, the same sampling, and the same noise level as that estimated (assuming a Gaussian distribution). We pro- posed two complementary methodologies based on binary deci- sion making relying on Monte Carlo simulations intended for detectability evaluation of periodic signals in the presence of high level additive noise. For each of them we gave both frequentist and Bayesian approaches. The latter supplied the most consistent and also optimistic results. The application of these methodologies supplied us consistent results (see Table 2) and lead us to conclude that the inferred orbi- tal period modulation phenomenon may be real. An additional and important conclusion was that the actual amplitude of the signal is lower (2.0 min) than that estimated by us through the usual data modelling technique. It has direct implications for a realistic phys- ical interpretation of the observed modulation phenomenon. At this moment, in case of ER Vul it is difficult to give a definite an- swer to this problem. Both considered mechanisms, light-time ef- fect due to the presence of an unseen companion (a low mass star or even a relatively high mass brown dwarf), and the magnetic activity cycles of the component stars seem plausible. Obviously, the results of the present study of the orbital period variability of ER Vul need further observational confirmation. New, high quality timing and radial velocity data will be necessary for this reason. Long time base and precise photometric data with high temporal resolution are also required in order to investigate the possible influence of light curve distortions and rapid variability (see Section 1) on the estimation of the times of minima, and thus, their influence on the shape of O–C curve. Furthermore, the devel- opment of a new model to explain orbital period modulation by magnetic cycles occurring in close binary systems whose both components are active (see also Lanza, 2005) will be extremely important for discriminating between different possible mechanisms. Acknowledgments This research has made use of NASA’s Astrophysics Data System Bibliographic Services, SIMBAD database [SIMBAD4 rel. 1.196 – 2012.07.16], operated at CDS (Strasbourg, France), the Lichtenk- encker-Database of the BAV operated by the Bundesdeutsche Arbe- itsgemeinschaft für Veränderliche Sterne e.V. (BAV), and ‘‘O–C Gateway. Database of times of minima (E) and maxima (RR)’’ of Var- iable Star Section of Czech Astronomical Society. References Applegate, J.H., 1992. ApJ 385, 621. Appourchaux, T., 2011. arXiv:1103.5352v2 [astro-ph.IM]. Axelsson, S., 1999. The base-rate fallacy and its implications for the difficulty of intrusion detection. In: Proceedings of the 6th ACM Conference on Computer and Communications Security (CCS ’99). Kent Ridge Digital Labs, Singapore, pp. 1–7. Bilir, S., Ak, T., Soydugan, E., Soydugan, F., Yaz, E., Filiz Ak, N., Eker, Z., Demircan, O., Helavci, M., 2008. AN 329, 835. Borkovits, T., Elkhateeb, M.M., Csizmadia, Sz., Nuspl, J., Bíró, I.B., Hegedüs, T., Csorvási, R., 2005. A&A 441, 1087. Çakirli, Ö., _ Ibanog˘lu, C., Fraska, A., Catalano, S., 2003. A&A 400, 257. Deeming, T.J., 1975. Ap&SS 36, 137. Dryomova, G., Perevozkina, E., Svechnikov, M., 2005. A&A 437, 375. Duemmler, R., Doucet, C., Formanek, F., Ilyin, I., Tuominen, I., 2003. A&A 402, 745. Ekmekçi, F., Özeren, F.F., Ak, H., 2002. AN 323, 31. Fawcett, T., 2006. Pattern Recogn. Lett. 27, 861. Gunn, A.G., Doyle, J.G., 1997. A&A 318, 60. Hall, D.S., 1976. The RS CVn binaries and binaries with similar properties. In: Fitch, W. (Ed.), Multiple Periodic Variable Stars. Reidel, Dordrecht, pp. 287–348. Hall, D.S., 1989. SSR 50, 219. Hall, D.S., 1990. Period changes and magnetic cycles. In: _ Ibanog˘lu, C. (Ed.), Active Close Binaries. Kluwer Academic Publishers, Dordrecht, pp. 95–119. Harmanec, P., Bozˇic´ , H., Thanjavour, K., Robb, R.M., Ruzˇdjak, D., Sudar, D., 2004. A&A 415, 289. _ Ibanog˘lu, C., Evren, S., Akan, M.C., Tunca, Z., Keskin, V., 1993. A&A 269, 310. Kalimeris, A., Rovithis-Livaniou, H., Rovithis, P., 1994. A&A 282, 775. Kjurchieva, D.P., Marchev, D.V., Zola, S., 2003. A&A 404, 611. Lanza, A.F., 2005. MNRAS 364, 238. Lanza, A.F., Rodonò, M., 1999. A&A 349, 887. Lanza, A.F., Rodonò, M., Rosner, R., 1998. MNRAS 296, 893. Liao, W.-P., Qian, S.-B., 2010. MNRAS 405, 1930. Liteanu, C., Rîca˘, I., 1985. Optimizarea proceselor analitice. Editura Academiei Republicii Socialiste România, Bucuresßti. Mennella, V., 1990. A&A 234, 203. Milano, L., Mancuso, S., Vittone, A., D’Orsi, A., Marcozzi, S., 1986. Ap&SS 124, 83. Nanouris, N., Kalimeris, A., Antonopoulou, E., Rovithis-Livaniou, H., 2011. A&A 535, 1–11. Northcott, R.J., Bakos, G.A., 1967. AJ 72, 89. Oláh, K., Budding, E., Kim, H.-I., Etzel, P., 1994. A&A 291, 110. Percy, J.R., Hosick, J., Leigh, N.W.C., 2003. PASP 115, 59. Percy, J.R., Jakate, S.M., Matthews, J.M., 1981. AJ 86, 53. Percy, J.R., Ralli, J.A., Sen, L.V., 1993. PASP 105, 287. Pohl, E., Evren, S., Tümer, O., Sezer, C., 1982. IBVS, No. 2189. Pop, A., 1999. IBVS, No. 4801. Pop, A., 2005. On the orbital period modulation of the eclipsing binary system Y Leonis. In: Sterken, C. (Ed.), The Light-Time Effect in Astrophysics. ASP Conf. Ser., pp. 263–268, vol. 335. Pop, A., 2007. Romanian Astron. J. 17, 35. Pop, A., Ba˘rbosu, M., 2006. Romanian Astron. J. 16, 179–183. 34 A. Pop, C. Vamos ß / New Astronomy 23-24 (2013) 27–35

Pop, A., Vamosß, C., 2007. Romanian Astron. J. 17, 101. Pop, A., Vamosß, C., 2012. NewA 17, 667. Pop, A., Liteanu, V., Moldovan, D., 2003. Ap&SS 284, 1207. Pop, A., Turcu, V., Marcu, A., 2011. Ap&SS 333, 17. Pop, A., Vamosß, C., Turcu, V., 2010. AJ 139, 425. Press, W.H., Rybicki, G.B., 1989. ApJ 338, 277. Qian, S.B., 2001. Possible magnetic activity cycles of four chromospherically active binaries: ER Vul, UV Psc, AR Lac and BH Vir. In: Brekke, P., Fleck, B., Gurman, J.B. (Eds.), Proc. IAU Symposium, 203, Recent Insights into the Physics of the Sun and Heliosphere: Highlights from SOHO and Other Space Missions, pp. 437– 440. Qian, S., Liu, Q., Yang, Y., 1998. Ap&SS 257, 1. Roberts, D.H., Lehár, J., Dreher, J.W., 1987. AJ 93, 968. Samolyk, G., 2012. JAAVSO 40, 975. Scargle, J.D., 1982. ApJ 263, 835. Shkolnik, E., Walker, G.A.H., Rucinski, S.M., Bohlender, D.A., Davidge, T.J., 2005. AJ 130, 799. Srivastava, R.K., Padalia, T.D., Srivastava, J.B., 1991. Ap&SS 182, 281. Sturrock, P.A., Scargle, J.D., 2009. ApJ 706, 393. Szebehely, V., Zare, K., 1977. A&A 58, 145. Welsh, W.F., 1997. Measuring variability in the presence of noise. In: Maoz, D., Sternberg, A., Leibowitz, E.M. (Eds.), Astronomical Time Series, Astrophysics and Space Science Library. Kluwer Academic Publishers, Dordrecht, pp. 171–174, vol. 218. Wilson, R.E., Raichur, H., 2011. MNRAS 415, 596. Yuan, J., Qian, S., 2007. ApJ 669, L93. Zeinali, F., Edalati, M.T., Mirtorabi, M.T., 1995. IBVS, No. 4190. Zwintz, K., Weiss, W.W., Kuschnig, R., Gruber, R., Frandsen, S., Gray, R., Jenkner, H., 2000. A&AS 145, 481. A. Pop, C. Vamos ß / New Astronomy 23-24 (2013) 27–35 35