key: cord-0264663-4e3ntoy8
authors: Jetsu, L.
title: Say hello to Algol's new companion candidates
date: 2020-05-27
journal: nan
DOI: 10.3847/1538-4357/ac1351
sha: be78512fdb9ef0d8c2cd3498af1de96a437a3ff3
doc_id: 264663
cord_uid: 4e3ntoy8

Constant orbital period ephemerides of eclipsing binaries give the computed eclipse epochs (C). These ephemerides based on the old data can not accurately predict the observed future eclipse epochs (O). Predictability can be improved by removing linear or quadratic trends from the O-C data. Additional companions in an eclipsing binary system cause light-time travel effects that are observed as strictly periodic O-C changes. Recently, Hajdu et al. (2019) estimated that the probability for detecting the periods of two new companions from the O-C data is only 0.00005. We apply the new Discrete Chi-square Method (DCM) to 236 years of O-C data of the eclipsing binary Algol ($beta$ Persei). We detect the tentative signals of at least five companion candidates having periods between 1.863 and 219.0 years. The weakest one of these five signals does not reveal a"new"companion candidate, because its $680.4 pm 0.4$ days signal period differs only $1.4 sigma$ from the well-known $679.85 pm 0.04$ days orbital period of Algol~C. We detect these same signals also from the first 226.2 years of data, and they give an excellent prediction for the last 9.2 years of our data. The orbital planes of Algol~C and the new companion candidates are probably co-planar, because no changes have been observed in Algol's eclipses. The 2.867 days orbital period has been constant since it was determined by Goodricke (1783).

The oldest preserved historical document of the discovery of a variable star is the Ancient Egyptian papyrus Cairo 86637, where naked eye observations of Algol's eclipses have been recorded into the Calendar of Lucky and Unlucky days (Porceddu et al. 2008; Jetsu et al. 2013; Jetsu & Porceddu 2015; Porceddu et al. 2018) . Montanari re-discovered its variability in the year 1669. Goodricke (1783) determined the orbital period P orb = 2. d 867 of this eclipsing binary (EB). The close orbit eclipsing stars are Algol A (B8 V) and Algol B (K2 IV). Curtiss (1908) discovered the 1. y 863 wide orbit third companion Algol C (K2 IV). Direct interferomet-ric images of these three members have been obtained (e.g. Zavala et al. 2010 ; Baron et al. 2012) .

Periodic long-term changes occur in the observed (O) minus the computed (C) primary eclipse epochs of EBs. The most probable causes are a third body (e.g. Li et al. 2018 ), a magnetic activity cycle (e.g. Applegate 1992) or an apsidal motion (e.g. Borkovits et al. 2005) . Hajdu et al. (2019) searched for third bodies in a large sample of 80 000 EBs. They detected 992 triple systems from the O-C data, and only four candidates that may have a fourth body. Their fourth body detection rate was 4/80 000 = 0.00005. Recently, Jetsu (2020, hereafter Paper I) formulated the new Discrete Chi-Square Method (DCM). He applied DCM to the O-C data of XZ And, and detected the periods of a third and a fourth body.

In Algol, the mass transfer from the less massive Algol B (0.8m ) to the more massive Algol A (3.7m ) arXiv:2005.13360v4 [astro-ph.SR] 2 Nov 2021 should cause a long-term P orb period increase (Kwee 1958) , which should have been observed as quadratic long-term O-C changes (Kiseleva et al. 1998 ). Longterm P orb increase or quadratic O-C changes have not been observed in Algol since its period was determined 238 years ago. However, its orbital period modulation does cause negative and positive O-C changes. The short-term low amplitude O-C changes follow 1. y 863 orbital motion cycle of Algol C, while the high amplitude O-C changes follow 30 y and 200 y quasi-periodic activity cycles (Applegate 1992) . The physical origin of period changes is not fully understood, because Algol's puzzling O-C diagram contains unknown signals and trends (e.g. Frieboes-Conde et al. 1970; Applegate 1992) . We apply DCM to Algol's O-C data, because this method can detect many signals superimposed on unknown trends. Kim et al. (2018) note that their TIDAK database O-C ephemerides "cannot be used for the prediction of future times of the primary or secondary minima." These ephemerides are determined by eliminating linear or quadratic trends from the available O-C data (Kreiner et al. 2001) . They usually need to be re-determined when new data are obtained. Although the O-C changes caused by a third body are strictly periodic, the predictions usually fail to separate aperiodic trends from periodic signals (e.g. Bours et al. 2014; Lohr et al. 2015; Song et al. 2019) . Furthermore, the detection rate of third bodies from O-C data is extremely low (e.g. Hajdu et al. 2019) . Against this background, it is totally unexpected that we can detect numerous periods in Algol's O-C data, as well as predict its O-C changes.

The epochs of the observed light curve minima give the observed (O) values. We obtained the n = 2238 observed eclipse epochs of Algol from the 2018 version of TIDAK database (Kim et al. 2018) . These eclipse epochs have been determined by hundreds of astronomers during the past two centuries. The nights when these eclipses could be observed were known beforehand. Every eclipse lasted eight hours. Both dimming and brightening took four hours. The probability for a negative or positive mid eclipse epoch error was the same, because the eclipse light curve was symmetric. It is therefore probable that the observational errors follow a Gaussian distribution, the epoch values contain no observational trends, and the observational errors are not heteroskedastic. Naturally, the accuracy of these data improves towards modern times, because the observational techniques have improved. We study only the primary minimum epochs when the dimmer Algol B eclipses the brighter Algol A. Therefore, we reject all fourteen secondary minima, because they occur P orb /2 after the primary minima. We analyse only the remaining n = 2224 primary minima between November 12th, 1782 and October 18th, 2018. These data are given in Table A4 (∆T = 86171 d = 236 y ). We obtain the computed (C) epoch values from the TIDAK database ephemeris HJD 2445641.5135 + 2.86730431E.

(1)

This ephemeris predicts that all Algol's primary eclipses occur at multiples HJD 2445641.5135 + E × P orb , where P orb = 2. d 86730431 is the orbital period of Algol and E is an integer number. This constant orbital period ephemeris "model" is quite accurate, because all O-C values are between −0. d 24 and +0. d 15 during 236 years. Out of all 2224 estimates, only 197 have an error estimate, and none of the 1236 first ones. However, this does not mean that these values without error estimates are unreliable or inconsistent. The error estimates are available for only about 10% of data. These are all new observations after the year 1921. The range of these known errors is between 0. d 0002 and 0. d 013. The most accurate TIDAK database O-C values have four decimals. Since the errors are not known for over 90% of observations, we use arbitrary errors σ i = 0. d 00010 for all O-C values. These arbitrary numerical values do not influence our results, because we use the same weight w i = σ −2 i for every observation, and we compute the DCM test statistic z from the sum of squared residuals (Eq. 8). We will also show that a weighted DCM search, where the O-C data accuracy improves towards modern times, does not alter our results (Sect. 5.4) .

We also analyse shorter subsamples of all data (Table  A5 ). In Sect. 5.2, we apply DCM to the first 226.2 years of all data (First226 y -data). This gives us a prediction for the last 9.2 years of all data (Last9 y -data). In Sect. 5.2, the same DCM procedures are also applied to the first 185.5 years of all data (First185 y -data), and the last 50 years of all data (Last50 y -data).

The Discrete Chi-Square Method (DCM) notations for the data are y i = y(t i ) ± σ i , where t i are the observing times and σ i are the errors (i = 1, 2, ..., n). The time span of data is ∆T = t n − t 1 . The mid point of data is t mid = t 1 + ∆T /2.

We analyse these data with DCM, which can detect many signals superimposed on arbitrary trends. Detailed instructions for using the DCM python code were given in the appendix of Paper I. In this current study, we provide all necessary information for reproducing our DCM analysis of Algol data. 1 DCM model is g(t) = g(t, K 1 , K 2 , K 3 ) = h(t) + p(t).

(2)

It is a sum of periodic and aperiodic functions

model non-linear, because all free parameters are not eliminated from all partial derivatives ∂g/∂β i . If thē β I frequencies are fixed to constant known tested numerical values, the model becomes linear, because all partial derivatives ∂g/∂β i no longer contain any free parameters. In this case, the solution for the remaining second group of free parameters,β II = [B 1,1 C 1,1 , ..., B K1,K2 , C K1,K2 , M 0 , ..., M K3 ], is unambiguous. We refer to this type of models and their free parameter solutions, when we use the concepts "linear model" and "unambiguous result". DCM model residuals are

For every combinationβ I = [f 1 , f 2 , ..., f K1 ] of tested frequencies, we compute the DCM test statistic

from the sum of squared residuals R = n i=1 2 i of a non-weighted linear model least squares fit. We use this non-weighted test statistic, because the errors for the data are unknown.

The global periodogram minimum is at z min = z(f 1,best , f 2,best , ..., f K1,best ),

where f 1,best , f 2,best , ..., f K1,best are the frequencies of the best DCM model for the data. Every scalar value of this z periodogram is computed from K 1 frequency values. For example, the K 1 = 2 periodogram could be plotted like a map, where f 1 and f 2 are the coordinates, and z = z(f 1 , f 2 ) represents the height. However, a graphical presentation for K 1 ≥ 3 is impossible, because it requires more than three dimensions. In Paper I, we solved this problem by presenting only the following one-dimensional slices of the full periodograms z 1 (f 1 ) = z(f 1 , f 2,best , ..., f K1,best )

z 2 (f 2 ) = z(f 1,best , f 2 , f 3,best , ..., f K1,best )

z 3 (f 3 ) = z(f 1,best , f 2,best , f 3 , f 4,best , ..., f K1,best ) (10)

z 4 (f 4 ) = z(f 1,best , f 2,best , f 3,best , f 4 , f 5,best , f K1,best )

z 5 (f 5 ) = z(f 1,best , f 2,best , f 3,best , f 4,best , f 5 , f K1,best )

z 6 (f 6 ) = z(f 1,best , f 2,best , f 3,best , f 4,best , f 5,best , f 6 ).

In the above K 1 = 2 map analogy, z 1 (f 1 ) would represent the height at f 1 coordinate when moving along the constant line f 2 = f 2,best that crosses the global minimum z min . DCM determines the following h i (t) signal parameters We determine the DCM model parameter errors with the bootstrap procedure (Efron & Tibshirani 1986; Efron & Tibshirani 1994) . During each bootstrap round, we select a random sample¯ * from the residuals¯ of the DCM model (Eq. 7). Each i can be chosen as many times as the random selection happens to favour it. This gives the artificial bootstrap data sample

DCM model for eachȳ * sample gives one estimate for every model parameter. For each particular model parameter, its error estimate is the standard deviation of all estimates obtained from allȳ * bootstrap samples. We have already used this same bootstrap procedure in our TSPA-and CPS-methods (Jetsu & Pelt 1999; Lehtinen et al. 2011 ). Finally, we note that our bootstrap procedure can not assess the bias in the y i input data, which first contaminates the i values, and then also the i and y i values. We use the Fisher-test to compare any pair g 1 (t) and g 2 (t) of simple and complex models. Their number of free parameters (η 1 < η 2 ), and their sums of squared residuals (R 1 , R 2 ) give the test statistic

Our null hypothesis is H 0 : "The complex model g 2 (t) does not provide a significantly better fit to the data than the simple model g 1 (t)."

Under H 0 , the test statistic F R has an F distribution with (ν 1 , ν 2 ) degrees of freedom, where ν 1 = η 2 − η 1 and ν 2 = n − η 2 (Draper & Smith 1998) . The probability for F R reaching values higher than F is called the critical level Q F = P (F R ≥ F ). We reject the H 0 hypothesis, if

where γ F is the pre-assigned significance level. It represents the probability of falsely rejecting H 0 when it is in fact true. The H 0 rejection means that we rate the complex g 2 (t) model better than the simple g 1 (t) model. The Q F critical level becomes smaller when F R increases. In other words, the H 0 hypothesis rejection probability increases for larger F R values. The basic idea of the Fisher-test is simple. The sum of complex model residuals R 2 decreases when the η 2 number of free parameters increases. When the complex model has more η 2 free parameters, the first (R 1 /R 2 − 1) term increases F R (Eq. 11), but at the same time the second (n−η 2 −1)/(η 2 −η 1 ) penalty term decreases F R . In conclusion, this second penalty term prevents overfitting.

The key ideas of DCM are 

3. The f 1 > f 2 > ... > f K1 grid combination of the best DCM model minimizes the z test statistic.

4. The bootstrap method gives the error estimates for all model parameters.

5. All different K 1 , K 2 and K 3 order nested models are compared using the Fisher-test, which reveals the best one of all models (Draper & Smith 1998; Allen 2004 ).

In short, DCM applies the following robust and well tested statistical approaches: Linear least squares fits (Idea 1), χ 2 and R test statistic (Idea 2), Dense tested frequency grids (Idea 3), Bootstrap utilizing residuals (Idea 4) and Fisher-test comparison of nested models (Idea 5).

The caveats of DCM are 1. DCM is designed for periodicity detection, but it gives no direct significance estimates for these detected periodicities. In this sense, DCM resembles our former TSPA-and CPS-methods (Jetsu & Pelt 1999; Lehtinen et al. 2011) . DCM utilizes indirect Fisher-test significance estimates to identify the best model among all tested models, but it gives no significance estimates for the detected periodicities of this best model. We will later discuss our indirect significance estimates, especially in connection with the look-elsewhere effect (Sect. 6.6).

2. The best frequency combination can be missed if the tested grid is too sparse (Idea 3). However, an adequately dense tested frequency grid eliminates the possibility for this kind of an error. The caveat is that denser grids require more computation time.

For example, all three signal z 1 , z 2 and z 3 periodograms for the original data are continuous and display no abrupt jumps, because the periodogram values for all close tested frequencies correlate (see Fig. A6 ). Since the frequencies of the minima of all these periodograms are accurately determined, there is no need to test an even denser grid (i.e. more trials), because this would not alter the final result of the non-linear iteration (Paper I: Eq. 18). In other words, the detected period values would no longer change, if we increased the number of of tested frequencies (Paper I: n L and n S trials). Since DCM gives no direct significance estimates for the detected periods (Caveat 1), there is no need to determine the number of independent trials, like for example the number of independent tested frequencies (e.g. Jetsu & Pelt 2000, their Eq. A.1).

3. If the grid of each tested f 1 > f 2 > ... > f K1 frequency contains n f values, the total number of tested frequency combinations is proportional to ∝ n K1 f . For example, it took about one month for an ordinary PC to compute the four signal DCM model 4,2,1 search, and to analyse its twenty bootstrap samples (Table A8 , model M=4).

4. Some DCM models are unstable because they are simply wrong models for the data. For example, a wrong p(t) trend order K 3 , or a search for too many K 1 signals, can cause such instability. In this paper, we denote such unstable models with "Um". We denote the two signatures of such unstable models with "Ad" = Dispersing amplitudes = Amplitudes and/or amplitude errors disperse. "If" = Intersecting frequencies = At least two model frequencies are too close to each other.

We give list all our symbols in Table A6 . Both of the above instabilities were defined in Paper I (Sect. 4.3), where a typical example of the wildly oscillating signals was also shown in Fig. 6 of Paper I.

DCM tests all reasonable alternative linear models for the data, and determines the unambiguous results for the best values of their free parameters. This brute numerical approach finds the best model among all alternative models. DCM "works like winning a lottery by buying all lottery tickets" (Paper I).

The light-time travel effect (LTTE) caused by a third body is (Irwin 1952) . This relation gives EB orbit around the common centre of mass of all three stars. The orbit parameters are the semimajor axis ([a] = AU), the orbital plane inclination ([i] = rad), the eccentricity of orbit (e), the longitude of periastron ([ω] = rad), the true anomaly ([ν] = rad) and the amplitude of light-time travel effect

which is half of the peak to peak amplitude A of the observed O-C changes ([A] = d) .

We compute the true anomaly from the Fourier expansion where

is the mean anomaly (Mueller 1995; Roy 2005 

where m 1 and m 2 are the masses of EB (Wolf et al. 1999; Zasche & Wolf 2007; Manzoori 2016; Esmer et al. 2021 ). The semi-major axis of the third body orbit is

where a = 173.15(A/2)/ sin i. For circular third body orbit, the suitable O-C curve DCM model order is K 2 = 1, the pure sinusoid (Eq. 13: e = 0). For an eccentric e > 0 third body orbit, the O-C curve is not a pure sinusoid, and the suitable DCM model order is K 2 = 2 (Hoffman et al. 2006) .

Here, we present separately the DCM period search results for all data (Sect. 5.1), First226 y -data (Sect. 5.2) and First185 y -data (Sect. 5.3). We also make some additional experiments (Sect. 5.4).

In Table A7 , the Fisher-test is used to compare the results for all data in twelve separate DCM period searches between P min = 6000 d and P max = 80000 d . These models have one, two or three signals (K 1 = 1, 2 or 3). The third body orbits can be eccentric (K 2 = 2 ≡ e > 0). The alternative tested p(t) trends are K 3 = 0, 1, 2 or 3. Table A7 contains many notations "-", because it makes no sense to compare the same pair of models twice, nor to compare the model to itself. The total number of compared pairs is (12 × 11)/2. For example, the Fisher-test comparison of the one signal M=1 and M=2 models gives a large test statistic value F = 2821. The critical level Q F of this F value falls below the computational 2 accuracy of 10 −16 (Table A7 : Q F <10 −16 ). This means that the linear K 3 = 1 trend model 1,2,1 is absolutely certainly a better model than the constant K 3 = 0 trend model 1,2,0 . The upward arrow "↑" indicates this result. Note that Table A7 contains numerous "Q F <10 −16 " cases, where the identification of the better model is absolutely certain.

All column M=10 arrows point upwards (↑), and all line M=10 arrows point leftwards (←) in Table A7 . Hence, this stable M=10 model is better than all other eleven alternative models. This best DCM model 3,2,1 for all data is a sum of K 1 = 3 signals having an order K 2 = 2, and a linear K 3 = 1 trend. We use this K 3 = 1 linear trend in all analysis of original data. The meaning of this linear trend is discussed later (Sect. 6.4, Eq. 27). We will also show that all data contains only three K 2 = 2 order signals between 8000 and 80000 days (Sect. 5.1.2).

Four of the twelve models are unstable "Um" (Table  A7 : M= 3, 5, 8 and 9). There are three models, where the detected period exceeds ∆T time span of data (Table  A7 : M= 2, 3 and 7). They are denoted with the symbol "Lp" = Leaking period = At least one detected period exceeds ∆T time span of data.

In Table A7 , we compared (12 × 11)/2 pairs of models against each other. The better model in each pair was identified with the Fisher test: the complex model above "↑", or the simple model on the left "←".

The structure of our next Table A8 is more complicated, because we squeeze all DCM eccentric orbit search results for all data into this single table. We search for periods between 8000 and 80000 days. The left side of this table gives the detected periods and amplitudes. The right side gives the Fisher-test comparison results. For example, the one signal M=1 model period and amplitude are P 1 = 88183 d ± 816 d and A 1 = 0. d 313 ± 0. d 004. The next six "-" notations for this M=1 model mean that it has no other periods P 2 , P 3 or P 4 , nor amplitudes A 2 , A 3 or A 4 . Fisher-test comparison between this one signal model 1,2,1 (M=1) and the two signal model 2,2,1 (M=2) gives an extreme test statistic value F = 183. The critical level Q F <10 −16 confirms that M=2 model is certainly the better one in this pair of models. Comparison of M=1 model to M=3 and M=4 models gives the same result.

For the next M=2, 3 and 4 models, the number of detected periods and amplitudes increases one by one. The number of Fisher-tests decreases one by one, because it is unnecessary to test the same pair of models twice "-", nor to compare any model to itself "-".

The periods and amplitudes for one, two and three signal models are consistent (Table A8 : M=1-3). When we detect a new signal, we re-detect the same old earlier signal periods and amplitudes for models having less signals. The one signal M=1 model shows a leaking period "Lp", because the P 1 = 88183 d period exceeds the ∆T = 86171 d time span of data. The two and the three signal M=2 and M=3 models are stable, but the M=4 model is not "Um".

The one-dimensional z 1 (f 1 ), z 2 (f 2 ), z 3 (f 3 ) and z 4 (f 4 ) periodogram slices (Eq. 10) of M=4 model are shown in Fig. A4 . The transparent diamonds denote locations of the red z 1 (f 1 ), the blue z 2 (f 2 ) the green z 3 (f 3 ) and the yellow z 4 (f 4 ) periodogram minima. These minima are clearly separated.

The four signal M=4 model is unstable, because it suffers from the amplitude dispersion "Ad" effect. The periodograms of this model do not betray this effect (Fig. A4 ), but the exceedingly high amplitude green h 3 (t) and yellow h 4 (t) signals do (Fig. A5) . The errors of both A 3 and A 4 amplitudes are large. The P 4 = 55172 d period is about two times longer than the P 3 = 26846 d period. DCM exploits the anti-phase sum of these two dispersing high amplitude signals for modelling all data.

The stable three signal M=3 model is a better model for all data than the failing unstable "Um" four signal M=4 model. Fisher-test reveals with an absolute certainty of Q F <10 −16 that this three signal M=3 model is also better than the M=1 model or the M=2 model (Table A8 : two times "↑" in Col 8).

Model M=3 periodogram minima are also clearly separated (Fig. A6 , lower panel). When all three periodograms are plotted in the same scale, the two z 1 (f 1 ) and z 2 (f 2 ) periodogram minima appear to be shallower than the z 3 (f 3 ) periodogram minimum, because the high amplitude h 3 (t) signal dominates in this three signal M=3 model (Fig. A6 , upper panel). This 79999 d period h 3 (t) signal has a much bigger impact on the sum of squared residuals R than the two lower amplitude 20358 d period h 1 (t), and 24742 d period h 2 (t) signals. This three signal M=3 model is shown in Fig. A7 . The level of residuals, denoted by blue dots, is stable and there are no trends. Each h j (t i ) signal

is also shown separately (Fig. A8 ). The red h 1 (t) and the blue h 2 (t) curves show two minima and two maxima, but the green large amplitude h 3 (t) curve shows only one minimum and one maximum. It takes about one month for an ordinary PC to compute the results for the four signal M=4 model, as well as to analyse at least twenty bootstrap samples (Table  A8 : model 4,2,1 ). The computation of five signal model would take several months. "Fortunately", there is no fourth or fifth signal between 8000 and 80000 days in all data, because the M=4 model is unstable "Um". The three signal M=3 model is the best model for all data. Therefore, we can search for additional periods shorter than 8000 days from the M=3 model residuals.

Since the M=3 model residuals contain no trends, we analyse them by using K 3 = 0 models having a constant p(t) level. The period search between 500 d and 8000 d gives two new periods 680. d 4 and 7290 d (Table A8 , model M=6). In the three signal M=7 model, the periods P 2 = 7124 d and P 3 = 7698 d give

which is equal to the time span ∆T = 86171 d of all data (Table A8 : M=7). In other words, the difference between the real P 2 = 7124 d and the spurious P 3 = 7698 d period is one round during ∆T . Our symbol for this type of spurious periods is "Sp" = Spurious period = Unreal periods caused by data time span and real periodicity. Therefore, we reject the M=7 model, and the best model for residuals is the M=6 model. In this analysis of residuals, DCM again consistently re-detects the same periods and amplitudes of earlier models having less signals. Model M=6 periodograms, and the model itself, are shown in Figs. A9 and A10. The two last 680. d 4 and 7290 d signals detected from the residuals are shown in Fig. A11 . As expected of a real O-C signal, both curves have only one minimum and one maximum. These two signals are 44.8 and 41.0 times weaker than the strongest first detected 79999 d signal.

For the original data, DCM detects simultaneously the three signals signals and the trend of M=3 model. For the residuals, the same applies to the two signals and the trend of M=6 model. In this sense, DCM differs from the "pre-whitening" technique, which requires that the trend must be determined and removed before even one signal at the time can be detected (e.g. Reinhold et al. 2013) . This "pre-whitening' technique, which applies the Discrete Fourier Transform (DFT), was compared to DCM in Paper I (Sect. 6).

We conclude that DCM detects five signals from all data (n = 2224). The full model for all Algol's O-C data is the sum of the M=3 model for the original data, and the M=6 model for the residuals (Table A8 ). Our notation for this sum model 3,2,1 + model 2,2,0 of two models in Table A8 is simply the "M=3+6 model". This model is denoted with the green continuous line in Figs. 1ab. Its standard deviation of residuals is 0. d 011. We also give a ten year prediction for Algol's O-C changes after our last observation on Oct 18th, 2018 (Fig. 1b) . 

In our appendix, we show that if an eccentric orbit O-C curve has a period p, then this curve is a sum of two circular orbit O-C curves having periods p and p/2. For this reason, the DCM period search results obtained for circular orbits in this section can be used to check the eccentric orbit results presented earlier in Sect. 5.1.2, and vice versa (Table A14) .

For third body circular orbit, the correct DCM model h i (t) signal order is K 2 = 1 (Eq. 13: e = 0). We fix the p(t) trend to K 3 = 1, and search for the correct K 1 number of circular orbit sinusoidal signals in all data. Two alternative approaches are tested. We will show that both approaches give the same results.

In the first alternative approach, we search for one, two, three and four sinusoidal circular orbit signals having periods between 8000 and 80000 days in all data (Table A9 ). The one signal M=1 model is stable. The two, three and four signal M=2, M=3 and M=4 models are unstable "Um", because they all suffer from dispersing amplitudes "Ad". The largest periods ("Lp") in these three models exceed the all data time span

From the M=4 model residuals, we detect the fifth sinusoidal signal period 10175 d (Table A9 : M=5). The next M=6 model is unstable "Um", and it is also rejected with the Fisher-test criterion (Eq. 12).

DCM detects signatures of five sinusoidal signals having periods longer than 8000 days. Therefore, we search for shorter periods from the M=5 model residuals. This reveals three additional sinusoidal M=9 model signals (Table A9 ). The next four signal model M=10 is rejected with the Fisher-test criterion (Eq. 12).

In our first alternative approach, the best circular orbit model is the M=4+5+9 model (Table A9) .

Our typical number of tested periods is n L = 80 in the long search, and n S = 40 in the short search. We use these dense grids to eliminate the "trial factor" error (Sect. 3: Caveat 2). Computation time is proportional to ∝ n K1 L and ∝ n K1 S . For larger number of signals, these dense tested grids of ours take a long time to compute. For example, the computation of four signal model for all data, and its twenty bootstrap samples, takes about one month for an ordinary PC.

In the second alternative approach we also search for circular orbit periods between 8000 and 80000 days. However, we reduce the computation time dramatically by testing only n L = 30 and n S = 8 frequencies. In this case, an ordinary PC can perform the six signal DCM search in about one week. Unlike in the first alternative approach, we do not need to search for the fifth and sixth signal from the four signal model residuals. We can perform the five and the six signal DCM search directly to all original data. The four, five and six signal circular orbit model results for all original data are given in Table A10 . All M=1, 2 and 3 models suffer from amplitude dispersion "Ad", as well as from leaking periods "Lp", because their largest detected periods exceed ∆T . Model M=3 also suffers from intersecting frequencies "If". We reject it with the Fisher-test criterion (Eq. 12). The best circular orbit model for all original data is the five sinusoidal signal M=2 model. The M=2 model periodogram is shown in Fig. A12 . The periodogram minimum of the largest P 5 = 120740 d period is real, because the violet z 5 (f 5 ) curve in the lower panel turns upwards at smaller tested frequencies (i.e. periods larger than ∆T ). The M=2 model itself is shown in Fig. A13 .

From model M=2 residuals, we find two periods shorter than 8000 days (Table A10 : M=5). We reject model M=6, because the periods P 1 = 7034 d ± 148 d and

Hence, the spurious "Sp" period P 1 is connected to the real period P 2 and the time span ∆T = 86171 d of all data.

The second alternative approach best circular orbit model is the M=2+5 model (Table A10) .

We compare the results of our two alternative approach circular orbit DCM analyses in Table A13 . All results are consistent. The periods and amplitudes agree within their error limits. We detect the same five longer sinusoidal signal periods from the original data, and the same two shorter period sinusoids from the residuals. We get these consistent results even after dramatically reducing the number of tested frequencies. Hence, these two analyses not suffer from the "trial factor" effect (Sect. 3: Caveat 2). The dispersing amplitudes "Ad" or the leaking periods "Lp" do not either mislead this analysis.

The eccentric orbit DCM search results for subsample First226 y -data are given in Table A11 . The one signal M=1 model and two signal M=2 model suffer from leaking periods "Lp". The stable three signal M=3 model is the best one for the original data, because the four signal M=4 model is unstable "Um". For the M=3 model residuals, the best model is M=6 model. We reject model M=7, because the relation [P −1 2 − P −1 3 ] −1 = 81963 d ± 13594 d reveals that the third P 3 = 7757 d period is a spurious "Sp" period connected to the real period P 2 = 7078 d and the time span ∆T = 82602 d of data.

The best model for First226 y -data is the M=3+6 model (Table A11) . This model is shown in Fig. 2 . It gives an excellent prediction for the next nine years of Last9 y -data (Fig. 2b) . The standard deviation of prediction residuals is only 0. d 0078 (n = 50). It is smaller than the standard deviation 0. d 011 of the predictive M=3+6 model residuals (n = 2174). However, the larger errors of the older observations can explain this contradiction. The main conclusion is that our nine years prediction succeeds. The eccentric orbit DCM search results for the shortest subsample First185 y -data are given in Table A12 . The one signal M=1 model is stable. The two and three signal M=2 and M=3 models are unstable (Table A12: "Um"). The best model for First185 y -data is the stable four signal M=4 model. For the M=4 model residuals, the stable M=6 model is the best one, because the M=7 model is unstable "Um". The best M=4+6 model for First185 y -data is shown in Fig. 3 . Our fifty years prediction succeeds only for the first few years (Fig. 3b) . However, this is no surprise, because the time span of predictive data is only ∆T = 67680 d = 185 y . For this reason, the longest and the strongest detected predictive signal period is P 4 = 62992 d = 172 y (Table A12 : M=4). This high amplitude signal determines the long-term prediction trend for Last50 y -data. We have already shown that the correct period for this long-term trend would be 219 y (Ta-ble A8: M=3, Table A11 : M=3). The short 185 y time span of First185 y -data prevents the detection this correct 219 y period. The correct 219 y signal trend turns upwards slower than the wrong 172 y signal trend. This is the simple reason for the failure of our fifty years prediction for First185 y -data.

The Last50 y -data prediction error for shows a peculiarity that seems to defy the laws of statistics. First, the ±3σ prediction error increases, as one would expect ( Fig. 3b: green dotted lines) . Surprisingly, this prediction error then begins to decrease, and the prediction becomes very accurate close to HJD 2450000. After this, the prediction error begins to increase again. This peculiarity certainly requires an explanation.

The reason of this peculiarity could already be inferred from the black interference curve in Fig. A3 (lowest right panel: P 1 = 24771). The scatter of g(t) interference curve is not the same at all phases. In this particular case, this scatter increases close to the maxima, but it decreases close to the minima. The largest and the smallest scatter coincides with the phases when the first time derivative fulfillsġ(t) = 0.

However, the above mentioned effects in Fig. A3 are caused by interference of only two signals, while the peculiar error limit effect in Fig. 3 occurs in the M=4+6 model sum of six signals. We show this model for twenty bootstrap samples in Fig. A14 (red dotted curves). The scatter of these curves increases when the predictive data ends at the dotted black vertical line. However, all dotted red curves converge close to the vertical continuous black line at HJD 2450000. After this line, they diverge again. Before this line, the data shows an increasing trend, but the positive slope is decreasing ( Fig.  3a : red circles). A suitable model would beġ(t) > 0 andg(t) < 0. After this line, this slope is still positive, but it is increasing. Now the suitable model would bė g(t) > 0 andg(t) > 0. This means that there is a turning pointġ(t) = 0 close this HJD 2450000 epoch, where thë g(t) sign changes from negative to positive. The second derivative sign change of any function forces this function to change its direction twice. This M=4+6 model turning point forces the bootstrap model solutions to converge. This simple effect explains why the prediction error increases, decreases, and again increases ( Our turning point hypothesis would explain the gap in O-C data close to HJD 2450000 ( Fig. A14 : vertical continuous line). There are no such gaps in Algol's modern O-C data, not even during the two World Wars. TIDAK database contains only four O-C values between HJD 2448288 and HJD 2449988 (≡ 4.6 years). Even today, one of these four is still marked "unpublished" (1997, Drozdz: HJD 2449317.4171 ). Close to the above mentioned turning point, the O-C data did no longer support the well established expected long-termġ(t) > 0 andg(t) < 0 trend. Perhaps for this reason, the contradictory new data was not published at that time. Only when the newġ(t) > 0 andg(t) > 0 trend was securely established, the continuous flow of supporting O-C observations began again.

We conclude that, except for the first few years, our Last50 y -data prediction fails. However, our turning point epoch prediction HJD 2450000 is excellent.

We divide all original data into two parts. Both halves are too short for the detection of the long 219 years period. This hampers their period analysis. In the first low accuracy half, we detect only one signal of about 137 years. From the more accurate second half, we detect four signals of 1.86 30.9, 39.7 and 103.3 years. The shortest one is equal to the orbital period of Algol C.

We also test two alternatives, where the weights of observations increase linearly. In two alternative experiments, the weights are doubled or quadrupled during the time span of all data. In both cases, the five strongest signals detected from the weighted data are identical to those detected from non-weighted data (Table A8 , M=3+6 model).

The eccentric orbit analysis indicates that all data contains five signals (Table A8, M=3+6) . Here, we argue that the correct number of signals may also be six. We use bold letters p 1 , p 2 , p 3 , p 4 , p 5 and p 6 for the periods of these signals (Table 1) . This notation helps the readers to separate these six periods from the numerous other P 1 , P 2 , ..., P 6 , p, p 1 , p 2 , p 3 and p periods. We use the tentative names Algol C, Algol D, Algol E, Algol F, Algol G and Algol H for the objects possibly connected to these periods. The corresponding peak to peak amplitudes are A 1 , A 2 , A 3 , A 4 , A 5 and A 6 .

Our six signal argument relies on two tables. The first table compares the eccentric and circular orbit analysis periods for all data (Table A14 ). The second table compares the periods detected in three different samples: All data, First226 y -data and First185 y -data (Table A15) .

In our Appendix, we apply DCM to simulated O-C data (Eq. 13). We show that the following four different effects are encountered when the O-C data contains one period p, or two periods p 1 and p 2 .

"Correct-p": DCM detects the correct period p.

"Half-p": DCM detects the spurious period p/2. (Table A12 ). Otherwise as in Fig. 1a. (b) Prediction for Last50 y -data. Otherwise as in Fig. 1b. "Double-p": DCM detects the spurious period 2p.

"Interference-p ": DCM detects the spurious period p caused by p 1 and p 2 interference (Eq. A7).

The "Half-p" and "Double-p" effects can mislead DCM analysis of low eccentricity O-C curves, which resemble pure sinusoids.

There is only one minimum and one maximum in the real O-C curve caused by the LTTE of a single third body. This third body can approach and recede only once during one orbital period p. Hence, the O-C "p interference" curves having two minima and two maxima can not be caused by one body alone, but they may indicate the presence of more than one body.

In the next Sects. 5.5.1-5.5.5, we illustrate one p 1 , p 2 , p 3 , p 4 , p 5 and p 6 signal at the time, how the above mentioned four effects can explain all eccentric and all circular orbit DCM period search results.

The circular orbit signal period P c,7 = 120740 d ± 41002 d differs about ±1σ from the eccentric orbit period p 6 = P e,5 = 79999 d ± 1216 d (Table A14) . Hence, the circular and eccentric orbit analyses give the same correct p 6 period ("Correct-p" effect).

This p 6 period is two times longer than the next circular orbit period P c,6 = 42422 d ± 640 d ("Half-p" effect). The p 6 = 219 y signal curve shows only one minimum and one maximum ( Fig. A8 : lowest panel green curves), because the two strongest circular orbit P c,7 and P c,6 signals are "in phase". These results confirm that DCM succeeds in detecting the p and p/2 regularities illustrated in Fig. A1 and Table A3 .

DCM detects the p 6 = 219 y signal in all data and First226 y -data (Table A15 ). The too short First185 ydata time span prevents the detection of the p 6 period. Therefore, the largest detected P 4 = 62992 d ±2499 d period differs more than ±3σ from p 6 .

We use an amplitude estimate A 6 = A e,5 = 0. d 287 ± 0. d 005 for this p 6 = 219 y signal (Table A14) The connection between the eccentric orbit p 5 = P e,4 = 24742 d ±141 d signal and the circular orbit P c,5 = 24747 d ±872 d signal is definitely the "Correct-p" effect (Table A14 ). The "Half-p" effect certainly connects this p 5 signal also to circular orbit P c,4 = 12294 d ±109 d signal.

However, two questions need to be answered. Why does the p 5 = 66 y .4 signal show two minima and two maxima ( Fig. A8 : mid-panel blue curves)? This is impossible for any single third body eccentric orbit. Why are the A c,5 and A c,4 amplitudes of the two circular orbit P c,5 and P c,4 signals practically equal (Table A14) ?

The easiest answer to both questions would be that the p 5 = 66. y 4 and p 4 = 33. y 7 signals represent two separate independent signals, which are "off-phase". Their "Interference-p " effect could induce the two unequal minima and two unequal maxima of the blue O-C curve (Fig. A8) , which resembles the black interference curve in Fig. A3 . In this case, the circular orbit P c,4 = 12294 d ±109 d signal could represent a real fourth independent p 4 = 33. y 7 signal.

The p 5 = 66. y 4 signal is detected in all data and First226 y -data (Table A15 ). This p 5 signal is not detected in the shortest First185 y -data sample, but the p 4 = 33. y 7 signal is. We conclude that the p 5 = 66. y 4 and p 4 = 33. y 7 signals are most probably two independent real signals.

The amplitudes of the circular orbit P c,5 and P c,4 signals give our A 5 = A c,5 = 0. d 018 ± 0. d 002 and A 4 = A c,4 = 0. d 018 ± 0. d 001 amplitude estimates for the p 5 and p 4 signals (Table A14) .

Here, we have shown that the eccentric orbit p 5 = 66. y 4 signal may arise from the "Interference-p " effect of two circular orbit p 5 = 66. y 4 and p 4 = 33. y 7 sinusoids. Later, we will present an alternative explanation ( Fig.  A15: Configurations 2 and 3) .

None of the eccentric orbit periods is close to the circular orbit period P c,3 = 10144 d ± 30 d = 27 . y 8 ± 0. y 1 (Table A14) . However, the "Double-p" effect certainly connects this P c,3 period to the eccentric orbit period P e,3 = 20358 d ± 128 d . This P e,3 signal shows two maxima and two minima ( Fig. A8 : lower panel red curves). These two equal maxima and two equal minima are symmetric. This kind of symmetry is detected in our simulations of low eccentricity spurious double sinusoids (Table A2: "Dp"≡"Double-p" effect). Therefore, the P c,3 period probably represents a real signal p 3 = 27. y 8.

The eccentric orbit P e,3 = 20358 d signal is detected in all data, First226 y -data and First185 y -data (Table  A15 ). This means that DCM detects the p 3 ≈ P e,3 /2 signal in all these three different samples.

Our amplitude estimate for this p 3 = 27. y 8 signal is (Table A14) .

In this section, we have shown that the eccentric orbit P e,3 = 56. y 0 signal probably represents the "double wave" of the p 3 = 27. y 8 signal. We will later present an alternative explanation ( The eccentric orbit p 2 = P e,2 = 7269 d ±29 d signal and the circular orbit P c,2 = 7395 d ± 37 d signal are certainly connected (Table A14 : "Correct-p" effect).

Like any real third body O-C curve, this p 2 = 20. y 0 signal shows only one minimum and one maximum (Fig. A11: lower panel blue curves). DCM detects this p 2 = 20. y 0 signal in all data and First226 y -data (Table A15 ). In shortest First185 y -data sample, this p 2 period may be connected to its double period P 3 = 15429 d ± 222 d (Table A15 : "Double-p" effect).

Our amplitude estimate for this p 2 = 20. y 0 signal is A 2 = A e,2 = 0. d 007 ± 0. d 001 (Table A14) . The eccentric orbit and circular orbit DCM searches give the same p 1 = 680. d 4 ± 0. d 4 = 1. y 863 ± 0. y 001 signal (Table A14 : "Correct-p" effect).

DCM detects this p 1 = 1. y 863 signal from all three samples (Table A15) . Like any real O-C curve, this signal shows only one minimum and one maximum ( Fig.  A11 : higher panel red curves).

We use A 1 = A e,1 = 0. d 0064 ± 0. d 0007 (Table A14) . This signal is discussed later in greater detail (Sect. 6.3). 5.5.6. Two weakest signals DCM detects indications of two additional weaker signals P c,2 = 2986 d ± 39 d (Table A13) and P 2 = 3387 d ± 17 d (Table A15 ). They could be separate signals, because their ±3σ error limits do not overlap. They are 0.80 and 0.45 weaker than the weakest detected p 1 = 1. y 863 signal. We can not confirm whether these two weakest signals are real or spurious.

6. DISCUSSION Applegate (1992) mechanism can not explain the numerous strictly periodic O-C signals of Algol, because quasi-periodic activity cycles are never regular. Apsidal motion follows only one period. LTTE of Algol's companion candidates could cause these numerous strictly periodic cycles. Assuming circular orbits, we use m 1 = 3.7m and m 2 = 0.8m (Zavala et al. 2010 ) to compute the m 3 mass and the a 3 semi-major axis estimates for these tentative companion candidates (Table  1) . These approximate mass and semi-major axis estimates are obtained by assuming that each candidate is a "third" component. The effects of other candidates inside the orbit of the "third" component are ignored in Eqs. 18 and 19.

We call the eclipsing Algol A and Algol B pair the central eclipsing binary (cEB). Algol C is called a wide orbit star (WOS), as well as all other new tentative companion candidates. We use the same hierarchial system diagrams as Tokovinin (2021) .

Our first hierarchial system diagram shows the circular orbit i = 90 o inclination case of Table 1 (Fig. A15 : Configuration 1). The eight members in this configuration are cEB and six WOSs. The orbital periods WOS candidates are between 1.863 and 219.0 years. The most massive (m 3 = 2.50m ) companion candidate Algol H is also the most distant one (a 3 = 44.7AU). The four other WOS candidates are low mass stars (0.23m ≤ m 3 ≤ 0.43m ). The closest m i=90 3 = 1.16m companion candidate has an orbital period p 1 = 680. d 4±0. d 4, which is close to the known orbital period P orb = 679. d 85 ± 0. d 04 of Algol C (Zavala et al. 2010 ). We will discuss this probable detection of Algol C later in Sect. 6.3.

Our second hierarchial system diagram shows one alternative for Configuration 1 (Fig. A15: Configuration  2) . The seven members are cEB and five WOSs. We have already shown that the sum of "off-phase" sinusoidal p 5 = 66. y 4 and p 4 ≈ p 5 /2 = 33. y 7 signals can cause the p 5 = 66. y 4 period double wave (Sect. 5.5.2). However, a single p 5 = 66. y 4 long-period binary can cause a similar effect, if the masses of its members are unequal. These unequal masses could also explain the two unequal maxima and minima of the blue O-C curve in Fig. A8 . The red lines in our Configuration 2 diagram show this hypothetical long-period p 5 = 66. y 4 binary having an orbital period p 6 = 219. y 0 around the barycentre of the whole system (Fig. A15) .

Our third hierarchial system diagram is a minor modification of Configuration 2 (Fig. A15: Configuration  3) . The seven members are, again, cEB and five WOSs. Now we take the five periods of M=3+6 model as such. Signal 66. y 4 is not separated into two signals (Sect. 5.5.2). We use the full P e,3 = 55. y 8 signal period, not the half of this period (Sect. 5.5.3). This P e,3 = 55. y 8 signal could also represent a long-period binary, where the masses of both components are approximately equal. In Configuration 3, the two longperiod p 5 = 66. y 4 and P e,3 = 55. y 8 binaries orbit each other during p 6 = 219. y 0. This may be the most stable one of our three configuration alternatives, because the cEB and the remaining two inner orbit WOSs would only weakly perturb the two hypothetical long-period binaries, and vice versa. This type of quintuple binary systems have been discovered (e.g. Zasche & Uhlař 2013, their Fig. 2 of V994 Her).

In binaries, the radial velocity observations can reveal the presence of a third body, like in the discovery of Algol C (Curtiss 1908) . For nearby hierarchial systems, combined astrometric orbit and radial velocity observations can be used to solve their detailed structure (e.g. Tokovinin 2021). When Hajdu et al. (2019) searched for WOSs from the O-C data of 80 000 EBs, they detected 992 systems having one WOS, but only four systems possibly had two WOSs. Our DCM analysis of Algol's O-C data suggests the presence of five or six WOSs. These O-C data can not reveal a lot about the structure this hierarchial system, not even the exact number of stars ( Fig. A15: Configurations 1, 2 or 3) . However, we can give some ideas that may help in the detection of Algol's WOS candidates. Here, we assume that the candidate orbits are circular and their orbital plane inclinations are i 3 = 90 o (Table 1 ). In this Configuration 1, the observed maximum and minimum radial velocities of WOSs candidates are (2) Correct hierarchial system alternative is Configuration 1 (Fig. A15) . (3) All orbits are circular. (4) Every candidate can be treated as a "third body". In other words, effects of other candidates inside "third body" orbit can be ignored in Eqs. 18 and 19. where v 0 = 4.0 km/s is Algol's radial velocity (Wilson 1953) . The angular distance between Algol and its WOSs changes constantly. We compute these angular distance changes in Algol's cEB frame of rest. At the O-C curve minima and maxima, the largest distance changes are

during a time interval ∆t ≤ p 3 /2. For longer time intervals, we use ∆t = p 3 /2 which gives ∆a max = 2a 3 . The smallest

distance changes coincide with the O-C curve mean level. This relation holds for t 0 ≤ p 3 . For longer time intervals, we use ∆t = p 3 which gives ∆a min = 2a 3 . The proper motion of Algol is µ 0 = 2.49 mas/y (van Leeuwen 2007). The minimum and maximum proper motion of each candidate is

where µ c = ∆a max (∆t = 1 y ) is the maximum proper motion during one year. Note that µ min = 0 for every candidate, because their µ c > µ 0 . We emphasize that our ∆a min and ∆a max estimates refer to the candidate distance changes with respect to cEB, while our µ min and µ max estimates refer to the proper motion of all members in the sky.

All parameters of Eqs. 21-26 are given in Table A16 . The two estimates for ∆a min and ∆a max are computed for observations spanning 5 or 20 years. This information is useful for future searches of our Algol's member candidates.

In December 2020, the latest third Gaia data release (DR3) confirmed no certain detections ±4" around Algol, and only one certain ±40" detection. In their analysis of Gaia DR3 data, Torra et al. (2020) note that "most problems come from the bright sources and the strange image profiles." They rejected 8159.3 million bright sources, 158.0 million very bright sources and 4066.7 million odd window profiles. Algol is definitely "too bright". Its brightness profile is constantly changing due to the movement of the known members Algol A, Algol B and Algol C, let alone due to the primary and secondary eclipses. Therefore, Gaia could not have measured the positions and movements of the objects in our Table A16 .

Algol H candidate would be easiest to detect, because its distance from the cEB is the largest. This most massive candidate is very probably also the brightest candidate. At the moment, its O-C curve is close to the mean level ( Fig. A8 : left hand lowest panel green h 3 (t) curve). Hence, Algol H would be close to its projected maximum a 3 = 1569 mas distance from cEB. The cEB is receding from us because its O-C values are increasing for the next fifty years. Currently, Algol H would be approaching us at its minimum radial velocity v min = −2 km/s (Table A16 ). The distance changes between cEB and Algol H would be small, only ∆a min = 4 or 64 mas during the next 5 or 20 years.

Direct interferometric images have been obtained of Algol A, Algol B and Algol C (e.g. Zavala et al. 2010; Baron et al. 2012 ). If Algol B and Algol C really are less massive than our distant Algol H candidate, why did the earlier interferometric imaging not reveal the presence of this massive candidate? Firstly, this Algol H candidate is about 20 times further away from the cEB than Algol C, which means that the area of interferometric imaging should have been about 20 × 20 = 400 larger. Secondly, this Algol H candidate could be a longperiod binary, where both members are much less massive and much dimmer than a single 2.50m star (Fig.  A15: Configurations 2 and 3) . One or two members of this long-period binary could be an evolved object, like a white dwarf. Thirdly, Zavala et al. (2010) and Baron et al. (2012) applied a three star model. Algol H contribution to their modelled total flux would have remained constant, because its position did not change during their observations (Table A16 ). We conclude that using an over 400 times larger imaging area, and a model of at least four stars, may lead to the interferometric detection of this distant Algol H candidate. The detection of the other four less massive candidates with this technique is much more challenging (Table 1 : i = 90 o , 0.23m ≤ m 3 ≤ 0.43m ). However, the ∆a min and ∆a max values of these less massive candidates show that their movements are easier to detect even during shorter periods of observations (Table A16) . Powell et al. (2021) studied the sextuple-eclipsing binary system TIC 168789840 with the speckle interferometry technique. They could resolve this hierarchial system of three eclipsing binaries. Their estimate for the outer period in this hierarchial system was about 2000 years. Algol is about twenty times closer to us than TIC 168789840 (d ≈ 570pc). The orbital period of our Algol H candidate is about 200 years. Hence, it might be possible to detect Algol H with the speckle interferometry.

DCM detects the weakest p 1 = 680. d 4 ± 0. d 4 signal in all three samples: all data, First226 y -data and First185 y -data. This p 1 signal is 44.8 times weaker than the strongest p 6 signal (Table 1) . DCM detects this weakest p 1 signal although it is buried under the interference of five stronger p 2 , p 3 , p 4 p 5 and p 6 signals, and a linear p(t) trend. The period of this p 1 signal differs only 1.4σ from the known orbital period P orb = 679. d 85 ± 0. d 04 of Algol C (Zavala et al. 2010) . This indicates that all other five detected stronger signals are real periodicities, but it does not irrefutably prove this idea. Our O-C data contains 127 rounds of Algol C around Algol AB, and this orbit is known to be stable (Zavala et al. 2010; Baron et al. 2012; Jetsu et al. 2013 ). Our lower limit for the mass of Algol C (Table 1 : i = 90 o and 1.2m ) is smaller than the interferometric estimates by Zavala et al. (2010, i = 83 . o 7 ± 0. o 1 and 1.5 ± 0.1m ) and Baron et al. (2012, i = 83 . o 66 ± 0. o 03 and 1.76 ± 0.15m ). This indicates that not even DCM can retrieve the full amplitude of this weak Algol C sig-nal when it is buried under five stronger signals and a linear trend.

All detected signals are strictly periodic, because they are also detected in the 9.2 years shorter subsample First226 y -data. Except for the p 2 and p 5 signals, the other four signals are also detected in the fifty years shorter subsample First185 y -data. This apparent absence of these two p 2 and p 5 signals in First185 ydata could be explained by the "Half-p" and "Double-p" effects (Table A15 ). However, strict periodicity alone does not prove that Algol's hierarchial system is stable.

The perturbations of WOS can cause periodic cEB orbital plane changes (Soderhjelm 1975, Eq. 27) . Such long-term orbital plane changes with respect to the line of sight may even stop the eclipses completely, or at least reduce the depth of eclipses, like in the case of AY Mus (Soderhjelm 1974). However, the cEB orbital plane is stable for Ψ = 0 o or 90 o , where Ψ is the angle between cEB and WOS orbital planes. This is the case for Algol C, the only currently known WOS of Algol (Baron et al. 2012, Ψ = 90 . o 20±0. o 32). No changes have been observed in the eclipses of Algol in modern times, and these events were most probably also observed over three thousand year ago (Jetsu et al. 2013) . This is possible only if all WOSs have Ψ = 0 o or 90 o . If the orbital planes of all WOS are co-planar, then all WOSs must have Ψ = 90 o , because this is the known case for Algol C. If all WOS orbit were not co-planar, this would certainly reduce the stability of this system, and perhaps also weaken or stop the observed eclipses.

The mass transfer from the less massive Algol B to the more massive Algol A should increase the orbital period (Kwee 1958, Eq. 5) . The numerous published mass transfer rate estimates range from 10 −13 m yr −1 to 10 −7 m yr −1 (Jetsu et al. 2013, Sect. 4 ). However, no regular long-term period increase has been observed since Goodricke (1783) discovered Algol's periodicity. All WOSs can also perturb the cEB by other physical mechanisms, like the Kozai effect (Kozai 1962), or the combination of Kozai cycle and tidal friction (Fabrycky & Tremaine 2007) . Against this background, our linear K 3 = 1 trend result for p(t) is surprising (Sect. 5.1.1). For 236 years, Algol's orbital period has been constant

where P 0 = 2. d 86730431 (Eq. 1) and M 1 = 0.1278 is p(t) coefficient for M=3 model in Table A8 . This causes the linear O-C change of 0. d 256 in Fig. A7 (upper panel: dotted line). It also means that LTTE effects alone can explain all observed O-C changes. No additional effects, like the quadratic K 3 = 2 trend caused by mass transfer, are needed to explain these O-C data.

In the future, long-term integrations may confirm the dynamical stability of this system. Currently, even the exact number of WOS candidates remains unknown, because three different hierachial system diagrams can explain the detected WOS periods (Fig. A15: Configurations 1, 2 and 3) . For any WOS period p 3 , the correct m 3 , e 3 , a 3 , i 3 , ω 3 and Ψ 3 initial value combinations for the long-term integrations are also unknown. Therefore, our O-C data can not give an unambiguous solution for this stability problem. Whether or not this system is stable, we can determine the p 3 periods that are observed today.

We admit that an unambiguous identification of all individual signals from the interference sum of numerous signals is not always possible. One example is the p 5 and p 4 signal identification in Sect. 5.5.2. However, this whole identification problem is irrelevant from the predictability point of view. The sum of identified signals is equal to the sum of unidentified signals. Both alternatives give the same prediction.

The linear and quadratic EB ephemerides can not predict the exact epochs of future eclipses (e.g. Kreiner Fig. 1 ). Different O-C subsets can give different periods, but this does not mean that there is something wrong with the period search methods themselves, like DCM. Our 9.2 years O-C prediction for Algol is based on First226 y -data (Fig. 2) . Strict periodicity can explain why this prediction succeeds. Predictability is impossible without strict periodicity. This prediction would fail, if even one of our detected signals were not strictly periodic, or if the K 3 = 1 linear p(t) trend were wrong.

Our next fifty years prediction is based on First185 ydata. Except for the first few years, this long-term prediction fails (Fig. 3) . The reason for this failure is simple. The longest 172 y period detected from First185 ydata is not correct. The short ∆T = 185 years time span of this sample prevents the detection of the correct signal period 219 y . This correct signal can be detected only from all data and First226 y -data. Together with the 0. d 26 trend p(t), this highest 0. d 29 amplitude dominating 219 y signal determines all long-term O-C predictions. The insignificant long-term trend contribution of all other weaker signals is always less than ±0. d 03, because the sum of their amplitudes is 0. d 06. Although our fifty years prediction for the O-C level fails (Fig.  3) , we get an excellent prediction for the turning point epoch at HJD 2450000 (Fig. A14) .

New O-C data after October 2018 can already be used to test our prediction for the next ten years (Fig. 1b) . These predictions should improve in the future, when all orbital period estimates become more accurate. Predictability should ultimately prove that all these signals are orbital periods. At the moment, we can not prove this. In the history of Astronomy, the seasons of the year posed a similar problem. Their one year periodicity was detected easily, but the reasons for it were understood much later: the orbit of the Earth around the Sun, and the tilted axis of Earth. However, it was possible to predict the seasons without understanding their origin. Our detected periods of Algol are certainly there, and for some reason or another they can be used to predict.

We test over thirty models having free parameters between η = 6 and 22 (Tables A7-A12 ). The total number of free parameters is even higher when the model for the original data is added to the model for the residuals. For example, the best M=3+6 model for all data has η = 17 + 11 = 28 free parameters (Fig. 1) . Our search for the correct model over a vast parameter space increases the probability for finding spurious apparently significant signals. This is called the "look-elsewhere effect" (e.g. Miller 1981; Bayer & Seljak 2020) . There are statistical methods that can account for the "lookelsewhere effect", and give direct significance estimates S for the periods of models having different degrees of freedom (e.g. Bayer & Seljak 2020, their Eq. 3.12) .

DCM applies Fisher-test to compare the significance of all pairs of simple and complex models. Fisher-test identifies the best model among all tested models (Eqs. 11 and 12). This approach does not account for the "look-elsewhere effect", because it gives no direct significance estimate S for the periodicities of this best model. Nevertheless, we can present several arguments indicating that the "look-elsewhere effect" has no significant impact on our results. 1. We apply the robust Fisher-test to compare any complex model having more signals than any simple model. We use the pre-assigned significance level γ F = 0.001 to reject the simple model (Eq. 12). This prevents over-fitting, because the probability that this best model selection fails is always smaller than one out of one thousand. In many cases, the extreme Q F <10 −16 critical levels confirm that the complex model is absolutely certainly Jetsu better than the simple model. This confirms that the data contain more signals than those present in the simple model. Our indirect Q F significance estimates confirm the presence of additional complex model periodicities, but they do not give us direct S significance estimates for these periodicities. Regardless of the "look-elsewhere effect", Fisher-test can confirm that the five signal M=3+6 model is the best model for all data.

2. The z periodogram values of close tested frequencies correlate and display no sudden jumps (see Sect. 3: Caveat 2). At some tested frequency grid density level, this means that the detected period values no longer depend on the number of tested periods (Paper I: n L and n S ). These unambiguous best period values are obtained from linear models. Increasing the number of tested periods does not change the values of these detected periods. Hence, the tested frequency grid density is not a trial factor effect ("look-elsewhere effect") that can change the five period values of our best M=3+6 model.

3. For all O-C data, we use Fisher-test to compare constant, linear, quadratic and cubic p(t) trends for one, two and three signal models (Table A7 ). The linear K 3 = 1 trend is the best one. This means that if the O-C data had been computed with the period 2. d 86732870 (Eq. 27), the best trend would have been the constant K 3 = 0 trend. After exploring numerous trend and signal combination alternatives in the vast free parameter space, we arrive at this simplest alternative: no trend at all in the O-C data! Although the "lookelsewhere effect" is certainly present, DCM detects this simplest trend alternative for our five signal M=3+6 model.

4. Our M=3+6 model prediction is excellent (Fig.  2 ). This indicates that the "look-elsewhere effect", or any other spurious effect, does not mislead DCM periodicity detection.

The time span of our data is "only" 236 years. Our biggest uncertainty is therefore the longest detected 219 years periodicity. It has been claimed that the Discrete Fourier Transform can sometimes detect clear signal periods slightly longer than the ∆T time span of data, "but with poor resolution" (Horne & Baliunas 1986) . The detection of periods close to ∆T depends strongly on the signal-to-noise ratio of the data. Such detections may not always succeed in our case, because we detect the 172 years period from the shortest sample of 185 years. This period is shorter than time span of this particular sample. We do not detect this "old" 172 years period from the longer samples of 226 and 236 years, but we do detect the "new" 219 years period. New additional O-C data may, or may not, confirm that this 219 years period of ours is correct.

The direct discovery of Algol H would solve the above problem for good. Eggen (1948) analysed Algol's O-C data. He arrived at an orbital period of 188.4 years for this hypothetical distant companion. Irwin (1952) estimated its orbital elements. We argue that this distant Algol H candidate may be currently found about 1.6 arc seconds away from the cEB, the eclipsing pair Algol A and Algol B.

As for other uncertainties, we can not determine the exact number of stars in this hierarchial system, but this does not prevent us from presenting an excellent 9.2 year prediction based on the first 226 years of O-C data (Fig. 2) . We admit that our longer fifty years prediction fails, because our 172 years period detected in the shortest 185 year sample is wrong (Fig. 3) . However, our turning point in this same prediction would explain the four years gap in the published O-C data around the year 1995 (Fig. A14 ). It will be interesting to see how well we can predict the future O-C data after October 2018 (Fig. 1b) .

The ephemerides of eclipsing binaries can be improved by removing linear or quadratic trends from the observed (O) minus computed (C) eclipse epochs (e.g. Kreiner et al. 2001; Kim et al. 2018) . However, even such improved ephemerides can not predict the exact epochs of future eclipses. The light-time travel effect of a third body causes strictly periodic predictable O-C changes (Irwin 1952) . The typical third and fourth body detection rates from O-C data are low, only 992/80 000 and 4/80 000, respectively (Hajdu et al. 2019) . Eclipse epoch predictions based on linear or quadratic trends, and light-time travel effects, usually fail because aperiodic trends mislead the detection of periodic signals (e.g. Bours et al. 2014; Lohr et al. 2015; Song et al. 2019) .

Considering this general background, it is unprecedented that our new Discrete Chi-square Method can detect five strictly periodic signals from 236 years of Algol's O-C data (Fig. 1a) . These tentative companion candidate orbital periods are between 1.863 and 219.0 years. One of these periods is definitely not a surprise, because our 680.4 ± 0.4 days period estimate for this weakest detected signal differs only 1.4σ from the well-known 679.85±0.04 days orbital period of Algol C. From our O-C data alone, we can not determine the exact number of companions in Algol's hierarchial system, or the stability of this system.

From the shorter 226.2 years subsample, we detect these same five above mentioned strictly periodic signals. They give an excellent prediction for the last 9.2 years of our O-C data (Fig. 2b) . Although it is impossible to detect the longest 219 year period from our shortest analysed subsample of 185 years, we can still predict the O-C data turning point epoch in the year 1995 (Fig. 3b) . This unexpected turning point event could explain the odd publication gap in the otherwise continuous modern O-C data of Algol.

We detect the linear O-C trend, which confirms that Algol's orbital period has not changed since it was discovered by Goodricke (1783) . The orbital planes of Algol C and the new other wide orbit star candidates are probably co-planar, because Algol's eclipses were observed already in Ancient Egypt (Jetsu et al. 2013; Jetsu & Porceddu 2015; Porceddu et al. 2018 ).

In the bigger picture, the predictions for complex nonlinear models rarely succeed. We give a prediction for the next decade of Algol's O-C changes after October 18th, 2018 (Fig. 1b) . These future O-C changes may prove that the abstract Discrete Chi-square Method approach works for complex non-linear models, and that Algol's data merely allowed us to check this. = 0, 45, 90, 135, 180, 225, 270 If the third body orbit is circular (e = 0), the suitable DCM model order is K 2 = 1, because the O-C curve is a pure sinusoid (Eq. 13: e = 0). If the third body orbit is not circular (e > 0), the O-C curve is not a pure sinusoid. In this case, the suitable DCM model order for these eccentric orbits is K 2 = 2 (Hoffman et al. 2006 ). Our notations for circular (e = 0) and eccentric (e > 0) orbit O-C curves are

These (O − C) e>0 and (O − C) e=0 curves have the same peak to peak amplitude A for any p, t p , e and ω combination (Eqs. 13 -17). Our notation for their difference curve is

having a peak to peak amplitude A diff . The amplitude ratio is

We also determine the phase differences

of two first minimum (t 1st.min ,t 2nd.min ) and maximum (t 1st.max ,t 2nd.max ) epochs of (O − C) diff curve. We simulate three cases of artificial O-C data (Table  A1 : Cases I, II and III). The simulated O-C values are computed for the real data time points t i from Table  A4 (n = 2224). We add 0. d 005 Gaussian random errors to these simulated O-C values. DCM period search for these simulated O-C data is performed between 8000 and 80000 days. We use the same period interval also in our DCM analysis of real data (Sects. 5.1-5.4)

In this section, we use p, A, t p , e and ω values of Case I (Table A1) . Our Fig. A1 shows all forty (O − C) e>0 and (O − C) e=0 curve pairs, as well as their (O − C) diff difference curves. We study only cases e ≤ 0.4, because our ν(t) Fourier expansion (Eq. 16) does not give the exact quantitative ν(t) values for higher eccentricities. However, our ν(t) estimates are sufficient for illustrating how eccentric orbit (O − C) e>0 curves deviate from purely sinusoidal circular orbit (O − C) e=0 curves.

The ∆A, ∆φ min and ∆φ max values for forty eccentric (O − C) e>0 curves are given above each panel of Fig.  A1 . When eccentricity e increases, the amplitude ratio ∆A increases. At the same time, the (O − C) diff curve symmetry decreases, because ∆φ min and ∆φ max values deviate more from 0.5. Both of these effects confirm that when eccentricity increases, the (O − C) e>0 curve deviates more from the pure (O − C) e=0 sinusoid. One symmetry remains: adding 180 o to ω reverses the ∆φ min and ∆φ max pair values.

In Case I, the correct one signal DCM model for simulated data has an order K 2 = 2 ≡ e > 0 (model 1,2,0 ). The number of signals (K 1 = 1) and the signal order (K 2 = 2) are both correct. The results for DCM search with this correct model are given in Table A2 . This table has the same structure as Fig. A1 . For example, the results for combination e = 0.05 and ω = 0 o are given in the upper left corner of both Table A2 and Fig. A1 .

DCM always detects the correct period p, because the ratio P 1 /p is close to unity for all forty e and ω combinations. The amplitude ratio A 1 /A is close to unity for lower eccentricities e ≤ 0.2. This ratio decreases for higher eccentricities. Yet, even in these cases the amplitude ratio is A 1 /A ≥ 0.95. The inaccuracy of our ν(t) Fourier expansion (Eq. 16) may partly explain this A 1 /A ratio decrease. DCM can certainly detect the correct simulated signal period p = 45976 d and amplitude A = 0. d 0994. Our abbreviation for this correct p period detection is "Correct-p" effect.

For eccentricities close to e = 0, the (O −C) e>0 curves for P 1 = p and P 1 = 2p periods are nearly identical. We use the abbreviation "Dp" to highlight all P 1 values for lower eccentricities e ≤ 0.1 (Table A2 ). In these cases, the spurious double period P 1 = 2p detection is possible, if the grid of tested frequencies is too sparse. The probability for detecting this spurious P 1 = 2p period would of course decrease, if our chosen simulated data error 0. d 005 were smaller. We call this spurious 2p period detection "Double-p" effect.

Here, we analyse again the same one signal simulated eccentric orbit (O − C) e>0 data of Case I, but our two signal DCM model 2,1,0 is wrong. The number of signals (K 1 = 2) and the signal order (K 2 = 1) are both wrong. In other words, we make the false assumption that the one signal eccentric orbit (O − C) e>0 curve is a sum of two circular orbit (O − C) e=0 curves. The results for this wrong model analysis are given in Table A3 . Note that this table also has the same structure as Table A2 and Fig. A1 .

The correct period P 2 = p is always detected, because the P 2 /p ratio is very close to unity for all forty e and ω combinations. The period of weaker detected signal is always P 1 ≈ P 2 /2 ≈ p/2. Furthermore, the accuracy of this approximation increases when e increases! Both of the p and p/2 periods are certainly detected at larger eccentricities e ≥ 0.2. The A 2 /A 1 amplitude ratio of these p and p/2 signals decreases for higher eccentricities. This happens at the expense of P 2 signal, because A 2 /A decreases about 10% when eccentricity increases Table A2 . Case I: Correct model results. Simulated (O − C)e>0 data signal period is p = 45976 d . Signal peak to peak amplitude is A = 0. d 0994. For different e and ω combinations, one signal DCM model1,2,0 search detects periods P1 and peak to peak amplitudes A1. Abbreviation "Dp" denotes Double-p effect cases, where spurious period P1 ∼ 2p may be detected, if tested frequency grid is too sparse. from e = 0.05 and 0.40. All these effects are also illustrated in Fig. A1 . For nearly circular orbits e ≤ 0.10, the A 2 /A 1 signal amplitude ratio is between 19 and 47. We use the abbre-viation "Hp"' to highlight the cases, where the detection of weaker spurious p/2 period signal requires a denser tested frequency grid (Table A2 ). Our abbreviation for this spurious p/2 period detection is Jetsu   Table A3 . Case I: Wrong model results. Simulated (O − C) e>0 signal period is p = 45976 d . Signal peak to peak amplitude is A = 0. d 0994. For different e and ω combinations, two signal DCM model2,1,0 search detects signals having periods P1 and P2, and peak to peak amplitudes A1 and A2. Abbreviation "Hp" highlights the Half-p effect cases, where detection of weaker P1 ∼ p/2 signal requires a denser tested frequency grid. Some ω values can eliminate the symmetry of the (O − C) e>0 curve even at these low e ≤ 0.10 eccentricities, like the e = 0.10 and ω = 45 o combination (O − C) e>0 curve that shows no "Hp" effect.

It is important to realize that every real eccentric orbit (O − C) e>0 curve can be presented as a sum of purely sinusoidal circular orbit (O − C) e=0 curve and a nearly sinusoidal (O − C) diff curve. The respective periods of these curves are p, p and ∼ p/2. All these three curves are "in-phase", and therefore the eccentric orbit (O − C) e>0 sum curve has only one minimum and one maximum.

In Case II, the simulated data contains a sum of two sinusoidal circular orbit (O − C) e=0 signals having periods p 1 = 12295 d and p 2 = 46159 d (Fig. A2) . The other parameters can be found from Table A1 (Case II). The higher amplitude p 2 signal dominates over the lower amplitude p 1 signal. These red and blue (O − C) e=0 curves, and their black (O − C) 1+2 interference curve, are shown in Fig. A2 .

In this Case II, the correct circular orbit model is DCM model 2,1,0 . This model has the correct number of signals (K 1 = 2) and the correct order (K 2 = 1). DCM detects the correct simulated P 1 = 12286 d ± 18 d and P 2 = 46122 d ± 57 d signal periods, as well as the correct amplitudes A 1 = 0. d 0170 ± 0. d 0004 and A 2 = 0. d 1019 ± 0. d 0003. In short, DCM succeeds in detecting both simulated circular orbit (O-C) e=0 signals.

Here, we analyse Case II simulated data using the wrong eccentric orbit one signal DCM model 1,2,0 . Both the number of signals (K 1 = 1) and the model order (K 2 = 2) are wrong. We detect P 1 = 46400 d ± 81 d period signal having a peak to peak amplitude A 1 = 0. d 1015 ± 0. d 0005. Since the p 2 = 46159 d period of the stronger signal dominates in the black (O − C) 1+2 interference curve of Fig. A2 , this detected P 1 period is close to, but slightly larger than, the p 2 period. Our DCM search result for P 1 is confirmed by the distance between the black (O − C) 1+2 interference curve minima, which is indeed longer than the distance between the dominating blue (O − C) 2 curve minima (Fig. A2) .

In Case III, the simulated sinusoidal (O − C) e=0 signal periods are p 1 = 12304 d and p 2 = 25274 d (Fig.  A3) . The signal amplitudes are nearly equal (Table A1 : Case II). The correct model for these simulated data is the DCM model 2,1,0 , which searches for the sum of two circular orbit (O − C) e=0 curves (K 1 = 2, K 2 = 1). DCM detects the correct P 1 = 12322 d ± 20 d and P 2 = 25259 d ± 89 d signals, as well as the correct amplitudes A 1 = 0. d 0190 ± 0. d 0004 and A 2 = 0. d 020 ± 0. d 0003. Again, DCM succeeds in detecting both simulated circular orbit (O − C) e=0 signals.

A.5. Case III: Wrong model analysis Finally, the same simulated Case III data is analysed by using the wrong eccentric DCM model 1,2,0 . In other words, we search for only one eccentric orbit (O − C) e>0 signal when the data contains two circular orbit (O − C) e=0 signals. DCM detects a signal having P 1 = 24771 d ± 34 d and A 1 = 0. d 0319 ± 0. d 0005. The simulated p 1 and p 2 signals' interference period is

where k = ±1, ±2, ... is the phase difference during p . In this particular case, k = 1 gives p = 23976 d . This black double wave (O − C) 1+2 curve is shown in Fig. A3 . DCM detects this "correct" interference signal period p , which is repeated through out the whole data. We call this spurious interference period p detection "Interference-p " effect.

The black p interference (O − C) 1+2 curve shows two minima and two maxima, because the red p 1 period and the blue p 2 period sinusoids are "off-phase" (Fig. A3) . Therefore, this black (O − C) 1+2 curve can not represent a real eccentric orbit (O − C) e>0 curve. Table A7 . All data trend (Sect. 5.1.1): One, two and three signal (K1 = 1, 2 and 3) models (M) have eccentric third body orbits (K2 = 2 ≡ e > 0). Trend orders K3 = 0, 1, 2 and 3 for p(t) are compared. Fisher-test is used to compare DCM search results between 8000 and 80000 days. Notations are "↑ " ≡ complex model above is better than left side simple model, and "←" ≡ left side simple model is better than complex model above. Parameters are F = Fisher test statistic and QF = critical level. Control file is dcm.dat. Unstable models are denoted with "Um". They have dispersing amplitudes "Ad", or intersecting frequencies "If". Some models have leaking periods "Lp" larger than ∆T = 86171 d . Best DCM model for all data is linear trend K3 = 1 order model3,2,1.

Period analysis: All original data = 1hjdAlgol.dat

M Model model1,2,1 model1,2,2 model1,2,3 model2,2,0 model2,2,1 model2,2,2 model2,2,3 model3,2,0 model3,2,1 model3,2,2 model3,2,3 dcm.dat Table A14 . All data: comparison of eccentric and circular orbit results. Cols 1-3. Eccentric orbit results (Table A8) . Cols 4-6. Circular orbit results (Table A10 ). Col 7. Connection between eccentric and circular orbit periods. Col 8. Effects are explained in Sect. 5.5. Eccentric and circular orbit periods are denoted with subscripts "e" and "c", respectively. Table A8 : Eccentric e > 0 ≡ K2 = 2 Table A10 : Circular e = 0 ≡ K2 = 1 Connection Effect M=3 Pe,5 = 79999 ± 1216 Ae,5 = 0.287 ± 0.005 M=2 Pc,7 = 120740 ± 41002 Lp Ac,7 = 0.6 ± 0.5 Ad Pe,5 ≈ 1 × Pc,7 Correct-p M=2 Pc,6 = 42422 ± 640 Ac,6 = 0.08 ± 0.01 Pe,5 ≈ 2 × Pc,6 Half-p M=3 Pe,4 = 24742 ± 141

Ae,4 = 0.029 ± 0.001 M=2 Pc,5 = 24247 ± 872 Ac,5 = 0.018 ± 0.002 Pe,4 ≈ 1 × Pc,5 Correct-p M=2 Pc,4 = 12294 ± 109

Ac,4 = 0.018 ± 0.001 Pe,4 ≈ 2 × Pc,4 Half-p M=3 Pe, 3 = 20358 ± 128 Ae,3 = 0.013 ± 0.001 M=2 Pc,3 = 10144 ± 91 Ac,3 = 0.0097 ± 0.0004 Pe,3 ≈ 2 × Pc,3 Half-p M=6 Pe,2 = 7269 ± 29

Ae,2 = 0.007 ± 0.001 M=5 Pc,2 = 7395 ± 37 Ac,2 = 0.0061 ± 0.0006 Pe,2 ≈ 1 × Pc,2 Correct-p M=6 Pe,1 = 680.4 ± 0.4

Ae,1 = 0.0064 ± 0.0007 M=5 Pc,1 = 680.7 ± 0.5 Ac,1 = 0.0057 ± 0.0009 Pe,1 ≈ 1 × Pc,1 Correct-p Table A15 . Eccentric orbit results for three samples. Cols 1-3. All data eccentric orbit results (Table A8 ). Cols 4-6. First226 y -data eccentric orbit results (Table A11 ). Cols 7-9. First185 y -data eccentric orbit results (Table A12 ). Col 10. First185 y -data period effect explained in Sect. 5.5. 

Understanding Regression Analysis

Tidal Evolution and Oscillations in Binary Stars

Applied Regression Analysis

An Introduction to the Bootstrap

Simultaneous Statistical Inference

Orbital motion Soderhjelm, S. 1974, Information Bulletin on Variable Stars

Their colours are red (z1(f1)), blue (z2(f2)), green (z3(f3)) and yellow (z4(f4)). Open diamonds denote locations of best frequencies

Best periods are at P1 = 20358 d , P2 = 24742 d and P3 = 79999 d . Otherwise as in Fig

Periodograms for residuals of M=3 model of all data (Fig. A7, blue dots): Two signal model periodograms (Table A8: M=6)

Figure A10. Model for residuals of M=3 model for all data (Fig. A7, blue dots): Two signals have periods P1 = 680. d 4 and P2 = 7290 d (Table A8, M=6 model). Otherwise as in Fig. A5

Signals in residuals of M=3 model for all data (Fig. A7, blue dots): Two signals yi,j (Eq. 20) have periods P1 = 680

Best periods are at P1 = 10144 d , P2 = 12294 d , P3 = 24247 d , P4 = 42422 d and P5 = 120740 d . Otherwise as in Fig

Signal periods are P1 = 10144 d , P2 = 12294 d , P3 = 24247 d , P4 = 42422 d and P5 = 120740 d . Otherwise as in Fig

We thank Dr. Chun-Hwey Kim for sending us the TI-DAK database O-C data of Algol. We also thank Dr. Sara Beck, Dr. Lindsay Ward, Dr. Gerard Samolyk, Dr. Stella Kafka and Dr. Nancy Morrison, who helped us in finding the O-C data of Algol from the Lichtenknecker Database of the BAV. This work has made use of NASA's Astrophysics Data System (ADS) services and the data from the European Space Agency (ESA) mission Gaia. We thank Linux Specialist Markus Minkkinen and Dr. Sebastian Porceddu from Center for Information Technology (University of Helsinki). Their computer support during the Covid-19 pandemia crisis enabled us to complete this work.