Unseen companions of variable stars (either pulsating stars or eclipsing close binary systems, or planets around pulsars) can be detected through the well-known Light Travel Time Effect (LTTE). It consists of the modulation of the variable star period induced by the orbital motion around the barycentre of the system (variable star + unseen companion(s)). As a consequence, the occurrence of the corresponding extremum light times is also modulated. The study of the LTTE has been tackled by different authors, like Woltjer (1922), De Sitter (1933), Irwin (1952, 1959), Kopal (1959, 1978) and Borkovits & Hegedüs (1996), Konacki & Maciejewski (1996, 1999), Pop (1998, 1999, 2000). The common goal of their approaches was to estimate the values of the orbital elements for the corresponding system.

Being given a timing data set and choosing the mathematical model which rules the LTTE, the values of the orbital elements are frequently estimated through the Levenberg-Marquardt least-squares minimization method. However, its application requires first/initial guesses of the values of the orbital elements. When the complexity of the investigated period variability phenomenon is relatively high, the identification of the set of initial values which leads towards the optimal fit, can be a real challenge. The solution to this problem, which was adopted by several authors, was to explore a large number of initial guesses (e.g., Hinse et al. 2012a, 2014b; Lee et al. 2014; Borkovits et al. 2015; Esmer et al. 2021).

The aim of the present paper is to review the problem of the estimation of first guesses of orbital elements relying on the parameters of an ephemeris containing a truncated Fourier series (e.g., Pop 1996) and taking into account the Fourier decomposition of the LTTE (see Kopal 1959, 1978; Borkovits & Hegedüs 1996; Konacki & Maciejewski 1996; Pop 1998, 1999, 2000). Some new formulae are derived and different practical aspects are also considered in order to account for the presence of time gaps in the analysed data as well as for the presence of the additive observational noise.

The occurrence of the extremum (minimum or maximum) light times of a variable star which has an unseen companion (stellar or planetary) is described by the well-known ephemeris (e.g., Kopal 1978)

where $t_0$ is the initial epoch, $P_p$ is the period of the observed variable star, and $n$ is the variability cycle number. The last term represents the time delay required to the light to travel across the projection of the radius vector $r_n$ of the variable star orbit onto the line of sight. In the formulae given below, we used the following notations: $a\sin i$ is the projected semi-major axis of the absolute orbit of the barycentre of the variable component around the barycentre of the hypothetical binary or triple system, $i$ is the inclination of the normal to the orbital plane with respect to the line of sight, $c$ is the speed of light, $e$ is the orbital eccentricity, $\omega$ is the periastron longitude, $T_p$ is the time of periastron passage, and $P_s$ is the orbital period.

When assuming more unseen companions, the LTTE model can also be applied, but only by neglecting the gravitational perturbations between companions (e.g., Hinse et al. 2012b, 2014a).

The constant term ${\tau }_0{F}_0(e)\sin \omega$ which appears in Eq. (18) is actually included in the observed times of minimum light. In fact, within the frame of our approach, the estimation of initial guesses of the orbital elements relies on the equivalence between the LTTE term in Eq. (1) and the truncated Fourier series appearing in Eqs. (2) and (18) [see also Hinse et al. (2012a)].

Note that the LTTE ephemeris, as derived by Kopal (1978) (Eq. (2)), can be written as Eq. (18), but with the following notation [see Eqs. (12) and (13)]:

In this section we present some methods for estimating the values of the orbital elements of the (hypothetical) absolute orbit of the variable star around the system’s barycentre; these elements are: ${\tau }_0$ or $a\sin i$, $e$, $\omega$, $P_s$, and $T_p$. The value of the orbital period can be immediately calculated from the formula $P_s=P_p/f_0$, in which, the values of $P_p$ and $f_0$ resulted from the fitting of the model (2) or (18) to the observational data.

Kopal’s (1978) approach kept the first two periodic terms from the truncated Fourier series appearing in the LTTE ephemeris (Eq. (2)) and after that, formulae for the orbital elements were derived. He also made some approximations in the formulae of the $F_k(e)$ and $G_k(e)$ coefficients:

Therefore, the orbital elements became:

Tests performed by Borkovits & Hegedüs (1996) indicated that reliable results could be obtained with these formulae for .

(i) taking into account an arbitrary number $(K)$ of amplitudes in a truncated Fourier series of type (18);

(ii) using higher accuracy estimates of the $F_k(e)$ and $G_k(e)$ coefficients, which rely on Bessel functions [Eq. (4); see also Konaki & Maciejewski (1996)] or on the equivalent power series in $e$ (Pop 1998, 2000).

The estimation of the values of $e$ and $\omega$ relies on the remark that amplitude ratios are functions depending only on the orbital elements [see also Kopal (1978) and Konaki & Maciejewski (1996, 1999)]

Let us denote $R_{k1}^{obs}={\tau }_k^{obs}/{\tau }_1^{obs}$, the respective amplitudes being those estimated from the Fourier decomposition of the ‘observed’ $O-C$ curve or by fitting one of the models (2) or (18) to the analysed timing data. As proposed by Pop (1998, 1999, 2000), the values of $e$ and $\omega$ can be estimated by minimizing the objective function

using, e.g., a simple algorithm of exhaustive search. Note that the above amplitude ratios were considered by Pop by analogy with the structural parameters of the light curves of pulsating stars, defined on the basis of their Fourier amplitudes (e.g., Simon & Lee 1981).

We note that:

(a) due to the fact that $A_k$ $(k=1,2,\mathrm{\dots },K)$ are periodic functions with a period of $\pi$, the ambiguity between the values $\omega$ and $\omega +\pi$ has to be removed;

(b) according to our numerical tests, the objective function $\delta$ [see Eq. (35)] is much more sensitive to the variation of $e$ compared to the one with respect to $\omega$ (see also Figs. 4 and 5 below). Therefore, we proposed another approach for estimating the value of $\omega$ [see Eq. (39) below].

Regarding the amplitude ratios from Eq. (34), $R_{21}^{calc}$, $R_{31}^{calc}$, $R_{41}^{calc}$, and the objective function, $\delta$, from Eq. (35), we add the following remarks.

In Fig. 1 we plotted the ratios $R_{21}$, $R_{31}$, and $R_{41}$ (in descending order) for all values of $e$ and $\omega$ [$e\in (0,1)$, $\omega \in [0,\pi )$]. One can observe that $R_{k1}(e,\omega )=R_{k1}(e,\pi -\omega )$, $k=2,3,\mathrm{\dots },K$, therefore the graphical representations of $R_{k1}$ are symmetric with respect to the plane $\omega =\pi /2$. Konacki & Maciejewski (1996, 1999) had a similar plot approach, for two of these quantities, stating that these ratios are, ”for small eccentricities (…), with good accuracy, functions of $e$ only”, remark similar to the one above in (b). We can add that Fig. 1 shows the surfaces of $R_{21}$, $R_{31}$, and $R_{41}$ related to the LTTE induced by a purely Keplerian orbital motion.

Fig. 2 displays the locus of the points of coordinates $(R_{21},R_{31},R_{41})$, for all $e$ and $\omega$, corresponding to a purely Keplerian LTTE. It can be used in differentiating the LTTE mechanism when compared to other arbitrary periodic phenomena and in determining how close a hypothetical movement is to a Keplerian one [see also Fig. 1 of Pop (1999)].

In order to illustrate our approach, we generated a data set of LTTE timing data using the orbital elements for the hypothetical companion of the W UMa type binary system V524 Monocerotis (He et al. 2012; Bulut & Aşkın 2022), superposed on a Gaussian noise with a relatively low standard deviation (see Fig. 3). Due to the fact that the number of the available observed data is relatively small (66), and the orbital cycle coverage is not rich, we increased the number of artificially generated data (200), keeping the same time base. In addition, the data sampling is also irregular and contains random length gaps.

In Fig. 4 we plotted the corresponding objective function, $\delta$, as given in Eq. (35), for the amplitude ratios $R_{k1}^{obs}$ for the simulated data set.

Fig. 5 shows the contour plot of the isolines of the objective function given in Eq. (35) for the simulated LTTE data. As we already mentioned above, for small values of $e$, the objective function $\delta$ is less sensitive to the variation of $\omega$, while for larger values of $e$, $\delta$ becomes more and more sensitive with respect to $\omega$ (see also Fig. 4).

As previously indicated, in order to consider the investigation of a period variability phenomenon possibly caused by a LTTE, we have to fit a truncated Fourier series to the available timing data. Therefore, it is useful to take into account the run of the respective $O-C$ curve, and also to analyse its amplitude spectrum in order to obtain further information. The following situations may be encountered:

(i) Most frequently the timing data displays time gaps which limits the order of the fitted Fourier series. Thus, according to (a) the local weather conditions experienced by the respective observers, and, sometimes to (b) the historical context (see, e.g., the case of the two World Wars), the available timing data can contain more or less wide time gaps. In order to avoid the oscillation of the fitted curve inside these gaps, a simple condition has to be fulfilled [similar to that proposed by Pop et al. (2004)]: half of the minimum dimensionless period $P_{0min}$ $(P_0=1/f_0=P_s/P_p)$ corresponding to the frequency of the highest harmonic which will be taken into account $K_{max}{f}_0$ has to be equal to the width of the largest time gap measured in cycle numbers $(\mathrm{\Delta }{n}_{max})$, i.e.,

which becomes

(ii) According to the observational technique and its location (ground or space based), the analysed timing data may contain a more or less amount of additive observational noise.

(iii) In the case of the Fourier expansion of the elliptical motion, the large orbital eccentricity values yield large amplitudes of the higher order harmonics (see Konacki & Maciejewski 1996). On the other hand, according to the level of the observational noise, only (very) low order harmonics can be observed in the amplitude spectrum of the detrended $O-C$ residuals. In case of lower eccentricity orbits, the second harmonic $(3f_0)$ and even the first one $(2f_0)$ may be significantly contaminated by the contribution of the observational noise. As a consequence, by taking too many harmonics, the fitted model will tend to describe the fluctuations caused by the noise. Clearly, in such a situation, the attempt to estimate the values of the orbital elements might be seriously affected.

An elegant solution to this problem may be the use of the Lanczos smoothing factors, already applied by Meylan & Burki (1986) in the case of the radial velocity curves of some pulsating stars. The Lanczos factors are in fact weights which multiply the estimated amplitudes of a truncated Fourier series [see, e.g., Arfken (1970), Hamming (1986)]. For our model, in Eq. (18) they are given by the formula

As it is well-known, these factors have values less than 1 and, for $k=K$, ${\sigma }_K(K)=0$, which means that the highest harmonic taken into account is completely eliminated. Thus, by reducing the harmonics amplitudes as well as their number, the Lanczos factors perform a smoothing of the Fourier fit. Therefore, a predictable consequence might be an underestimation of the orbital eccentricity value.