© ESO 2019

1. Introduction

In a binary system consisting either of two stars, or a star and a planet, tidal forces act on the non-mass point component(s). These tidal forces lead to several phenomena, including a decrease of the semi-major axis (tidal decay), a decrease in the orbital eccentricity (circularization), a decrease of the rotational rate (synchronization), and a precession of the orbit in eccentric cases, in other words, the semi-major axis of the orbit rotates. This latter phenomenon is called apsidal motion or perihelium or periastron precession. Its timescale is much shorter (10–1000 yr) than the timescale of the other listed effects (10–1000 Myrs). The periastron point of an exoplanet’s orbit is at an angle ω from the tangential line of the sky (Fig. 1). We point out here that if apsidal motion is present then this angle is time-variable, resulting in transit timing variations (TTVs) as well as radial velocity (RV) variations (see also Ferrero et al. 2013; Schmitt et al. 2016; Rauw et al. 2016) and its rate is governed by the Love numbers.

The fluid Love numbers k_i,Love (i = 2, 3, 4...) of celestial bodies measure the reaction of a body to perturbing forces. Assuming the body is in hydrostatic equilibrium, they are a direct, but complicated function of the internal radial density distribution. In the simplest two-layer, core-mantle models with polytropic equation of state, k_2,Love ~ f(M_core∕M_total) where M_core and M_total are the core and total masses of the planet, respectively (Becker & Batygin 2013). The full derivation of the Love numbers depend on the exact internal density-profile and is given in, for example, Padovan et al. (2018) Kellermann et al. (2018), Kramm et al. (2011), Kopal (1959) and Love (1911). The goal of this paper is to develop a model that includes the apsidal motion effect in the study of RV and TTV data and to fit this model to the available WASP-18 data, to estimate the planetary fluid Love number. The presence of apsidal motion requires an eccentric orbit, and the non-circular nature of the WASP-18Ab orbit is discussed in Sect. 6.2.

2. Causes of apsidal motion

The argument of periastron changes because of general relativistic effects at a rate of (Einstein 1915): (1)

and as a result of the non-sphericity and rotation of the components (star and planet), which produce rotational potential changes and tidal interaction (Sterne 1939): (2)

Here N denotes the “Newtonian”-term, n is the mean motion, 2k₂ is the 2nd order fluid Love number, a is the semi-major axis, e is the eccentricity, R_star and R_planet are the radii of the star and planet, respectively and M_star and M_planet are the masses of the star and planet, respectively. P_orb, P_rot,star, and P_rot,planet are the orbital period and the rotational periods of the star and planet, respectively.

Equation (2) is valid if the obliquity of the star is close to zero – Rossiter-McLaughlin measurements can provide this information.

The planet’s contribution to apsidal motion is about 100 times bigger than that of the star (Ragozzine & Wolf 2009). On the left-hand side the apsidal motion is measurable in principle from TTVs or from RVs as we show hereafter.

The eccentricity is known from RV analysis and/or from occultation observations. RV analysis can even provide the mass ratio for SB1 systems if we have a good estimate for the stellar mass from isochrone fits or from asteroseismology and the secondary’s mass is very small with respect to that of the primary. The period and thus the mean motion, the fractional radii R_i ∕a are known from transit light curve analysis. The rotational period of the star can be estimated accurately enough from its true radius – provided by isochronefitting, asterosoismology or by analysis of its emission and distance by for example, Gaia – and this, combined with the V_esin i value of the star, yields the stellar rotational period. Out-of-transit light curve modulations can also be utilized to estimate P_rot although this technique is limited to stars with significant surface features (spots and/or plages). Then, k₂ of the star can be obtained from theoretical calculations (we use the tables of Claret 2004) and observational results show that such calculations agree reasonably well with the observed values and they vary in a very limited range (Torres et al. 2010). Thus, if we assume a rotational rate for the exoplanet, for example, that it is synchronous or in a 3:2 resonance (like Mercury) with the orbital period, the Love number of the planet can be determined. We show that the results are in fact not very sensitive to the planetary rotational rate. We note that the argument of periastron can change because of a perturbing third body, as well as by magnetic interaction between the star and planet.

We also note that the Love-number is defined by for example, Becker & Batygin (2013) as (3)

while Sterne (1939) used the so-called apsidal motion constant: (4)

where η₂ is the solution of Radau’s equation for j = 2 at the surface of the star or the planet. The similar notations in previous works are unfortunate, and therefore we use the word Love in the index where we mean the Love-number and we use simple k₂ for the apsidal motion constant.

Fig. 1 Illustration of the apsidal motion. At time t0 the apsis (semi-major axis of the orbit) is oriented at an angle of ω0 (right). At time t it moved to the direction characterized by (left). Traditionally, the argument of periastron is measured from East to the direction of the observer, however, the radial velocity (and the TTV) equations are invariant if we measure it from West. This causes no difficulties in the measurements of , and the East–West discrepancy can be resolved via time-series done via interferometry, astrometry or direct imaging.

3. TTV and radial velocity variations under apsidal motion

We consider the apsidal motion as a linear process in time (5)

and in time (its second derivative is taken to be zero) where t is the time.

An OXY Z coordinate system is defined in such way that its origin O is in the common center of the mass (CMC), X and Y are oriented to east and north in the tangential plane of the sky, and axis Z connects the observer and the star and it is oriented from the origin away from the observer. The Z-coordinate of the star can be calculated step-by-step as follows: (6)

where P is the anomalistic period (time between two consecutive periapsis passages). This period is constant because it appears in Kepler’s third law and we do not consider mass loss from the star and the planet, nor changes in the semi-major axis due to a perturber or tidal forces. (Change in a has much longer timescale than changes in ω owing to tidal forces.) We note that transit observations give the mid-transit time T_tr,N at cycle N and that is connected to the anomalistic period as (up to first order in eccentricity, cf. Gimenez & Garcia-Pelayo 1983, Eq. (19) and Csizmadia 2009): (7)

and for the occultation (8)

One can see that the transit period P_tr = T_tr,N+1 − T_tr,N will show oscillations with the apsidal motion period (denoted by U usually, ) and with amplitude of ~eP∕π in the first-order approximation.

The mean anomaly M of the star is (9)

where t₀ is fixed and t_p is the periastron passage time at the epoch.

Let t₀ be the epoch of a well observed transit (in practice, we use a value determined from several transits observed in a short time-window which can be more precise than just a single observation). At that time, the transit occurs when the true anomaly is (10)

ω₀ is the argument of the periastron at the transit epoch. ϑ₀ is a small correction due to the fact that the mid-transit time occurs not at conjuction but at the smallest sky-projected distance of the star and the planet (see Csizmadia 2018; Gimenez & Garcia-Pelayo 1983; Martynov 1973; Kopal 1959): (11)

where the upper and lower signs are valid for the transit and the occultation, respectively, and i denotes the orbital inclination. The eccentric and mean anomalies of the point where the mid-transit occurred are (12) (13)

M and M₀ are actual meananomalies and its value at the epoch. Then one gets from Eq. (9) that (14)

The eccentric and true anomalies are (15) (16)

where E and v are the eccentric and true anomalies. If a₁ is the semi-major axis of the star around the CMC and r₁ its actual distance from CMC, then (17)

Concerning the Z-coordinate: (18)

where d is the distance to the star from Sun and V_γ is the systematic radial velocity of the system to the Sun.

The radial velocity is the time-derivative of the Z-coordinate plus an apparent tidal component, V_tide, which is non-negligible in WASP-18 system: (19)

With and V_tide = 0 we get back the usual radial velocity equation free of apsidal motion and ellipsoidal components. We used the usual notation for the RV half-amplitude (20)

while the second term in Eq. (19) comes from the time-derivative of ω. The component denoted by V_tide stems from the fact that the star is ellipsoidally distorted because of the tidal potential of the planet, therefore we see smaller projected stellar area at conjuctions and larger at quadratures which causes line shape distortions. This is reflected in the observed radial velocities and must be taken into account (Kopal 1959; Arras et al. 2012). A detailed formulation of this effect is given in Sect. 6.2 by Eqs. (22)–(28). The model we fit to the observations in Sect. 5 is described by Eq. (19).

4. Sensitivity analysis

4.1. TTVs

Transit times can presently be determined to a precision of around 20 s with space-based measurements (Csizmadia 2010). This precision is at least a factor of three worse for ground-based observations of shallow transits. Substitutingthe known system values of WASP-18Ab (see Table 1) and assuming a Love number of WASP-18Ab close to Jupiter’s k_2,Love = 0.34 then we get that the periastron precession rate should be ° per day¹. By virtue of Eq. (7) we determine that the corresponding TTV-amplitude is about eP∕π ~ 3.7 min and the corresponding sinusoidal TTV-curve has a period of ~266 yr. Hence, the annual rate of change in the first ten years is about 0.01 s yr⁻¹. Thus, it is not surprising that the TTV method did not detect it yet. It also predicts that the observed transit period should be highly stable for a long time. This is confirmed by Wilkins et al. (2017) who found that the transit period of WASP-18Ab is stable to less than 6.6 × 10⁻⁷ days since the discovery (cf. the substitution to Eq. (7) which yields much smaller variation in the period than their upper limit). Since TTVs likely do not detect the apsidal motion in this system because the eccentricity is small, the timing precision is low and the observational window is not long enough, we investigated the opportunities provided by RVs.

4.2. RVs

Due to the apsidal motion the periastron is shifting and a phase shift will occur which is observable. This shift causes a velocity difference between the actually observed and the one predicted with constant ω, that is, with . We illustrate this difference in Figs. 2 and 3 where we show the difference in radial velocities between zero and non-zero -values for different time-periods. Our estimate shows that for WASP-18Ab this should be on the order of 100 m s⁻¹ difference at its maximum after 2000 days of observations which is about eight times bigger than the typical observational error bars of one RV-point (Fig. 3). This shows that RV variations are much more sensitive to the apsidal motion than TTVs because they are more accurate.

Fig. 2 Evolution of the maximum deviation of the observable and the predicted RV-measurements. This maximum value was selected as the maximum deviation over an orbital cycle. The observable RV curve was calculated with the parameters given in Table 1. The big change is due to the fact that after a while we assume completely opposite phase for the RV-fit if we neglect the apsidal motion.

5. Data, model fit and results

We took the 123 published RV observations of Hellier et al. (2009), Triaud et al. (2010), Albrecht et al. (2012) and Knutson et al. (2014). We excluded the points between orbital phases − 0.1 to 0.1, meaning theones obtained during primary transit, because we did not fit the Rossiter–Mclaughlin effect. Therefore, we used 54 RV points for the subsequent analysis.

These RV observations span 1849 days (=5.06) yr – long enough to suggest that apsidal motion should be observable in this dataset (Figs. 2 and 3). We carried out four fits with different approaches:

• M1:

  Model I, we fit as a freeparameter and we use the RV data points only without using the transit and occultation timing data. We calculated the Love number from the derived apsidal motion rate.

• M2:

  Model II, we fit as a free parameter and we use the RV data points and the transit and occultation timing data simultaneously. This approach allows us to constrain better the argument of periastron and . As in M1, we calculated the Love number from the measured value of .

• M3:

  Model III, we considered the apsidal motion constant (half the Love number) of the planet k_2,planet as a free parameter and we fitted Eqs. (2) and (19) together to all data.

• M4:

  Model IV, where we fixed and k_2,star = k_2,planet = 0 for comparison purposes – in this model there is no apsidal motion present.

All timings in this paper are in barycentric Julian date, in the barycentric dynamical time system (BJD_TDB).

Fig. 3 Zoom to the first five years of Fig. 2. The typical averaged error bar of the available observations (see Sect. 5) is noted by the thick vertical line. The curve shows the maximum deviation during a revolution from an RV-curve without apsidal motion over time.

5.1. Fit of RV data only (M1)

In Model 1, we fitted only the RV data. The free parameters were: V_γ, K, , – to avoid correlations between parameters e and ω₀ (Albrecht et al. 2011) –, , P and four RV-offset values between the different instruments. t₀ was fixed at the epoch given by Hellier et al. (2009). We fitted the model outlined in Sect. 3, Eq. (19) by using a Genetic Algorithm-based Harmony Search optimizer (HS, Csizmadia 2018) with 1000 individuals in the population. We tried RV-jitter values of 0 and 10 m s⁻¹ and by a bisection method we adjusted the RV-jitter – fixed during the HS-run – until we reached , N_obs being the number of RV-observations.

The results of the HS were refined by running ten chains of MCMC, with each chain consisting of one hundred million steps. A thinning-factor of 1000 was applied. The marginalized likelihoods were determined and the peak of that distribution defined the final solution. The uncertainties were estimated by calculating the 68% confidence levels from the marginalized likelihoods.

The correlation plots of the parameters are also produced from the MCMC chains, and we calculated the Pearson correlation coefficients for each parameter pair, as well as the Gelman–Rubin statistic to check possible correlation and convergence problems. The Gelman–Rubin convergence test, denoted by R, showed that R < 1.1 for all parameters, indicating that convergence was reached (Sect. 5.3).

To get the Love-number, the rotational period of the star can be estimated from the known R_star = 1.36 ± 0.06 R_⊙ (Triaud et al. 2010) and from the V_e sin i_* = 11.5 km s⁻¹ (Hellier et al. 2009), and assuming that stellar rotational axis has i_* = 90° (suggested by Rossiter–Maclaughlin measurements performed by Triaud et al. 2010 and Albrecht et al. 2012 of the sky-projected obliquity of the system which are consistent with zero), we get P_orb∕P_rot = 0.17184 ± 0.01345.

The results of the fit are given in Table 2. The apsidal motion rate – comprising relativistic and Newtonian terms – is degree per day showing a tentative 1.6σ detection of apsidal motion based purely on the RV data. For the most likely synchronous planetary rotation case we found k_2,planet = 0.48 ± 0.37 and thus k_2,Love = 0.96 ± 0.74. The uncertainties on the system parameters were propagated when we calculated the Love-number from Eq. (2) for different rotational statuses of the planet, listed in Table 2. That is why the uncertainty range of the Love-number is so high.

5.2. Fit of RV+TTV data with (M2)

In Model 2 we combined the RV-data and the available transit and occultation timing data. Therefore we used the inclination taken from Southworth et al. (2009) to get the value of ϑ in Eqs. (7)–(8) for the transit and occultation times. Its error bar was propagated to the fit by drawing the inclination value from a Gaussian distribution.

We took the four mid-occultation times observed by Spitzer and nine primary transit times observed from the ground (WASP, TRAPPIST, see Triaud et al. 2010; Nymeyer et al. 2011; Maxted et al. 2013) but excluded the less precise amateur measurements published on the ETD page². We also did not use the transit timing of McDonald & Kerins (2018) which is based on sparsely sampled HIPPARCOS data spread over more than three years. Their error bar on the mid-transit time is about 15 min while the expected O–C value of the whole effect we search for is about 3.7 min. We iteratively solve Eqs. (7), (8), and (11) to predict the time of each transit and occultation. By this usage of the known timing of primary and secondary transit events we can constrain our RV-model further, especially the value of the argument of periastron. TTV and RV data were fitted simultaneously. The results are shown in Table 2. We found degree per day which is a circa 3.2σ detection of the apsidal motion. In other words, the inclusion of transit and occultation timing data increased the significance level of the detection because the argument of periastron and is better constrained by the timing data. This yielded k_2,planet = 0.32 ± 0.16, k_2,Love = 0.64 ± 0.32.

5.3. Fitting the RV+TTV data with k_2,planet

Another possibility to fit the joint RV+TTV data set is to use k_2,planet as a free parameter and to calculate the relativistic and the Newtonian apsidal motion rates via Eqs. (1) and (2) for the fit. The advantage of this formalism is to use the posterior distribution for error estimation.

The fitted parameters were in this case: anomalistic period P, V_γ, K, k_2,planet, four RV offsets to take account of the instrumental offsets between the five instruments used, and , . The epoch of transit was taken from Hellier et al. (2009) and fixed. Any uncertainty or error of it was corrected by ω₀. We minimized the χ² given by (21)

where the last eight terms are priors on the semi-major axis a, the ratio of the orbital period and the stellar rotation period P_orb∕P_rot, the scaled semi-major axis R_star∕a, the radius ratio R_planet∕R_star, effective temperature of the star T_eff, the stellar k_2,star, the stellar mass M_star and observed transit period. The model values (index m) were drawn with gaussian distribution centered at their literature values with their one sigma standard deviation (Table 1). We note that they are priors and not fitted parameters, so the model has only ten free parameters and eight priors. The priors are used this way to propagate the error bars of the stellar parameters, too.

The lastterm in the priors expresses that the observed transit period of Hellier et al. (2009) is related to the transit period derived from Eq. (7), neglecting the tiny terms. This decreases slightly the correlation between k₂ and anomalistic period.

One can notice that most of the correlations (Fig. 4) between the parameters are negligible except three pairs: k₂ –period, k₂ –e cos ω and e cos ω–period. This can be understand easily: in the term cos(v + ω) of Eq. (19) one can approximate v = nt + 2e sin nt + ⋯ and therefore we will have and therefore the term causes the correlation between the anomalistic period and k_2,planet. This correlation can be broken down by more eccentric orbits, more observations and/or longer MCMC-chains. The correlation does not mean that we can measure only the ratio of k₂ ∕P because the anomalistic period P appears in other terms without a relation to k₂. As it is clear from Fig. 4, the best solutions are concentrated onto a small part of the period-k₂ diagram and therefore the correlation just increases the uncertainty of the finally derived parameters but it still allows us to constrain the Love number of the planet.

Because of the five years of the data coverage and of the phase-shifts caused by apsidal motion it is difficult to visualize the fits. Therefore, we decided to show the residuals of the fits of the RV, transit and occultation data. These can be seen in Figs. 5–7 and the resulting parameters of the fits are listed in Table 2. Notice that the found RV-jitter is in perfect agreement with Albrecht et al. (2012). Considering the transit and occultation data in Figs. 6 and 7, one can speculate that such an improvement in TTVs might be in accordance with the prediction of Iorio (2011), who suggest that a shift of the occultation will be observable within a decade for ultrashort-period planets.

We also carried out a fit without any apsidal motion and without apparent tidal RV-term (M4) and the residuals of that fit are compared to the ones of M3 in Figs. 5–7 and in Table A.1. As one can see, the fit with apsidal motion+apparent tidal term provides a better residual curve than without it, preferring M3 model over M4.

Fig. 4 Correlation plots and probability histograms for WASP-18AB RV+TTV fit. Red and orange areas denote the 1σ region, different gray areas denote the 2σ and 3σ regions. Numbers over the panels are the values of the Pearson-correlation coefficients between the parameters. The solid and dashed lines show the most probable solutions and the 1σ limits in thehistograms. Only a subset of the parameters are shown.

6. Discussion

Models 2 and 3 give highly consistent values for the Love-number and for the apsidal motion rate, in other words there is no significant difference between fitting the apsidal motion constant and fitting the apsidal motion itself directly. However, the error propagation in Model 3 is more robust, and so this is our preferred solution. The following discussion is therefore based on the results of Model 3 (Sect. 5.3, Table 2 and Figs. 4–7).

6.1. Love number of WASP-18Ab

Assuming a synchronous orbit for the planet P_rot,planet = P_orb, we derive an apsidal motion constant and a Love number for WASP-18Ab: and . Its significance is ~ 3.3σ. The apsidal motion rate is about 0.0087 ± 0.0033° per day (Table 2), so its significance means a strong tentative detection at 2.6σ level. This means an apsidal motion period yr. When we decreased the error bars of the stellar parameters, the error bar on the apsidal motion rate decreased as well – the biggest contribution to the error bar stems from the scaled semi-major axis. That clearly shows that the uncertainties of the known stellar and system parameters dominates the error bars in the determination of k_2,planet. Any future improvement in these parameters is encouraged. We repeated the fit described in Sect. 5.3 with other rotational rates of the planet, too. If the planet rotational-orbital periods are in 3:2 resonance, then the apsidal motion constant is . If the planet rotates very slowly () then it is . In the unlikely event that the planet rotates three times faster than the orbital period, . (The escape velocity from the planet governs the maximum rotation rate and this means that WASP-18Ab can rotate a maximum of approximately 24 times faster than its orbital revolution. However, the tidal forces synchronize the orbit quite fast.)

Thus, although some trend can be seen, the planet’s rotational speed does not influence the k₂ -estimation significantly in the realistic cases and within the present level of accuracy. This is because in Eq. (2) the second, tidal term is about fifteen times bigger in case of synchronous rotation than the first, rotational term and therefore it dominates the expression. Such moderate values of k₂ are suggestive of a more homogeneous density distribution of a massive hot Jupiter, and the uncertainty range allows us to say that WASP-18Ab has a similar Love number as Jupiter in our own solar system (Ni 2018). It is beyond the scope of this paper to establish which giant planet interior models are compatible with the determined Love number.

Fig. 5 RV-residuals of the model fit. x-axis: the index of the RV points used for the fit, in order of increasing time. y-axis: the residual of the fit in m s−1. The dashed black curve represents the residuals of the fit without any apsidal motion because we forced (including all relativistic and classical tidal terms), as well as Vtide = 0.0 m s−1. The solid red curve represents the residuals of the model with apsidal motion (see Eq. (19)). A significant improvement in the residuals can be seen. Table A.1 helps to identify the observational point via its index.

Fig. 6 Residuals of the model fit on the transit observations. We note that the transit observations were taken from the compilation of Maxted et al. (2013). x-axis: the index of the transit observations used for the fit, in order of increasing time (not the cycle number!). y-axis: the residual of the fit in minutes. The open black points represent the residuals of the fit without any apsidal motion because we forced (including all relativistic and classical tidal terms), as well as Vtide = 0.0 m s−1. The red filled circles represent the residuals of the model with apsidal motion (see Eq. (19)). For sake of clarity, the red symbols were shifted horizontally by 0.15 units. Table A.1 helps to identify the observational point via its index. No improvement can be seen in the transit data with the more complicated apsidal motion model which is contrary to the case of RVs (see Fig. 5.) where we could see significant improvement. This is in accordance with the expectation of Sect. 3.

Fig. 7 Residuals of the model fit on the occultation observations. Notice that the occultation observations were taken from the compilation of Maxted et al. (2013), too. x-axis: the index of the occultation observations used for the fit, in order of increasing time (not the cycle number!). y-axis: the residual of the fit in minutes. The meaning of symbols is the same as in Fig. 6. Table A.1 helps to identify the observational point via its index. A small improvement can be seen in the last occultation measurement with the complicated apsidal motion model relative to the simple model.

6.2. Various views about the eccentricity of WASP-18Ab

WASP-18Ab was discovered by Hellier et al. (2009) who reported an eccentricity of e = 0.0093 ± 0.0029 based on only nine RV points obtained by CORALIE. According to Eq. (7) of Zakamska et al. (2011), who estimated the expected precision as a function of number of data points, K-amplitude and observational error bars, the uncertainty on the eccentricity in Hellier et al. (2009) is underestimated by a factor of two. However, Zakamska et al. (2011) did not study how the joint fit with transit light curves can increase the precision in eccentricity.

Arras et al. (2012) called attention to the fact that the ellipsoidal shape distortion of the star, caused by the tidal interaction with the planet, introduces a distorted RV signal. They propose that one can observe a small apparent eccentricity with an argument of periastron of 270° even if the orbit is really circular, and they predicted an RV distortion of ~30 m s⁻¹ for WASP-18Ab. Such spurious, mimicked eccentricities were also mentioned by Kopal (1959, Sect. V.4). In contrast to Arras et al. (2012) Eq. (V.1-22) of Kopal (1959) yields only ~2 m s⁻¹ of such RV distortion.

We find that Arras et al. (2012) account only for the effect of tides; they neglect the effect of stellar rotation, and therefore overestimate the change of shape of the star. When we set synchronous rotation for the star, then Kopal’s formula gave ~24 m s⁻¹ of tidally induced RV variations, close to the 30 m s⁻¹ predicted by Arras et al. (2012). But synchronous rotation is not the case for the host star WASP-18A. Therefore, one can propose to revise our views on small exoplanetary eccentricities when it was suspected that they are just apparent eccentricities and stem just from tidal effects. To help such analyses, we give the formula (only the dominating term), based on Kopal (1959) but using modern notation (valid for linear limb darkening with coefficient u₁, which is a good enough approximation here) which should be added to the radial velocity model (cf. Eq. (19)): (22)

with (23) (24) (25)

and (26) (27)

μ and K stand for micrometer and Kelvin, respectively, and lastly (28)

q is the planet-to-star mass ratio, λ is the effective wavelength of the observations (we set 0.6 μm), τ₀ describes the gravity darkening effect.

Nymeyer et al. (2011) obtain occultation light curves at 3.6 and 5.8 μm and from the phase shift and duration of the transit and the occultation they obtained e = 0.0091 ± 0.0012, a 7σ significant eccentricity. Maxted et al. (2013) analyzed the 3.6 and 4.5 μm occultation light curves of the system and they reported e = 0.003 ± 0.004 which is compatible with both a circular and the aforementioned eccentric orbits.

Knutson et al. (2014) obtain six new RV points and determined an eccentricity of 0.0068 ± 0.0027 and an argument of periastron of ω = 261.1 ± 7.4 degree. One can speculate that the periastron precession contributes to the error bar in the case of their fixed ω. The significance of their eccentricity is 2.5σ.

Our result: e = 0.0085 ± 0.0020 is a very strong detection of the eccentricity because we fitted a time-variable ω. Actually, it differs by only 0.2σ from the eccentricity Hellier et al.’s (2009) and it does by only 0.5σ and 0.3σ from Knutson et al. (2014) and Nymeyer et al. (2011), respectively³. We consider these differences negligible.

Our ω₀ = 257.27° ± 2.13° rules out the exact 270° argument of the periastron at the epoch which was fixed at T₀ = 2 454 221.48163 (recall that the epoch is taken from Hellier et al. 2009). We must note that the epoch was fixed in our solution and its uncertainty is maybe compensated by ω₀. We suggest obtaining new, very precise, high quality photometric and radial velocity measurements of WASP-18Ab’s primary transit(s) to carry out a new joint fit (cf. Southworth et al. 2009, 2010).

We propose that the eccentricity is real and it is not caused by tidally induced RV signals. This is supported by substitution to Kopal (1959)’s equation and even more by the Spitzer measurements of Nymeyer et al. (2011). We note that the observed transit and occultation lengths of Nymeyer et al. (2011) and Hellier et al. (2009) differ by 303 ± 84 s which is a 3.6σ difference. Such a difference in the durations is expected for an eccentric configuration.

6.3. A visual companion to WASP-18A

We searched the second data release (DR2; Gaia Collaboration 2018) from the Gaia mission (Gaia Collaboration 2016) for close companions to WASP-18A. This revealed an object 11.75 magnitudes fainter than WASP-18A in the Gaia G passband, at a separation of just under 30′′. The two stars have proper motions and parallaxes that are consistent with each other to within 1σ (Table 3). We therefore conclude that they are likely to be a bound pair. Based on the magnitude difference, and the absolute magnitudes tabulated by Pecaut & Mamajek (2013)⁴, we suggest that the companion is a late-M type star (approximately M6.5V). Combining the visual separation with the distance to WASP-18A (123.92 ± 0.37 pc), we find that the companion is physically separated from WASP-18A by 3519 ± 9 AU.

Assuming a mass of 0.1 M_⊙ for the companion star, we are able to estimate the radial velocity induced by the companion on WASP-18A. When we assume a circular, coplanar orbit with WASP-18Ab then the over five years of observations we can see a maximum velocity change of (29)

where Δt = 5.06 yr is the time window between the first and final RV-observations. A substitution yields V_max ~ 0.1 m s⁻¹ radial velocity drift which is negligible in this study. Via Kepler’s third law, the orbital period of the stellar companion around the primary star can be estimated tobe ~180 000 yr.

We also estimated the impact of the third body on the apsidal motion rate. To do so we used Eq. (12) of Borkovits et al. (2011), assuming a coplanar, circular outer orbit we have that (30)

which is also negligible.

Because of this companion we changed the planet name from the earlier-used WASP-18b to WASP-18Ab because it orbits around the brighter primary star, and the faint, newly-discovered companion star is named here WASP-18B.

7. Conclusions

In this study we show that there is tentative evidence for the presence of apsidal motion in the WASP-18Ab system. We used archival data spanning more than five years to demonstrate this.

When the origin of this apsidal motion is a combination of general relativistic effects and tidal interaction between the host star and theplanet, then it follows that the apsidal motion constant and the fluid Love number of WASP-18Ab are and , respectively. This result assumes a synchronous rotation for the planet (but not for the star) and the result is thus model and assumption dependent because the stellar fluid Love number is taken from theoretical calculations. We show that the Love number is largely insensitive to the actual rotational period of the planet because the tidal contribution to the apsidal motion is about fifteen times larger than that of the rotational modulation.

Such a Love number is compatible with the formerly known Love number of Jupiter in our own solar system (k_{2,Love-Jupiter} ~ 0.34) because of the large error bars, and compatible with k_{2,Love-Jupiter} = 0.53 proposed byNi (2018) based on Juno’s recent measurements. By substitution one can see that about 13% of the total apsidal motion stems from general relativity and 87% comes from the tidal interaction. It is not very likely that the observed and eccentricity arise from an apparent effect, such as that proposed by Arras et al. (2012). If we do not observe any periastron precession, then we must assume that the star and the planet both have zero Love number, or that the orbit is circular against so many findings, or that a third body neutralizes the effect of the apsidal motion, which is particularly unlikely. New RV observationscan help to solve this issue.

An unseen third body may also be the reason for the observed apsidal motion, but the lack of TTVs – which is not contradictory to the observed apsidal motion at the found range – and the missing RV trend excludes the presence of a third body with a significant impact on the results. We found a stellar companion to the transit host star, the mutual true separation at the epoch of Gaia observations is about 3519 AU. We estimated its effect on the radial velocity and periastron precession rate and it was found to be negligible. However, this companion may contribute to the excitation of the orbit causing the eccentricity and preventing circularization (e.g. Borkovits et al. 2011).

Thus, the most plausible result is that we can see signs of the apsidal motion generated by tidal interaction and by general relativity in the five-year long RV curve of WASP-18Ab – and the observed Love number also is in the plausible range. We feel, however, that new observations are needed to confirm or reject this result as we are observing a very small signal in a limited dataset.

Acknowledgements

H.H. and Sz.Cs. thank DFG Research Unit 2440: “Matter Under Planetary Interior Conditions: High Pressure, Planetary, and Plasma Physics” for support. Sz.Cs. acknowledges support by DFG grants RA 714/14-1 within the DFG Schwerpunkt SPP 1992, Exploring the Diversity of Extrasolar Planets. Sz.Cs. thanks the Hungarian National Research, Development and Innovation Office, for the NKFIH-OTKA K113117 and NKFI-KH-130372 grants. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. We thank an anonymous referee for his/her useful comments.

Appendix A: Additional table

References

All Tables

All Figures

Fig. 1 Illustration of the apsidal motion. At time t0 the apsis (semi-major axis of the orbit) is oriented at an angle of ω0 (right). At time t it moved to the direction characterized by (left). Traditionally, the argument of periastron is measured from East to the direction of the observer, however, the radial velocity (and the TTV) equations are invariant if we measure it from West. This causes no difficulties in the measurements of , and the East–West discrepancy can be resolved via time-series done via interferometry, astrometry or direct imaging.

In the text

	Fig. 2 Evolution of the maximum deviation of the observable and the predicted RV-measurements. This maximum value was selected as the maximum deviation over an orbital cycle. The observable RV curve was calculated with the parameters given in Table 1. The big change is due to the fact that after a while we assume completely opposite phase for the RV-fit if we neglect the apsidal motion.
In the text

	Fig. 3 Zoom to the first five years of Fig. 2. The typical averaged error bar of the available observations (see Sect. 5) is noted by the thick vertical line. The curve shows the maximum deviation during a revolution from an RV-curve without apsidal motion over time.
In the text

Fig. 4 Correlation plots and probability histograms for WASP-18AB RV+TTV fit. Red and orange areas denote the 1σ region, different gray areas denote the 2σ and 3σ regions. Numbers over the panels are the values of the Pearson-correlation coefficients between the parameters. The solid and dashed lines show the most probable solutions and the 1σ limits in thehistograms. Only a subset of the parameters are shown.

In the text

Fig. 5 RV-residuals of the model fit. x-axis: the index of the RV points used for the fit, in order of increasing time. y-axis: the residual of the fit in m s−1. The dashed black curve represents the residuals of the fit without any apsidal motion because we forced (including all relativistic and classical tidal terms), as well as Vtide = 0.0 m s−1. The solid red curve represents the residuals of the model with apsidal motion (see Eq. (19)). A significant improvement in the residuals can be seen. Table A.1 helps to identify the observational point via its index.

In the text

Fig. 6 Residuals of the model fit on the transit observations. We note that the transit observations were taken from the compilation of Maxted et al. (2013). x-axis: the index of the transit observations used for the fit, in order of increasing time (not the cycle number!). y-axis: the residual of the fit in minutes. The open black points represent the residuals of the fit without any apsidal motion because we forced (including all relativistic and classical tidal terms), as well as Vtide = 0.0 m s−1. The red filled circles represent the residuals of the model with apsidal motion (see Eq. (19)). For sake of clarity, the red symbols were shifted horizontally by 0.15 units. Table A.1 helps to identify the observational point via its index. No improvement can be seen in the transit data with the more complicated apsidal motion model which is contrary to the case of RVs (see Fig. 5.) where we could see significant improvement. This is in accordance with the expectation of Sect. 3.

In the text

Fig. 7 Residuals of the model fit on the occultation observations. Notice that the occultation observations were taken from the compilation of Maxted et al. (2013), too. x-axis: the index of the occultation observations used for the fit, in order of increasing time (not the cycle number!). y-axis: the residual of the fit in minutes. The meaning of symbols is the same as in Fig. 6. Table A.1 helps to identify the observational point via its index. A small improvement can be seen in the last occultation measurement with the complicated apsidal motion model relative to the simple model.

In the text