© The Authors 2024

1 Introduction

The understanding of planetary interiors constitutes one of the four major challenges of contemporary exoplanetary science (Schneider 2018). It is crucial for assessing their potential habitability, formation, and evolution. Interior processes play an essential role in creating and maintaining the physical conditions that are required to support life (see, e.g., Van Hoolst et al. 2019). They are useful for determining which planetary formation process is most plausible; namely, whether it is core accretion (Pollack et al. 1996 and references therein) or gravitational instability (Boss 1997 and references therein). According to the former, the accretion of planetesimals would result in the formation of a core comprised of heavy elements. On the other hand, planets forming through disc instability can rapidly accrete gas. These two scenarios lead to different core sizes. The size of the core is related to the second order fluid k_2,p Love number (Love 1911), as explained in Becker & Batygin (2013). Knowledge of the interior of a giant planet can help us distinguish between these two scenarios, in addition to its Love number, as we describe in this work.

Love numbers can also be used to break the known mass-radius-composition degeneracy of exoplanets (Baumeister et al. 2020). When an exoplanet is characterizable both by transit and radial velocity techniques, we can derive its radius and mass – and, thus, the planetary mean density. However, mass, radius and mean density are degenerate when it comes to determining the interior structure of bodies, such as the radial density, temperature, pressure, and composition profiles. We can find multiple solutions for the planetary interior producing the same planetary total mass and radius (see, e.g., Valencia et al. 2007; Wagner et al. 2011; Damasso et al. 2018). This kind of degeneracy can be decreased based on the strong assumption that the planet and its host star share the same metallicity (Dorn et al. 2015).

For Solar System planets, the degeneracy between mean density and interior structure can be reduced by in situ measurements of the gravitational field, as well as by seismic measurements, such as InSight (Banerdt et al. 2020), and by observations of moon motions. However, for exoplanets, gravitational moments cannot be measured. Hence, Love numbers were proposed as further observables in exoplanetary interior studies (see, e.g., Batygin et al. 2009; Ragozzine & Wolf 2009 and references therein), as has already been applied to eclipsing binary stars systems (see Russell 1928 and references therein). These authors measured the deformations and mass re-distribution inside the planet due to the tidal interaction with the host star. As shown by Baumeister et al. (2020), thanks to the knowledge of the second-order fluid Love number of the planet, we can better infer the distribution of possible thickness of each interior layer.

In this work, we employ a technique used by Csizmadia et al. (2019), which we further refined to determine k_2,p from periastron precession. This precession is derived from radial velocity (RV) measurements and transit and occultation mid-times. We apply it to the system WASP-19 based on an approach that requires the calculation of the secular evolution of the orbital elements, which is related to the mass distribution of the star and the planet through their second-degree Love numbers.

This paper is organized as follows. Section 2 and Appendix A present the theory of apsidal motion. Section 3 describes the system WASP-19 and provides evidence of a possible companion star in the WASP-19 system. In Sect. 4, we describe archival transit and occultation timing observations as well as the archival radial velocity investigations of the system which are supplemented by our recent RV study of the system. Section 5 contains our data analysis along with the search for a second planetary companion. Our conclusions are given in Sect. 6.

2 Apsidal motion

Apsidal motion, namely, the secular perturbation in the argument of periastron, can be detected either in RV or in TTV datasets (or their combination) to constrain the internal structure of planets. This type of analysis on RVs was introduced by Kopal (1959, 1978). Here, we make use of such equations, however, we have rewritten them in a way that allows them to be directly applied to the observed data sets. We describe their derivations more in detail in Appendix A.

The total rate of apsidal advance is the result of the three contributions:

(1)

The general relativistic term dω_GR/dt (Einstein 1915) is:

(2)

where G is the gravitational constant, M_⋆ the mass of the star, a and e the semi-major axis and eccentricity of the orbit, c the speed of light in vacuum, and P_a the anomalistic period.

The secular perturbations of the tidal and rotational components dω_tidal/dt and dω_rot/dt are derived in Appendix A and can be expressed more explicitly if approximated to the second order, as in Csizmadia et al. (2019):

(3)

and

(4)

In these equations, n is the mean motion, m_p and R_p the mass and radius of the planet, R_⋆ the radius of the star, and P_rot,p and P_rot,⋆ are the rotational periods of the two components. Also, k_2,p and k_2,⋆ are their second-order fluid Love number corresponding to the double the apsidal motion constant (Csizmadia et al. 2019):

(5)

both for the star (i = ⋆) and the planet (i = p). These two equations were obtained by Sterne (1939). Just as in Eq. (A.8), it is assumed that the rotational axis of both components are perpendicular to the orbital plane.

Regardless of whether the stellar and planetary rotation is direct or retrograde, the corresponding terms in Eq. (4) are positive, therefore they cause the apsidal line to advance.

Assuming the values of the orbital parameters of WASP-19Ab from literature and the Love number of the star from Claret (2004), we calculated the amplitude of the general relativity term and the stellar tidal and rotational contributions, as reported in Table 1. Moreover, a first estimate for the planet can be calculated assuming synchronous rotation for the planet and k_2,p between the values 0.01–1.5 (corresponding to the extreme cases of a mass point and a homogeneous body, respectively). The total contribution is between about 14.71–1689 ″d⁻¹ ≃ 0.0041–0.479° d⁻¹, respectively when assuming k_2,p=0.01 and k_2,p = 1.5. The apsidal motion period is therefore between ≃ 2.10 and ≃41.25 yr.

The equations in Appendix A show that the angle of periastron ω oscillates over one orbital period but also shows a gradual long-term change that can be retrieved in the RV observations. Due to the finite length of the typical datasets of observations (12 yr in the case of WASP-19Ab), in RV datasets, the change in the angle of periastron ω can be retrieved in the radial velocity of the star (see, e.g., Jackson et al. 2008 for more details on the timescales). In long baselines of observations, all orbital elements must be considered along with periastron precession. As described in Kopal (1978), they include the change in semi-major axis (tidal decay), eccentricity (circularization) and time of periastron passage. The semi-major axis and its eccentricity would diminish in value and the time of periastron passage is shifted in time.

3 WASP-19Ab

3.1 The system

The planet WASP-19Ab, hosted by the bright, active G-dwarf star (V = 12.3, Zacharias et al. 2012 and T_eff,⋆ ≃ 5568 K, Torres et al. 2012), was discovered through the transit method by the Wide Angle Search for Planets (WASP)-South observatory (Hebb et al. 2010). It is a hot Jupiter with mass m_p = 1.139 M_2/, inflated radius r_p = 1.410 R_2/ (Mancini et al. 2013) and orbital period P_orb ≃ 0.79 d.

So far, the system is among the best studied from the point of view of orbital and planetary parameters (e.g., Lendl et al. 2013; Mancini et al. 2013; Sedaghati et al. 2015, 2017) and atmosphere characterization (very pronounced water absorption, Huitson et al. 2013; Sedaghati et al. 2017, and titanium oxide detected in the transmission spectrum, Sedaghati et al. 2017). The proximity of WASP-19Ab to its host star results in significant tidal interaction that may govern the evolution of the planetary orbit, which makes the system a good candidate for our study.

Based on transit observations, certain studies have suggested the possibility of variations in the mid-transit time. Mancini et al. (2013) suggested a non-linear ephemeris. However, Petrucci et al. (2020) discarded orbital decay with 74 complete transit light curves spanning over ten years of observations. Moreover, as described more in detail in Sect. 5.1, Cortés-Zuleta et al. (2020) used transit timing variations (TTVs) to infer an upper mass limit for an additional candidate planetary companion in the system.

3.2 Companion star: WASP-19B

We searched the catalogue of the third Gaia data release (DR3) (Gaia Collaboration 2016, 2023) for possible co-moving companions to WASP-19A. We searched all objects within 4′ of the target, which corresponds to a physical separation of around 0.3 pc at the distance of WASP-19A, for those with compatible parallax (θ) and proper motions (PM). Only one object was found to have parallax and proper motion values in both the RA (α) and Dec (δ) directions that are compatible with those of the target to within 5σ. This object (Gaia DR3 ID = 5411725145916602496), which is 8.8 magnitudes fainter in the G-band, has a parallax that differs from that of WASP-19A by 1.6σ and PM_α and PM_δ values differing by just 0.1 σ and 0.6 σ respectively. This star is found in the southernmost part of the south-east region of WASP-19A, at a separation of 68″. The magnitude and BP-RP colour of the companion suggests that this star is a mid M-dwarf, with a best matching spectral type for the magnitude of M6V and for the colour M2V, using the updated table 1 of Pecaut & Mamajek (2013).

To estimate the probability of this putative companion arising from a chance matching of WASP-19 in parallax and proper motions, we searched a much larger area of sky in Gaia DR3, still centered on WASP-19, but with a radius of 2°. This produced 77 objects (including the companion; hereafter: ‘WASP-19B’) that match WASP-19 within 5σ in all of θ, PM_α, and PM_δ. This implies that there is an 8.6% chance of finding a matching object within our initial search radius of 4′, and just a 0.7% chance of finding an object as close as WASP-19B. Relaxing the 5σ limit to 3σ changes these probabilities to 2 and 0.16%, respectively. Furthermore, none of these additional matches, all of which are more than ten times further from WASP-19A than WASP-19B, demonstrate as good a match as WASP-19B. Mathematically, this is expressed as WASP-19B having the lowest value of , where for each parameter, X (where X stands for θ, PM_α, PM_δ) , where and are the values and uncertainties of parameter X for WASP-19 and its companion, respectively. We recognise that the uncertainty on the parallax of WASP-19B is relatively large. Even if WASP-19B is likely to be a true, gravitationally bound, companion to WASP-19A, some doubt remains as to whether it truly is a bound companion to WASP-19A. We adopted the planet name WASP-19Ab because it orbits around the brighter primary star. Details of WASP-19B are given in Table 2. With Eq. (9) of Csizmadia et al. (2019), we also determined that the effect of the companion star on the radial velocity of star A is negligible. Star B may be the cause of the presence of an eccentricity different from zero in the orbit on the transiting planet. While tidal interaction leads to circularization, the excitation of the orbit due to a distant companion might prevent it, as discussed in detail in Appendix B.

4 Observations

4.1 Transits and occultations

Mid-transit and occultation times used in this paper are listed in Tables C.1–C.2. We collected publically available data, including some amateur observations listed in the TRansiting ExoplanetS and CAndidates (TRESCA)² project database, part of the Exo-planet Transit Database (Poddaný et al. 2010). In particular, we included amateur data that were selected by Mancini et al. (2013) and other high-quality data showing a clear light curve without gaps. As detailed in the tables, we also re-analyzed some transits with the Transit and Light Curve Modeller (TLCM, Csizmadia 2020).

The time unit of each mid-transit and mid-occultation time was carefully checked and eventually converted into BJD_TDB (barycentric Julian date, in the barycentric dynamical time system). The difference between BJD_UTC and BJD_TDB depends on the number of leap seconds (Eastman et al. 2010). For our dataset, it is around one minute. If the time unit used in the reference paper was not clear, the main author was contacted and the issue was clarified; otherwise, the data were discarded and not reported here. We used the estimated mid-times for each primary transit and occultation to evaluate whether they exhibit deviations from a Keplerian orbit over a long baseline of observations (see Sect. 5.2 for more details) as well as to constrain the sidereal period of the planet.

4.2 RVs

We collected all available radial velocity measurements from 2008 onward and we present new RV observations (see ID 8 below). Information on the RV observations are listed in Table 3 and the complete datasets are shown in Appendix D.

Here, we list the main details on the datasets, with the labels corresponding to the discussion in Appendix D:

• ID 1: Hebb et al. (2010; discovery paper) reported 34 RV measurements obtained with the CORALIE spectrograph on the 1.2 m Euler telescope between 2008 and 2009. To the best of our knowledge and after contacting the authors of the paper, the dataset reported in the paper is in BJD_UTC. We present the data transformed to BJD_TDB;

• ID 2: Hellier et al. (2011) reported 3 RV points from the CORALIE spectrograph on the Swiss Euler 1.2 m telescope in 2009. The data were reported in BJD_UTC in Hellier et al. (2011), so they are transformed here into BJD_TDB;

• ID 3: Hellier et al. (2011) reported 36 points during and around more than one transit in 2010 using the HARPS spectrograph on the ESO 3.6 m at La Silla, taken in order to measure the Rossiter-Mclaughlin effect. As well as for ID 2, the data were reported in BJD_UTC in the original paper, so they have been transformed here into BJD_TDB;

• ID 4: Albrecht et al. (2012) reported 20 data points taken in 2010 as measurements of the Rossiter-McLaughlin effect with the Magellan II (Clay) 6.5 m telescope and the Planet Finder Spectrograph (PFS);

• ID 5: Knutson et al. (2014) reported 3 RV measurements obtained with the High Resolution Echelle Spectrometer (HIRES) on the 10 m Keck I telescope between 2008 and 2013;

• IDs 6,7: Sedaghati et al. (2021) reported 88 RV points taken with VLT/ESPRESSO between 2019 and 2020, used to derive the Rossiter-McLaughlin effect and characterize the atmosphere of the planet. The two different IDs are given because of some setup changes in ESPRESSO, as described by Sedaghati et al. (2021), which lead to the introduction of a small offset. Moreover, as clarified with the first author of that paper, ID 6 was reported in BJD_UTC in the original paper, while ID 7 in BJD_TDB, therefore the first one is transformed;

• ID 8: reported 19 RV points taken by us using the High Accuracy Radial velocity Planet Searcher (HARPS, Mayor et al. 2003) at the ESO La Silla 3.6 m telescope under the programme 0104.C-0849 (PI: Sz. Csizmadia). The observations were performed between 6 and 21 February 2020. We tried to obtain 1–3 points per night and we planned the observations to have a good phase coverage. The exposure time was 1200 s while the spectral resolution as 115 000. The objAB observing mode of HARPS was used, covering the 378–691 nm wavelength range. The measured signal-to-noise ratio (S/N) varied between 10.2 and 21.6 depending on weather and seeing conditions. The data were reduced by the standard HARPS pipeline available at La Silla Observatory and the results with their corresponding uncertainties are reported in Tables D.1.

All RV observations are shown in Fig. 1. As described in Sect. 5.2, the datasets have different offsets {which are subtracted in the figure. They are fitted in our model by including the term V_i,instr in Eq. (7). As for the transit and occultation mid-times, the time unit of the observation was transformed and uniformized to BJD_TDB when necessary.

Fig. 1 Radial velocity observations of WASP-19Ab, time vs. radial velocities, with corresponding observational errorbars. The different data sets (see Sect. 4) are represented by different colours. Between the datasets, there are significant offsets between different instruments and studies. They have been removed accordingly to the results in Table 4.

5 Data analysis

5.1 Long-term effects in RVs

In this section, we study the long-term effects affecting the RV curve. If present, they must be identified and subtracted before modeling and estimating the rate of periastron precession in Sect. 5.2.

First, we applied the l₁ periodogram (Hara et al. 2017) to the archive and newly acquired RV data (ID 8, see Fig. 2) to search for further companions in the system, in addition to the confirmed planet in the system. This method was developed specifically to find companions in RV datasets and can be used in a similar way to a Lomb–Scargle periodogram. To verify the significance of the peaks, we followed the re-sampling approach suggested by Hara et al. (2017): we randomly removed 10–20% of the data and re-computed the periodogram. Some peaks did not appear after the re-sampling, therefore, we did not consider them real. After the re-sampling, we identified only one peak at 0.78 d, as shown in the top plot of Fig. 2. It was retrieved in correspondence with the orbital period of WASP-19Ab, with a false alarm probability (FAP) of log₁₀(FAP) = −30.9 (FAP ≃ 10⁻³¹). The analysis shows no evidence of any additional periodic signal in the RV data.

We also re-computed the periodogram on the RV residuals after the subtraction of the signal of the first planet (see bottom plot in Fig. 2). After a re-sampling of the data, only a few small peaks are left. However, their FAP is ⩾0.2, meaning that any peaks are due to artefacts.

We estimated the sensitivity of our search for perturbing bodies in the data by performing a sensitivity study based on the transit timing variations. We calculated the maximum mass of a perturber in the system as a function of the maximum amplitude of the TTVs (see Eq. (2) in Holman & Murray 2005). Our analysis shows that the maximum allowed mass of the perturber, M_c, as a function of its orbital period, P_c, assuming different values of the orbital eccentricity (e = 0.0, 0.2, 0.4,0.6, and 0.8). The analysis shows that the region above 4000 M_2/ and period below 20 days is excluded by the TTV observations.

We also performed a sensitivity study based on the RVs. The maximum amplitude of the RV residuals after the subtraction of the signal of the orbit of planet b (with Sedaghati et al. 2021’s parameters) is around 52 m s⁻¹. Through Kepler’s third law, we estimated the maximum mass of the putative perturber. With an orbital period equal to 12 yr (the total length of the RV dataset), we obtained M_c ≃ 3.9 M_2/ and an orbital distance of ≃ 5 AU. We fitted the data with an additional Keplerian signal (beyond the signal of planet b) mimicking a hypothetical candidate for planet c. We injected a candidate and repeated the experiment two billion times with different random eccentric orbital parameters and fitted the new HARPS RV data. We varied the candidate parameters as follows: and between −1.0 and 1.0, the RV semi-amplitude K_c between 0.1 and 100 m s⁻¹ (to widely investigate the amplitude of the RV residuals with errorbars) and the orbital period between 0.1 and the length of the dataset (≃19 d for the first attempt with the only HARPS dataset and ≃12 yr for the second attempt with the whole dataset). Moreover, from Kepler’s third law, we set the following condition to ensure that the orbit would not collide with the star:

(6)

where G is the gravitational constant, while R_⋆ and M_⋆ are the radius and mass of the host star, respectively. With these constraints and varying the orbital parameters in the aforementioned ranges, the candidate mass range is found to be between ≃1.4 and ≃400 M_⊕.

The results of our simulation (see Fig. 3) suggest scenarios that satisfy the condition χ² ⩽ 500.The top plot shows the χ² of the model when including the proposed perturber, corresponding to five more free parameters: transit epoch, period, RV semi-amplitude, and (plotted as a function of the orbital period). Many of these candidates are discarded because the corresponding χ² is too high or the amplitude of the perturbing body on the RV curve would be over a dozen m s⁻¹, so they would definitely be visible in the RV data and in the periodogram. The red lines underline the periods corresponding with the stable orbital regions present in the stability map by Cortés-Zuleta et al. (2020). Their analysis allowed for the presence of a second planet either in 1:2, 5:3, 2:1, 5:2, or 3:1 resonance with planet b, corresponding to orbital periods of 0.39 d, 1.31 d, 1.58 d, 1.97 d, or 2.37, respectively. Some of these periods correspond with regions in the figure that are densely populated by our best candidates. The bottom plot shows the mass of the candidate plotted against its orbital period. Again, we highlight the stable orbit by Cortés-Zuleta et al. (2020) with red dashed lines.

The expected χ² would be around the number of degrees of freedom (19 total datapoints for dataset 8, of which 15 were fitted as out-of-transit points and 11 free parameters), while the lowest χ² we obtained during the simulations is ≃ 134 (clearly visible in the top plot of Fig. 3 at ≃10 d orbital period). We fit the HARPS data with a one-planet eccentric orbit model to compare it to the scenarios with a perturbed. Without adding any RV jitter, the χ² of this model is ≃35. Since the two-planet model has four more free parameters, by applying the Bayesian information criterion (BIC), we found no evidence of the presence of a perturber in the system with a period within the length of the dataset with this approach. For weak (respectively, strong) evidence (see Kass & Raftery 1995 for more details) that the two-planet model is preferred, a χ² value below 115 (respectively, 98) is needed. We also injected a candidate in the whole RV dataset (IDs from 1 to 8), but we obtained an even higher χ², resulting in an even more unlikely presence of a perturbed. Our results exclude it up to an RV amplitude of ≃1 m s⁻¹.

This simulation confirms the conclusion of Knutson et al. (2014) when looking for the acceleration of the system due to possible long-period companions. The cited paper does not find an acceleration that differs from zero by more than 3σ. In this way, we extended the orbital period range of investigation in Knutson et al. (2014) from between 0 and 1702 d up to ≃4200 d.

Since we did not find any additional companion and there is no evidence for long-term trends in the data, we proceeded with our analysis by introducing periastron precession in the fit.

Fig. 2 l1 periodogram of all RV data outside the transit (top). There is only one significant peak with false alarm probability ≃ 10−31. The periodogram of the residuals after the subtraction of the contribution of planet b (bottom). All residual peaks have false alarm probabilities of >20%.

Fig. 3 Properties of third bodies excluded by RV data. The plot shows the results for the simulated scenarios including a putative companion planet c of different masses at different periods. Each black dot represents one of the around 20 000 selected scenarios (see the text for more details). Top plot: χ2 value vs the orbital period of the perturber. The red lines correspond to the period of the stable orbits predicted by Cortés-Zuleta et al. (2020). The two horizontal dashed lines represent the χ2 required to have weak (upper line) or strong (lower line) evidence of a second planet, as discussed in the text. Bottom plot: mass and orbital period of the selected candidates.

5.2 Rate of apsidal motion

In this section, we describe the code used for the data analysis (by allowing for non-zero periastron precession) with the fitted parameters and the results.

We used the approach and idl code developed by Csizmadia et al. (2019) to analyze the RV data in parallel with mid-transit and occultation times. Since then, the code has been modified and updated. Now the Genetic Algorithm minimization is followed by a differential-evolution Markov chain Monte Carlo (DE-MCMC, Ter Braak 2006) algorithm to explore the parameter space. The values of the fitted parameters (see Table 4 for the list of such parameters) obtained by the Genetic Algorithm are used as a starting point for the DE-MCMC analysis to derive the posterior probability distributions. The Genetic Algorithm is initialized with a population set made of 200 individuals. In the DE-MCMC we run 90 chains. The first 6 × 10⁴ steps of each chain are discarded as part of the burn-in phase. They are followed by 10⁴ more DE-MCMC steps, which are used for the posterior distributions. The convergence is checked by the Gelman - Rubin test applied on the steps that follow the burn-in phase. At the end of the run, all parameters resulted in having R < 1.2 (Gelman & Rubin 1992).

The priors that are used on the fitted parameters are set from the results in the literature (mainly from Hellier et al. 2011) and noted in Table 4 along with the results. The prior on the eccentricity is set as the weighted mean and standard deviation of the results in the literature (Hebb et al. 2010; Anderson et al. 2010, 2013; Hellier et al. 2011; Mandell et al. 2013; Lendl et al. 2013; Zhou et al. 2014; Knutson et al. 2014). From this, the priors on and are calculated as a combination of the two parameters. The contribution on the χ² of the priors was calculated the same way as in Csizmadia et al. (2019).

We modeled the radial velocity of the host star with a Keplerian motion but with a time-variable ω:

(7)

with a secular advance of the longitude of periastron, , caused by the tidal and rotational potential (as expressed in Eq. (A.8)):

(8)

where i is the index of the instrument used and j is the index of the time when the observation was taken; V_γ is the systematic velocity of the barycenter of WASP-19Ab system relative to the barycenter of the Solar System. Different datasets may have different instrumental zero-point offsets, V_i,instr. Finally, K is the radial velocity half-amplitude:

(9)

In this expression, a_⋆ is the semi-major axis of the star around the common center of mass and i is the inclination between the orbital plane and the plane of the sky (perpendicular to the line of sight).

Here, δV is the non-orbital apparent contribution that arises from the rotation of the star around its own axis. The stellar and planetary shapes are distorted ellipsoids with the longest axis in the direction of the radius vector connecting them, as shown in Fig. 4. Their apparent area varies continuously during an orbital revolution, therefore the emitted light changes (i.e., the ellipticity effect). Moreover, in close binary systems, part of the radiation of each of the two components falls on the surface of the companion. It is absorbed and re-emitted (or scattered) in all directions, varying during the orbital phase (i.e., the reflection effect).

These two photometric effects in the RV data due to distortion are described by Kopal (1959) as a non-orbital contribution to the RV curve:

(10)

where the Love number of the star is assumed from theoretical calculations (we take it from Claret 2004). Then, P_rot,⋆ is the rotational period of the star around its axis and f₂ takes into account the limb darkening. For more details on f₂ see Kopal (1959) and Csizmadia et al. (2019).

Arras et al. (2012) already includes such calculations for precise RV analysis of WASP-18A, but neglects the effect of stellar rotation. In other words, a synchronous rotation of the star is assumed. However, for WASP-19A this is not the case. It is known to be active, showing a rotational modulation with P_rot,⋆ = 10.5 d (Hebb et al. 2010). Therefore, the term P_orb/P_rot,⋆ is not expected to be 1, but ≃ 0.79 d/10.5 d ≃ 0.08, introducing a factor of about 12 difference in the amplitude of such effect calculated by Arras et al. (2012) and Kopal (1959). The amplitude of this contribution is around 0.19 m s⁻¹, so it is negligible for WASP-19A if compared to the typical observational error-bars. The period of this effect is linked to the sine term of the true anomaly in Eq. (10) and, therefore, it is half of the orbital period.

Moreover, we added jitter terms for each spectrograph to the RV model (see, e.g., Baluev 2009 for a discussion of instrumental jitters). The jitter is composed of an instrumental and a stellar contribution. It is not possible to separate them. As a starting point, if present, we used the values proposed by the corresponding literature. In particular, 13 m s⁻¹ for Hebb et al. (2010)’s data (ID 1, as in Anderson et al. 2010), 14.1 m s⁻¹ and 6.9 m s⁻¹ for Coralie (ID 2) and HARPS (ID 3) Hellier et al. (2011)’s data, respectively, 20 m s⁻¹ for Albrecht et al. (2012)’s (ID 4) and 17.8 m s⁻¹ for Knutson et al. (2014)’s (ID 5) RVs. Then, we refined all the jitter values by minimizing the χ² as:

(11)

where O_i and C_i are the observed points (see Table D.1) and the (calculated) expected values from the model, respectively; σ_i represents the observed errors corresponding to O_i and n_data is the number of datapoints in the subset. This procedure was iteratively repeated until the condition of Eq. (11) was fulfilled. Due to the highly inhomogeneous set of data, this procedure allows for every RV dataset to have a weight in the fit proportional to the number of its data points.

We also added a jitter value for transits and one for occultations to account for inhomogeneities in the dataset and for a perturber with an orbital period longer than 12 yr, as described in Sect. 5.1. The fitted values for the jitter terms can be found in Table 4, along with all the fitted and derived parameters in our model.

The orbital eccentricity is a key issue for Love number studies: in circular orbits, apsidal motion is not observable. In addition, the argument of periastron is hard to measure if eccentricity is not well known. When the eccentricity is small, then its uncertainty range is often large and the investigators find that the orbit is compatible with a circular one. However, such results in the literature do not necessarily imply that the orbit is circular, despite many investigators setting a zero eccentricity (and not propagating the error bars of eccentricity to the corresponding result). It only means that an upper limit exists for eccentricity.

Moreover, a certain dataset may indicate the presence of an eccentric orbit with an accurate result. However, the sample size is not big enough or not aptly distributed to get the uncertainty ranges accurate enough. The results are imprecise and the real value of eccentricity is vague.

Although we collected more RV and TTV data than any previous study of the WASP-19Ab system, the datasets come from various instruments with different precisions and, thus, are non-homogeneous. The value of eccentricity may be affected by the different data accuracy, as some of the data sets are diluted by less precise measurement points. Because of these reasons, any previous estimate of eccentricity can be accurate but imprecise.

Applying a prior on the eccentricity in the fitting procedure can be a valid approach. The average and standard deviation of the previously reported eccentricities can better represent reality than any single determination from a limited dataset. Therefore, we calculated the weighted average of the previous eccentricity measurement and their weighted standard error using their uncertainty range as weights. Therefore, we performed a fit using this mean eccentricity as a prior with their weighted standard error as a Gaussian width.

In addition, in the calculation of the χ² of the fit, priors also on the stellar temperature, T_eff,⋆, and Love number, k_2,⋆, the mass ratio, q, the orbital inclination, i, the orbital over rotational period ratio, F_⋆, and the ratio, r_p/R_⋆, were added as described in Csizmadia et al. (2019).

Together with the RVs, we modeled mid-transit and mid-occultation times listed in Tables C.1–C.2 mainly to refine the determination of the first transit epoch and the period. To include periastron precession in the fit of predicted mid-times, we followed the approach of Gimenez & Garcia-Pelayo (1983) and Csizmadia et al. (2019).

Fig. 4 Star-planet system geometry when tidal interactions occur. The spherical shapes that star and the planet would have in the absence of tidal interactions are drawn in red and dark brown. The tidal bulges are drawn in orange and light brown, respectively. They arise on the stellar and planetary surfaces (respectively in orange and light brown) and rotate with the star itself and are elongated in the direction of the other body. The tidal lags are also drawn. The image is not to scale.

Fig. 5 χ2 values of the fit when varying in the model. The horizontal dashed line separates the statistically worse models (over the line) from the best models (below the dashed lines). The models below the dashed lines are statistically equivalent within 1σ (68% confidence level for the red line) and 2σ (95% confidence level for the green line). Note: the significantly larger χ2 at . The high amount of local minima is an indication of the complexity of the problem.

5.3 Presence of apsidal motion

Because of the large parameter correlations between the mean motion and periastron precession rate (see Csizmadia et al. 2019), as a first analysis, we calculated 80 different models with fixed values of , the so-called semi-grid analysis. This procedure also helps to get an initial range for the fitted value of . The limiting values were and + 0.45°d⁻¹. We also fixed the epoch to 2 454 775.3372 BJD_TDB, while all other parameters were free in the DE-MCMC analysis. The result is presented in Fig. 5, showing the assumed value of the periastron precession rate and the corresponding χ² of the fit. The figure clearly shows that the joint fit of the RV and TTV data is significantly worse if , namely, no periastron precession is present. A small amount of periastron precession improves the quality of the fit significantly. The figure shows that small negative periastron precession rates also improve the fit but still result in a statistically worse fit than the small positive precession values. From that figure and the statistical analysis, we conclude that small positive periastron precession rates are needed to have a good fit of the joint RV and TTV data of WASP-19Ab, and we are able to give a lower and an upper limit of . The result is refined in the next section.

6 Discussion of the results

6.1 Apsidal motion rate

In this section, we show and discuss the results of the analysis performed with the model described in Sect. 5. Contrarily to the semi-grid method (Sect. 5.3), here we fit along with the other parameters.

During the analysis, we noticed that different DE-MCMC chains converged to different values of . In particular, as it is shown in Fig. 6 and in Table 5, during the run multiple chains explored various positive and negative values of . Figure 6 clearly shows also the correlation between the anomalistic period and the periastron precession rate, as discussed in Csizmadia et al. (2019). It is due to the formula that links the anomalistic period P_a and the sidereal period P_s:

(12)

This formula that links the two quantities gives a between period and periastron precession rate, but no unique solution. The added value of including transit and occultation mid-times lies in the fact that they constrain the sidereal period, and therefore a range of values for P_a and (as shown in Fig. 6), while avoiding some other ranges.

A different chain convergence shows why it was important to run various chains. Only one chain would have converged to one of the possible values, excluding the other results or could have been trapped in a local minimum (the same local minima found by the semi-grid method approach; see Fig. 5).

The highest peak of (solution B) in the posterior distribution results around approximately +0.065°d⁻¹ (summarized in Table 5). The positive solutions C, D, and E can be discarded because of their high χ² values. Moreover, solution E would lead to a planetary Love number higher than 1.5, which is unphysical.

The negative periastron precession rate (solution A) would correspond either to a retrograde orbit (of the planet with respect to the rotation of the star) or to an unseen perturber. The retrograde orbit can be excluded from Rossiter-McLaughlin effect studies and TESS light curves (Wong et al. 2020). They confirm that the spin axis of the star is almost perpendicular to the orbital plane.

We investigated in detail the presence of a perturber in Sect. 5.1 and no evidence is found. Assuming a negative periastron precession amplitude as observed of approximately −0.064ºd⁻¹ = −229.68″d⁻¹, through Eqs. (12) and (16) in Borkovits et al. (2011), we can estimate the period of the body producing such a precession rate. We subtracted the general relativity contribution (calculated in Sect. 5.2 as 2.85″ d⁻¹) from these solutions and derived the orbital period of the perturber producing such an . For a circular orbit, a body with a mass in the range 0.1 M_⊕ and 20 M_2/ would have an orbital period below 10 d, as shown in Fig. 7. Such a planet would be retrievable in the RV data, while no long-term (over a 10-day orbital period) companion in this mass range can produce a negative periastron precession. Therefore the solutions yielding negative values can be discarded.

The results for the fitted parameters of the -positive solution (solution B) with lowest χ² are shown in Table 4, along with their errorbars. With 99 RV points (outside the transit, 203 in total), 162 mid-transit and 11 mid-occultation points, and 19 free parameters, the reduced χ² is 1.02. The observed value of the total contribution from the DE-MCMC analysis of the data is (see Table 4), corresponding to an apsidal motion period of .

Figure 8 shows the residuals of the transit and occultation mid-times with such a fit. Figure 9 shows the phase-folded RV curve and the residuals of our periastron precession fit.

We subtract the General Relativity term and the tidal and rotational components relative to the star from the observed . What is left in Eq. (1) are the tidal and rotational components due to the planet. They are related to the interior of the planet and show clearly the dependence on Love number (see Eqs. (4) and (4)). Then, k_2,p can then be estimated by inverting these equations.

The only unknown parameter is the rotational period P_rot,p of the planet itself. Since there is no measurement of the rotational period of the planet, we can only assume a value for P_rot,p. The two extreme cases are presented when:

• tidal forces lead to the disruption of the planet (P_rot,p = 0.14 d), then we obtain k_2,p = 0.07;

• the planet is non-rotating (P_orb/P_rot,p → 0), then we obtain k_2,p = 0.22.

The measured Love number of the planet must lie in between these two values. We can suppose that tidal forces synchronize the orbital and rotational periods quite fast. This means that P_orb = P_rot,p. We can then derive k_2,p from Eqs. (A.2) and (A.7): , with a precision of 15%.

For comparison, the most up-to-date values of the tidal Love number for planets and satellites in the Solar System and exo-planets are listed in Table 6. The value obtained for the fluid Love number of WASP-19Ab has the same order of magnitude as the other values derived for exoplanets.

The lower the value, the more matter is condensed toward the center of the body. The low Love number obtained, k_2,p = 0.20, is an indication of a high-mass core. This has to be tested and verified with simulations of planetary interior structure through the equations of state, however, this is beyond the scope of this paper.

Moreover, if hydrostatic equilibrium is assumed (whereby the surface of the body can be described as an equipotential surface of the tidal and rotational potentials), the second-order fluid Love number h_2p can also be calculated. This assumption is close to reality if the body is a giant gaseous planet. In this case, h_2,p = k_2,p+1 (Munk & MacDonald 1960), therefore h_2,p = 1.20 for synchronous rotation. This causes a small-shaped deformation that can be hard to observe photometrically (Hellard et al. 2019).

During the analysis, we noticed that the main contribution to the uncertainty of the fit is related to the transit fit. It is due to underestimated errorbars on the mid-transit times, which compromise the reliability of the results, increasing the χ² of the fit. This was already suggested by Mancini et al. (2013) regarding the 54 transits they analysed. Therefore, we increased the errorbars to 0.0005 d (approximately 43 s) if the reported uncertainties were lower than this threshold.

We also underline the importance of acquiring new light curves. Precise transit and, mostly, occultation mid-times are needed to break the degeneracy between periastron precession and the change of orbital period due to for instance companions or orbital decay.

Fig. 6 Posterior distributions and correlation between the anomalistic period and periastron precession rate where the chains converged. The middle plot shows the strong correlation between these two free parameters, as in Eq. (12). On the sides, are the corresponding histograms of the posterior distributions. As can be seen from the histograms on the sides, many chains converged to various positive and negative peaks, as predicted by the semi-grid method and reported in Table 5.

Fig. 7 Mass and orbital period of a second body that would cause a negative periastron precession rate of approximately −0.0638°d−1 (corresponding to the negative solution of the DE-MCMC analysis) on planet b, as deduced from Borkovits et al. (2011). See text for more details.

Fig. 8 Residuals of the periastron precession fit for transits and occultation mid-times, calculated with the sidereal period. Transit mid-times are plotted as blue circles, occultation mid-times as red circles. The errorbars include the jitters, as in Table 4.

Fig. 9 Data with errorbars, phase-folded with our fit including periastron precession, shown in the upper plot. The colours represent different datasets, as in Fig. 1 (see also Sect. 4.2 for more details). Residuals of the periastron precession rate fit are shown in the lower plot. The error-bars include the jitters, as in Table 4. In both plots, the in-transit points are not shown as they were not used in the fit.

6.2 Alternative explanations

We investigated alternative scenarios to the presence of periastron precession. In particular, we searched for distant massive companions with the approach of Knutson et al. (2014). The long-term acceleration of the system was fitted by adding a γ acceleration to the RV Eq. (7) as:

(13)

We ran the DE-MCMC algorithm in the case of a circular and an eccentric orbit. The posterior distribution of has median of (7.37 ± 5.02) 10⁻⁷m s⁻² for the circular case and (7.43 ± 4.89) 10⁻⁷ m s⁻² for the eccentric one. Both results are compatible with zero and reproduce the results of Knutson et al. (2014) of (7.52 ± 3.94) 10⁻⁷ m s⁻². Our results increase the confidence level of a non-detection of long-term acceleration from Knutson et al. (2014)’s 1.9–1.5σ.

We also investigated the tidal decay scenario, as it has been analysed and confirmed in the system WASP-12b (see respectively Patra et al. 2017; Yee et al. 2020; Turner et al. 2021; Wong et al. 2022). The observed TTVs can be due to tidal decay or periastron precession. We fitted the RV semi-amplitude K in Eq. (9) as a time-varying parameter:

(14)

Deriving Kepler’s laws with respect to time, we can see that a change in the RV semi-amplitude results in an opposite change in the semi-major axis.

For the circular orbit case, we obtained while for the eccentric case (−2.60 ± 2.73) × 10⁻⁸ m s⁻². As Rosário et al. (2022), we did not find significant evidence of orbital decay. Both in the case of a long-term velocity and of tidal decay, the results are compatible with zero and therefore considered as constant.

7 Conclusions

In this study, we review the theory on tidal interaction, as described by Kopal (1978) and reported with modern notation. Using assumptions on the Love number of the star, based on well-established theoretical models of its evolution and interior structure, the apsidal motion rate allows us to constrain the Love number k_2,p of the planet. We applied the results to a case study, namely, the system WASP-19Ab, by using archival mid-transit times, occultations, and radial velocity data, including previously unpublished RV observations made by HARPS. This dataset covers twelve years in total. First, we ruled out previous claims of additional planetary companions in the system as possible perturbers in the radial velocity curve that could hide the small contribution of the periastron precession motion. Moreover, we investigated the presence of a long-term acceleration in the RV dataset and of tidal decay. Both these scenarios have been rejected with respect to being compatible with a zero trend. Then, we used a model that accounts for an eccentric orbit and apsidal motion in order to fit the data. We detected the presence of apsidal motion in the planetary orbit and linked it to the Love number of the planet. By assuming synchronous rotation for the planet, we derived

The planetary fluid Love number is constrained in this system for the first time. There are few systems for which it has been measured. We take the opportunity to compare WASP-19Ab to them. Table 6 shows the fluid and tidal Love numbers of Solar system planets and bodies and exoplanets which are reported in the literature. The derived value of the fluid Love number k_2,p of WASP-19Ab is on the same order of magnitude as the one of other Jupiter-like planets (WASP-4b, WASP-18Ab, WASP-103b and WASP-121b). A low value suggests the presence of a dense core in the planet, similar to the case of Saturn, Neptune, and WASP-121b. Note that the Love numbers of WASP-12b and WASP-4b (derived by Campo et al. 2011 and Bouma et al. 2019, respectively) are intentionally left out of the table. In the first case, the TTVs might be explained by tidal decay (Yee et al. 2020 and reference therein). In the second case, the presence of a second, long-period planet in the system (Turner et al. 2022) further complicates the TTV analysis. The results of Bouma et al. (2020) were discussed further by Harre & Smith (2023) in the light of this second planet, making the determination of k_2,p uncertain.

No general conclusion on planetary physics can be drawn from this set of measured Love numbers because the sample is too small to be representative of the existing exoplanet population. We need observations of the Love number of more exoplanets.

We also characterize, for the first time, the stellar multiplicity of the system. We find a stellar companion with an apparent separation of 68″, corresponding to a physical distance of almost 20 000 AU, from the system under analysis. However, its expected effect on the RV of WASP-19A is negligible.

Acknowledgments

We acknowledge the support of DFG Research Unit 2440: “Matter Under Planetary Interior Conditions: High Pressure, Planetary, and Plasma Physics” and of DFG grants RA 714/14-1 within the DFG Schwerpunkt SPP 1992, Exploring the Diversity of Extrasolar Planets. L.M.B. would also like to thank Carl Knight for providing us with two transits of WASP-19Ab observed by the Ngileah Observatory, New Zealand, in May 2022. L.M.B. would also like to acknowledge H. Hußmann and F. Sohl for the useful discussion on the different kinds of Love numbers. The work is based on observations collected at the European Southern Observatory under ESO programme 0104.C-0849. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This work presents results from the European Space Agency (ESA) space mission Gaia (https://www.cosmos.esa.int/gaia). Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC) (https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA).

Appendix A Equations of apsidal motion

Fig. A.1 Geometry of the orbital system in the simple case of orbital inclination of i =90°. The star symbol and the open dot denote the positions of the star and the planet, respectively; ω is the angle of periastron, θ is the phase from superior conjunction, and v is the true anomaly; l.o.s. denotes the line of sight. In this configuration, θ equals 0 for transits and π for occultations.

In this Appendix, we rewrite the equations derived by Kopal (1959, 1978) with modern notation.

First, we define the stellar and planetary disturbing accelerations ℝ_i and 𝕊_i, where the subscript i represents either the star (i = ★) or the planet (i = p). ℝ_i acts in the direction of the radius vector and 𝕊_i in the plane of the orbit in the direction perpendicular to the radius vector and positive in the direction of motion.

The radial perturbation accelerations can be expressed as:

(A.1)

and

(A.2)

where G is the Gravitational constant, M_★ and m_p are the masses of the star and the planet, respectively, R_★ and r_p their radii and k_2,★ and k_2,p their second-order fluid Love number. Here, r is the radius vector:

(A.3)

where a is the semi-major axis of the orbit, e is the orbital eccentricity and υ is the true anomaly of the planetary orbit, as shown in Fig. A.1.

A tidal lag in latitude arises if the spin axis of the two bodies is not perpendicular to the orbital plane and in general is small or even zero. On the other hand, a lag in longitude arises from asyn-chronism between spin and orbital rotation. The angles ϵ_★ and ϵ_p represent the tidal lag of the star or the planet in astrocentric longitude (Efroimsky & Makarov 2013) as:

(A.4)

where P_orb is the orbital period and P_rot,i is the rotational period of the star or the planet. The function sign takes the value of +1, 0 or -1, depending on the sign of the expression in the square brackets. Q_i is the tidal quality factor. A high Q factor corresponds to a weak dissipation and vice-versa (see Kaula (1964) and reference therein).

The radius vector, r_ϵ, is defined as:

(A.5)

The difference between r and is due to a finite speed of progression of the radial tidal waves in viscous matter. The amplitude of tidal distortion is not maximum when the bodies are closest to each other, but at a later time. We consider a fluid case, which corresponds to ϵ_i=0.

In addition, we define the tangential perturbation accelerations:

(A.6)

for the star and

(A.7)

for the planet. Finally, ℝ = * + R_p and S = S* + S_p.

There is a third perturbing acceleration component beyond 𝕊 and ℝ. It is normal to the plane of the orbit and positive toward the north pole. However, its contribution is zero as it is proportional to the latitude component of the angle of tidal lag. It can arise only if the orbital plane is inclined with respect to the equatorial plane of the body (star or planet). We neglect this term because the Rossiter-McLaughlin measurements of WASP-19A show an equatorial orbit (see, e.g., Sedaghati et al. 2021).

With the introduction of the terms, 𝕊 and ℝ, the perturbation equations can be easily written following Kopal (1959, 1978). The temporal change of the argument of periastron measured from the ascending node is:

(A.8)

where n is the mean motion that can be defined through Kepler’s law or through the orbital period P_orb.

Figure A.2 shows the effect of the periodical and secular perturbations on ω, as in Eq. (A.8). The plots were made setting the physical and orbital parameters to those of the exoplanetary system WASP-19Ab (see Sect. 3) with the Love number of WASP-18Ab (Csizmadia et al. 2019) as a first guess. We present how it evolves over 1, 10, 50, 100, and 500 orbits.

Fig. A.2 Temporal change of the osculating orbital element ω(t) due to tidal and rotational potentials, as described by Eq. (A.8), assuming the orbital properties of the system WASP-19 (see Sect. 3.1) and the Love number of WASP-18Ab Csizmadia et al. (2019).

Appendix B Planetary eccentricity driven by WASP-19B

Here, we show that a companion star on a very eccentric orbit is able to cause an eccentric orbit of the planet in the WASP-19 system.

Borkovits et al. (2011) gave analytic expressions for the variations of the orbital elements up to the sixth order of the inner eccentricity for any arbitrary orientation of the orbits. However, we do not know the spatial orientation of the orbit of the companion star B nor the longitude of the node of the planet b. Therefore, we assume co-planar orbits for the planet and star B, and we also assume that their semi-major axes are aligned. In this specific case, the integration of Eq. (9) of Borkovits et al. (2011) gives the simple expression for the eccentricity of the planet.

(B.1)

where f (υ_B) is

(B.2)

and, A_L,

(B.3)

Here υ_B is the true anomaly of star B and it holds the time-dependence, P and M are the orbital periods around star A and the masses of the bodies and e is the eccentricity. Indices A and B stand for stars A and B respectively, while index b denotes the planet.

The lower period limit of star B around star A can be estimated assuming that P_B is minimum if star B is in apoastron. In this case, the observed separation between stars A and B corresponds to a distance equal to a_B(1 + e_B). For different values of e_B, we can determine the semi-major axis that yields the orbital period, P_B, via Kepler’s third law. In case of e_B = 0.995 a_B ≈ l0⁴ AU and P_B ≈ 1 Myrs (note: the periastron point is still located at 50 AU from star A).

Since we do not know the exact orbit of star B, we just estimate the planetary eccentricity. Taking the Taylor series of the exponential function in Eq. (B.1) and neglecting the time dependence, we obtain

(B.4)

Using the aforementioned numerical values and assuming M_B~ 0.4M_⊙, we get e_b ≈ 0.0011. This calculation leads to an order of magnitude estimation of whether any (or a specific) combination of orbital elements can be responsible for the observed eccentricity. A deviation from this estimate can be due to the unknown spatial orientation of the different orbits.

Appendix C Tables of transits and occultations

Appendix D Tables of RVs

References

All Tables

All Figures

	Fig. 1 Radial velocity observations of WASP-19Ab, time vs. radial velocities, with corresponding observational errorbars. The different data sets (see Sect. 4) are represented by different colours. Between the datasets, there are significant offsets between different instruments and studies. They have been removed accordingly to the results in Table 4.
In the text

	Fig. 2 l1 periodogram of all RV data outside the transit (top). There is only one significant peak with false alarm probability ≃ 10−31. The periodogram of the residuals after the subtraction of the contribution of planet b (bottom). All residual peaks have false alarm probabilities of >20%.
In the text

Fig. 3 Properties of third bodies excluded by RV data. The plot shows the results for the simulated scenarios including a putative companion planet c of different masses at different periods. Each black dot represents one of the around 20 000 selected scenarios (see the text for more details). Top plot: χ2 value vs the orbital period of the perturber. The red lines correspond to the period of the stable orbits predicted by Cortés-Zuleta et al. (2020). The two horizontal dashed lines represent the χ2 required to have weak (upper line) or strong (lower line) evidence of a second planet, as discussed in the text. Bottom plot: mass and orbital period of the selected candidates.

In the text

Fig. 4 Star-planet system geometry when tidal interactions occur. The spherical shapes that star and the planet would have in the absence of tidal interactions are drawn in red and dark brown. The tidal bulges are drawn in orange and light brown, respectively. They arise on the stellar and planetary surfaces (respectively in orange and light brown) and rotate with the star itself and are elongated in the direction of the other body. The tidal lags are also drawn. The image is not to scale.

In the text

Fig. 5 χ2 values of the fit when varying in the model. The horizontal dashed line separates the statistically worse models (over the line) from the best models (below the dashed lines). The models below the dashed lines are statistically equivalent within 1σ (68% confidence level for the red line) and 2σ (95% confidence level for the green line). Note: the significantly larger χ2 at . The high amount of local minima is an indication of the complexity of the problem.

In the text

Fig. 6 Posterior distributions and correlation between the anomalistic period and periastron precession rate where the chains converged. The middle plot shows the strong correlation between these two free parameters, as in Eq. (12). On the sides, are the corresponding histograms of the posterior distributions. As can be seen from the histograms on the sides, many chains converged to various positive and negative peaks, as predicted by the semi-grid method and reported in Table 5.

In the text

	Fig. 7 Mass and orbital period of a second body that would cause a negative periastron precession rate of approximately −0.0638°d−1 (corresponding to the negative solution of the DE-MCMC analysis) on planet b, as deduced from Borkovits et al. (2011). See text for more details.
In the text

	Fig. 8 Residuals of the periastron precession fit for transits and occultation mid-times, calculated with the sidereal period. Transit mid-times are plotted as blue circles, occultation mid-times as red circles. The errorbars include the jitters, as in Table 4.
In the text

Fig. 9 Data with errorbars, phase-folded with our fit including periastron precession, shown in the upper plot. The colours represent different datasets, as in Fig. 1 (see also Sect. 4.2 for more details). Residuals of the periastron precession rate fit are shown in the lower plot. The error-bars include the jitters, as in Table 4. In both plots, the in-transit points are not shown as they were not used in the fit.

In the text

	Fig. A.1 Geometry of the orbital system in the simple case of orbital inclination of i =90°. The star symbol and the open dot denote the positions of the star and the planet, respectively; ω is the angle of periastron, θ is the phase from superior conjunction, and v is the true anomaly; l.o.s. denotes the line of sight. In this configuration, θ equals 0 for transits and π for occultations.
In the text

	Fig. A.2 Temporal change of the osculating orbital element ω(t) due to tidal and rotational potentials, as described by Eq. (A.8), assuming the orbital properties of the system WASP-19 (see Sect. 3.1) and the Love number of WASP-18Ab Csizmadia et al. (2019).
In the text