Issue 
A&A
Volume 623, March 2019



Article Number  A45  
Number of page(s)  12  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201834376  
Published online  06 March 2019 
An estimate of the k_{2} Love number of WASP18Ab from its radial velocity measurements
Deutsches Zentrum für Luft und Raumfahrt, Institut für Planetenforschung,
Berlin,
Rutherfordstr 2,
12489, Germany
email: szilard.csizmadia@dlr.de
Received:
4
October
2018
Accepted:
29
November
2018
Context. Increasing our knowledge of the interior structure, composition, and density distribution of exoplanets is crucial to make progress in the understanding of exoplanetary formation, migration and habitability. However, the directly measurable mass and radius values offer little constraint on interior structure, because the inverse problem is highly degenerate. Therefore, there is a clear need for a third observable of exoplanet interiors. This third observable can be the k_{2} fluid Love number which measures the central mass concentration of an exoplanet.
Aims. The aims of this paper are (i) to develop a basic model to fit the longterm radial velocity and TTV variations caused by tidal interactions, (ii) to apply the model to the WASP18Ab system, and (iii) to estimate the Love number of the planet.
Methods. Archival radial velocity, transit and occultation timing data were collected and fitted using the model introduced here.
Results. The best model fit to the archival radial velocity and timing data of WASP18Ab was obtained with a Love number of the massive (~10 M_{Jup}) hot Jupiter WASP18Ab: k_{2,Love} = 0.62_{−0.19}^{+0.55}. This causes apsidal motion in the system, at a rate of ~0.0087 ± 0.0033°∕days ≊ 31.3 ± 11.8 arcsec day^{−1}. When checking possible causes of periastron precession, other than the relativistic term or the nonspherical shape of the components, we found a companion star to the WASP18 system, named WASP18B which is a probable M6.5V dwarf with ~0.1 M_{⊙} at 3519 AU distance from the transit host star. We also find that small orbital eccentricities may be real, rather than an apparent effect caused by the nonspherical stellar shape.
Key words: techniques: radial velocities / planets and satellites: interiors / planets and satellites: individual: WASP18b / methods: data analysis
© ESO 2019
1. Introduction
In a binary system consisting either of two stars, or a star and a planet, tidal forces act on the nonmass point component(s). These tidal forces lead to several phenomena, including a decrease of the semimajor 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 semimajor 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 timevariable, 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 twolayer, coremantle 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 densityprofile 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 WASP18 data, to estimate the planetary fluid Love number. The presence of apsidal motion requires an eccentric orbit, and the noncircular nature of the WASP18Ab 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 nonsphericity 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_{2} is the 2nd order fluid Love number, a is the semimajor 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 – RossiterMcLaughlin 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 lefthand 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_{e}sin i value of the star, yields the stellar rotational period. Outoftransit 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_{2} 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 Lovenumber is defined by for example, Becker & Batygin (2013) as (3)
while Sterne (1939) used the socalled apsidal motion constant: (4)
where η_{2} 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 Lovenumber and we use simple k_{2} for the apsidal motion constant.
Fig. 1 Illustration of the apsidal motion. At time t_{0} the apsis (semimajor 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 timeseries 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 Zcoordinate of the star can be calculated stepbystep 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 semimajor 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 midtransit time T_{tr,N} at cycle N and that is connected to the anomalistic period as (up to first order in eccentricity, cf. Gimenez & GarciaPelayo 1983, Eq. (19) and Csizmadia 2009): (7)
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 firstorder approximation.
The mean anomaly M of the star is (9)
where t_{0} is fixed and t_{p} is the periastron passage time at the epoch.
Let t_{0} be the epoch of a well observed transit (in practice, we use a value determined from several transits observed in a short timewindow which can be more precise than just a single observation). At that time, the transit occurs when the true anomaly is (10)
ω_{0} is the argument of the periastron at the transit epoch. ϑ_{0} is a small correction due to the fact that the midtransit time occurs not at conjuction but at the smallest skyprojected distance of the star and the planet (see Csizmadia 2018; Gimenez & GarciaPelayo 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 midtransit occurred are (12) (13)
M and M_{0} 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_{1} is the semimajor axis of the star around the CMC and r_{1} its actual distance from CMC, then (17)
Concerning the Zcoordinate: (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 timederivative of the Zcoordinate plus an apparent tidal component, V_{tide}, which is nonnegligible in WASP18 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 halfamplitude (20)
while the second term in Eq. (19) comes from the timederivative 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 spacebased measurements (Csizmadia 2010). This precision is at least a factor of three worse for groundbased observations of shallow transits. Substitutingthe known system values of WASP18Ab (see Table 1) and assuming a Love number of WASP18Ab close to Jupiter’s k_{2,Love} = 0.34 then we get that the periastron precession rate should be ° per day^{1}. By virtue of Eq. (7) we determine that the corresponding TTVamplitude is about eP∕π ~ 3.7 min and the corresponding sinusoidal TTVcurve has a period of ~266 yr. Hence, the annual rate of change in the first ten years is about 0.01 s yr^{−1}. 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 WASP18Ab is stable to less than 6.6 × 10^{−7} 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 nonzero values for different timeperiods. Our estimate shows that for WASP18Ab this should be on the order of 100 m s^{−1} difference at its maximum after 2000 days of observations which is about eight times bigger than the typical observational error bars of one RVpoint (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 RVmeasurements. 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 RVfit 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 RVcurve without apsidal motion over time. 
Results of the M1–4 fits.
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 ω_{0} (Albrecht et al. 2011) –, , P and four RVoffset values between the different instruments. t_{0} 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 Algorithmbased Harmony Search optimizer (HS, Csizmadia 2018) with 1000 individuals in the population. We tried RVjitter values of 0 and 10 m s^{−1} and by a bisection method we adjusted the RVjitter – fixed during the HSrun – until we reached , N_{obs} being the number of RVobservations.
The results of the HS were refined by running ten chains of MCMC, with each chain consisting of one hundred million steps. A thinningfactor 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 Lovenumber, 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^{−1} (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 skyprojected 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 Lovenumber from Eq. (2) for different rotational statuses of the planet, listed in Table 2. That is why the uncertainty range of the Lovenumber is so high.
5.2. Fit of RV+TTV data with (M2)
In Model 2 we combined the RVdata 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 midoccultation 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^{2}. 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 midtransit 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 RVmodel 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 ω_{0}. We minimized the χ^{2} given by (21)
where the last eight terms are priors on the semimajor axis a, the ratio of the orbital period and the stellar rotation period P_{orb}∕P_{rot}, the scaled semimajor 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_{2} and anomalistic period.
One can notice that most of the correlations (Fig. 4) between the parameters are negligible except three pairs: k_{2} –period, k_{2} –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 MCMCchains. The correlation does not mean that we can measure only the ratio of k_{2} ∕P because the anomalistic period P appears in other terms without a relation to k_{2}. As it is clear from Fig. 4, the best solutions are concentrated onto a small part of the periodk_{2} 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 phaseshifts 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 RVjitter 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 ultrashortperiod planets.
We also carried out a fit without any apsidal motion and without apparent tidal RVterm (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 WASP18AB 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 Pearsoncorrelation 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 Lovenumber 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 WASP18Ab
Assuming a synchronous orbit for the planet P_{rot,planet} = P_{orb}, we derive an apsidal motion constant and a Love number for WASP18Ab: 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 semimajor 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 rotationalorbital 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 WASP18Ab 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_{2} 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_{2} are suggestive of a more homogeneous density distribution of a massive hot Jupiter, and the uncertainty range allows us to say that WASP18Ab 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 RVresiduals of the model fit. xaxis: the index of the RV points used for the fit, in order of increasing time. yaxis: 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 V_{tide} = 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). xaxis: the index of the transit observations used for the fit, in order of increasing time (not the cycle number!). yaxis: 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 V_{tide} = 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. xaxis: the index of the occultation observations used for the fit, in order of increasing time (not the cycle number!). yaxis: 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 WASP18Ab
WASP18Ab 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, Kamplitude 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^{−1} for WASP18Ab. Such spurious, mimicked eccentricities were also mentioned by Kopal (1959, Sect. V.4). In contrast to Arras et al. (2012) Eq. (V.122) of Kopal (1959) yields only ~2 m s^{−1} 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^{−1} of tidally induced RV variations, close to the 30 m s^{−1} predicted by Arras et al. (2012). But synchronous rotation is not the case for the host star WASP18A. 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_{1}, which is a good enough approximation here) which should be added to the radial velocity model (cf. Eq. (19)): (22)
μ and K stand for micrometer and Kelvin, respectively, and lastly (28)
q is the planettostar mass ratio, λ is the effective wavelength of the observations (we set 0.6 μm), τ_{0} 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 timevariable ω. 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^{3}. We consider these differences negligible.
Our ω_{0} = 257.27° ± 2.13° rules out the exact 270° argument of the periastron at the epoch which was fixed at T_{0} = 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 ω_{0}. We suggest obtaining new, very precise, high quality photometric and radial velocity measurements of WASP18Ab’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 WASP18A
We searched the second data release (DR2; Gaia Collaboration 2018) from the Gaia mission (Gaia Collaboration 2016) for close companions to WASP18A. This revealed an object 11.75 magnitudes fainter than WASP18A 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)^{4}, we suggest that the companion is a lateM type star (approximately M6.5V). Combining the visual separation with the distance to WASP18A (123.92 ± 0.37 pc), we find that the companion is physically separated from WASP18A 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 WASP18A. When we assume a circular, coplanar orbit with WASP18Ab 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 RVobservations. A substitution yields V_{max} ~ 0.1 m s^{−1} 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 earlierused WASP18b to WASP18Ab because it orbits around the brighter primary star, and the faint, newlydiscovered companion star is named here WASP18B.
7. Conclusions
In this study we show that there is tentative evidence for the presence of apsidal motion in the WASP18Ab 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 WASP18Ab 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,LoveJupiter} ~ 0.34) because of the large error bars, and compatible with k_{2,LoveJupiter} = 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 fiveyear long RV curve of WASP18Ab – 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.
WASP18A and its stellar companion.
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/141 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 NKFIHOTKA K113117 and NKFIKH130372 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
Radial velocity, transit and occultation timing fit residuals of WASP18Ab system.
References
 Albrecht, S., Winn, J. N., Johnson, J. A., et al. 2011, ApJ, 738, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Albrecht, S., Winn, J. N., Johnson, J. A., et al. 2012, ApJ, 757, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Arras, Ph., Burkart, J., Quataert, E., & Weinberg, N. N. 2012, MNRAS, 422, 1761 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, J. C., & Batygin, K. 2013, ApJ, 778, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Borkovits, T., Csizmadia, Sz., ForgácsDajka, E., & Hegedüs, T. 2011 A&A, 528, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Claret, A. 2004, A&A, 424, 919 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Csizmadia, Sz. 2018, MNRAS, submitted [Google Scholar]
 Csizmadia, Sz., IllésAlmár, E., & Borkovits, T. 2009, New Astron., 14, 413 [NASA ADS] [CrossRef] [Google Scholar]
 Csizmadia, Sz., Renner, S., Barge, P., et al. 2010, A&A, 510, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Einstein, A. 1915, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin: PAW), 831 [Google Scholar]
 Ferrero, G., Gamen, R., Benvenuto, O., & FernándezLajús, E. 2013, MNRAS, 433, 1300 [NASA ADS] [CrossRef] [Google Scholar]
 Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gimenez, A., & GarciaPelayo, J. M. 1983, Ap&SS, 92, 203 [NASA ADS] [CrossRef] [Google Scholar]
 Hellier, C., Anderson, D. R., Collier Cameron, A., et al. 2009, Nature, 460, 1098 [NASA ADS] [CrossRef] [Google Scholar]
 Iorio, L. 2011, MNRAS, 411, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Kellermann, C., Becker, A., & Redmer, R. 2018, A&A, 615, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kopal, Z. 1959, Close Binary Systems, The International Astrophysics Series (London: Chapman & Hall), 1959 [Google Scholar]
 Knutson, H. A., Fulton, B. J., Montet, B. T., Kao, M., & Ngo, H. 2014, ApJ, 785, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Kramm, U., Nettelmann, N., Redmer, R., & Stevenson, D. J. 2011, A&A, 528, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Love, A. E. H. 1911, Some Problems of Geodynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Martynov, D. Y. 1973, in Eclipsing Variable Stars, IPST Astrophysics Library, ed. V. P. Tsessevich (New York: John Wiley & Sons) [Google Scholar]
 Maxted, P. F. L., Anderson, D. R., Doyle, A. P., et al. 2013, MNRAS, 428, 2645 [NASA ADS] [CrossRef] [Google Scholar]
 McDonald, I., & Kerins, E. 2018, MNRAS, 477, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Ni, D. 2018, A&A, 613, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nymeyer, S., Harrington, S., Hardy, R. A., Stevenson, K. B., & Campo, Ch. J. 2011, ApJ, 742, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Padovan, S., Spohn, T., Baumeister, P., et al. 2018, A&A, 620, A178 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Ragozzine, D., & Wolf, A. S. 2009, ApJ, 698, 1778 [NASA ADS] [CrossRef] [Google Scholar]
 Rauw, G., Rosu, S., Noels, A., et al. 2016, A&A, 594, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schmitt, J. H. M. M., Schröder, K.P., Rauw, G., et al. 2016, A&A, 586, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Southworth, J., Hinse, T. C., Dominik, M., et al. 2009, ApJ, 707, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Southworth, J., Hinse, T. C., Dominik, M., et al. 2010, ApJ, 723, 1829 [NASA ADS] [CrossRef] [Google Scholar]
 Sterne, T. E. 1939, MNRAS, 99, 451 [NASA ADS] [Google Scholar]
 Torres, G., Andersen, J., & Gimenez, A. 2010, A&ARA, 18, 67 [Google Scholar]
 Triaud, A. H. M. J., Collier Cameron, A., Queloz, D., et al. 2010, A&A, 524, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wilkins, A. N., Delrez, L., Barker, A. J., et al. 2017, ApJ, 836, L24 [NASA ADS] [CrossRef] [Google Scholar]
 Zakamska, N. L., Pan, M., & Ford, E. B. 2011, MNRAS, 410, 1895 [Google Scholar]
 Zhou, G., Bayliss, D. D. R., KedzioraChudczer, L., et al. 2015, MNRAS, 454, 3002 [NASA ADS] [CrossRef] [Google Scholar]
A new study by Ni (2018) predicts k_{2,Love} = 0.53 for Jupiter. This does not change the estimate significantly.
http://var2.astro.cz/ETD/. We remark that Zhou et al. (2015) carried out a groundbased K_{s}band occultation measurement on WASP18Ab which would extend our baseline, but, unfortunately, they did not publish the midtime of the occultation, only the transit depth. We also excluded the measurements of Wilkins et al. (2017) obtained by HST because there were long gaps in their transit light curves, decreasing the precision.
http://www.pas.rochester.edu/~emamajek/ EEM_dwarf_UBVIJHK_colors_Teff.txt
All Tables
Radial velocity, transit and occultation timing fit residuals of WASP18Ab system.
All Figures
Fig. 1 Illustration of the apsidal motion. At time t_{0} the apsis (semimajor 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 timeseries done via interferometry, astrometry or direct imaging. 

In the text 
Fig. 2 Evolution of the maximum deviation of the observable and the predicted RVmeasurements. 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 RVfit 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 RVcurve without apsidal motion over time. 

In the text 
Fig. 4 Correlation plots and probability histograms for WASP18AB 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 Pearsoncorrelation 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 RVresiduals of the model fit. xaxis: the index of the RV points used for the fit, in order of increasing time. yaxis: 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 V_{tide} = 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). xaxis: the index of the transit observations used for the fit, in order of increasing time (not the cycle number!). yaxis: 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 V_{tide} = 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. xaxis: the index of the occultation observations used for the fit, in order of increasing time (not the cycle number!). yaxis: 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.