Issue 
A&A
Volume 684, April 2024



Article Number  A78  
Number of page(s)  20  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202346852  
Published online  05 April 2024 
Evidence of apsidal motion and a possible comoving companion star detected in the WASP19 system
^{1}
Institute of Planetary Research, German Aerospace Center (DLR),
Rutherfordstrasse 2,
12489
Berlin, Germany
email: lia.bernabo@dlr.de
^{2}
ELKHSZTE Stellar Astrophysics Research Group,
6500
Baja,
Szegedi út Kt. 766, Hungary
^{3}
Zentrum für Astronomie und Astrophysik, Technische Universität Berlin,
Hardenbergstr. 36,
10623
Berlin, Germany
^{4}
Institut für Geologische Wissenschaften, Freie Universität Berlin,
12249
Berlin, Germany
^{5}
Thüringer Landessternwarte Tautenburg,
Sternwarte 5,
07778
Tautenburg, Germany
^{6}
Dipartimento di Fisica, Università degli Studi di Torino,
Via Pietro Giuria, 1,
10125
Torino, Italy
Received:
9
May
2023
Accepted:
8
January
2024
Context. Love numbers measure the reaction of a celestial body to perturbing forces, such as the centrifugal force caused by rotation, or tidal forces resulting from the interaction with a companion body. These parameters are related to the interior density profile. The nonpoint mass nature of the host star and a planet orbiting around each other contributes to the periastron precession. The rate of this precession is characterized mainly by the secondorder Love number, which offers an opportunity to determine its value. When it is known, the planetary interior structure can be studied with one additional constraint beyond the mass, radius, and orbital parameters.
Aims. We aim to redetermine the orbital period, eccentricity, and argument of the periastron for WASP19Ab, along with a study of its periastron precession rate. We calculated the planetary Love number from the observed periastron precession rate, based on the assumption of the stellar Love number from stellar evolutionary models.
Methods. We collected all available radial velocity (RV) data, along with the transit and occultation times from the previous investigations of the system. We supplemented the data set with 19 new RV data points of the host star WASP19A obtained by HARPS. Here, we summarize the technique for modeling the RV observations and the photometric transit timing variations (TTVs) to determine the rate of periastron precession in this system for the first time.
Results. We excluded the presence of a second possible planet up to a period of ~4200 d and with a radial velocity amplitude bigger than ≃ 1 m s^{−1}. We show that a constant period is not able to reproduce the observed radial velocities. We also investigated and excluded the possibility of tidal decay and longterm acceleration in the system. However, the inclusion of a small periastron precession term did indeed improve the quality of the fit. We measured the periastron precession rate to be 233_{−35}^{+25}″d^{−1}. By assuming synchronous rotation for the planet, it indicates a k_{2} Love number of 0.20_{−0.03}^{+0.02} for WASP19Ab.
Conclusions. The derived k_{2,p} value of the planet has the same order of magnitude as the estimated fluid Love number of other Jupitersized exoplanets (WASP18Ab, WASP103b, and WASP121b). A low value of k_{2,p} indicates a higher concentration of mass toward the planetary nucleus.
Key words: planets and satellites: dynamical evolution and stability / planets and satellites: interiors / planetstar interactions
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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 massradiuscomposition 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 redistribution 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 secondorder 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 midtimes. We apply it to the system WASP19 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 seconddegree Love numbers.
This paper is organized as follows. Section 2 and Appendix A present the theory of apsidal motion. Section 3 describes the system WASP19 and provides evidence of a possible companion star in the WASP19 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:
The general relativistic term dω_{GR}/dt (Einstein 1915) is:
where G is the gravitational constant, M_{⋆} the mass of the star, a and e the semimajor 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):
and
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 secondorder fluid Love number corresponding to the double the apsidal motion constant (Csizmadia et al. 2019):
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 WASP19Ab 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^{−1} ≃ 0.0041–0.479° d^{−1}, 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 longterm 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 WASP19Ab), 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 semimajor axis (tidal decay), eccentricity (circularization) and time of periastron passage. The semimajor axis and its eccentricity would diminish in value and the time of periastron passage is shifted in time.
Expected gravitational, rotational and tidal contribution to the periastron precession in the extreme Love number cases (masspoint and homogeneous planet) and the corresponding apsidal motion period, assuming synchronous rotation for the planet.
WASP19A and its stellar companion
3 WASP19Ab
3.1 The system
The planet WASP19Ab, hosted by the bright, active Gdwarf 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 WASP19Ab 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 midtransit time. Mancini et al. (2013) suggested a nonlinear 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ésZuleta 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: WASP19B
We searched the catalogue of the third Gaia data release (DR3) (Gaia Collaboration 2016, 2023) for possible comoving companions to WASP19A. We searched all objects within 4′ of the target, which corresponds to a physical separation of around 0.3 pc at the distance of WASP19A, 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 Gband, has a parallax that differs from that of WASP19A 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 southeast region of WASP19A, at a separation of 68″. The magnitude and BPRP colour of the companion suggests that this star is a mid Mdwarf, 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 WASP19 in parallax and proper motions, we searched a much larger area of sky in Gaia DR3, still centered on WASP19, but with a radius of 2°. This produced 77 objects (including the companion; hereafter: ‘WASP19B’) that match WASP19 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 WASP19B. 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 WASP19A than WASP19B, demonstrate as good a match as WASP19B. Mathematically, this is expressed as WASP19B 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 WASP19 and its companion, respectively. We recognise that the uncertainty on the parallax of WASP19B is relatively large. Even if WASP19B is likely to be a true, gravitationally bound, companion to WASP19A, some doubt remains as to whether it truly is a bound companion to WASP19A. We adopted the planet name WASP19Ab because it orbits around the brighter primary star. Details of WASP19B 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
Midtransit 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)^{2} project database, part of the Exoplanet Transit Database (Poddaný et al. 2010). In particular, we included amateur data that were selected by Mancini et al. (2013) and other highquality data showing a clear light curve without gaps. As detailed in the tables, we also reanalyzed some transits with the Transit and Light Curve Modeller (TLCM, Csizmadia 2020).
The time unit of each midtransit and midoccultation 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 midtimes 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.
Radial velocity datasets used in this work.
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 RossiterMclaughlin 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 RossiterMcLaughlin 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 RossiterMcLaughlin 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.C0849 (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 signaltonoise 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 midtimes, the time unit of the observation was transformed and uniformized to BJD_{TDB} when necessary.
Fig. 1 Radial velocity observations of WASP19Ab, 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 Longterm effects in RVs
In this section, we study the longterm 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_{1} 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 resampling approach suggested by Hara et al. (2017): we randomly removed 10–20% of the data and recomputed the periodogram. Some peaks did not appear after the resampling, therefore, we did not consider them real. After the resampling, 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 WASP19Ab, with a false alarm probability (FAP) of log_{10}(FAP) = −30.9 (FAP ≃ 10^{−31}). The analysis shows no evidence of any additional periodic signal in the RV data.
We also recomputed the periodogram on the RV residuals after the subtraction of the signal of the first planet (see bottom plot in Fig. 2). After a resampling 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^{−1}. 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 semiamplitude K_{c} between 0.1 and 100 m s^{−1} (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:
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 χ^{2} ⩽ 500.The top plot shows the χ^{2} of the model when including the proposed perturber, corresponding to five more free parameters: transit epoch, period, RV semiamplitude, and (plotted as a function of the orbital period). Many of these candidates are discarded because the corresponding χ^{2} is too high or the amplitude of the perturbing body on the RV curve would be over a dozen m s^{−1}, 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ésZuleta 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ésZuleta et al. (2020) with red dashed lines.
The expected χ^{2} would be around the number of degrees of freedom (19 total datapoints for dataset 8, of which 15 were fitted as outoftransit points and 11 free parameters), while the lowest χ^{2} 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 oneplanet eccentric orbit model to compare it to the scenarios with a perturbed. Without adding any RV jitter, the χ^{2} of this model is ≃35. Since the twoplanet 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 twoplanet model is preferred, a χ^{2} 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 χ^{2}, resulting in an even more unlikely presence of a perturbed. Our results exclude it up to an RV amplitude of ≃1 m s^{−1}.
This simulation confirms the conclusion of Knutson et al. (2014) when looking for the acceleration of the system due to possible longperiod 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 longterm trends in the data, we proceeded with our analysis by introducing periastron precession in the fit.
Fig. 2 l_{1} 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ésZuleta 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 nonzero 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 midtransit and occultation times. Since then, the code has been modified and updated. Now the Genetic Algorithm minimization is followed by a differentialevolution Markov chain Monte Carlo (DEMCMC, 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 DEMCMC analysis to derive the posterior probability distributions. The Genetic Algorithm is initialized with a population set made of 200 individuals. In the DEMCMC we run 90 chains. The first 6 × 10^{4} steps of each chain are discarded as part of the burnin phase. They are followed by 10^{4} more DEMCMC steps, which are used for the posterior distributions. The convergence is checked by the Gelman  Rubin test applied on the steps that follow the burnin 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 χ^{2} 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 timevariable ω:
with a secular advance of the longitude of periastron, , caused by the tidal and rotational potential (as expressed in Eq. (A.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 WASP19Ab system relative to the barycenter of the Solar System. Different datasets may have different instrumental zeropoint offsets, V_{i,instr}. Finally, K is the radial velocity halfamplitude:
In this expression, a_{⋆} is the semimajor 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 nonorbital 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 reemitted (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 nonorbital contribution to the RV curve:
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_{2} takes into account the limb darkening. For more details on f_{2} see Kopal (1959) and Csizmadia et al. (2019).
Arras et al. (2012) already includes such calculations for precise RV analysis of WASP18A, but neglects the effect of stellar rotation. In other words, a synchronous rotation of the star is assumed. However, for WASP19A 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^{−1}, so it is negligible for WASP19A if compared to the typical observational errorbars. 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^{−1} for Hebb et al. (2010)’s data (ID 1, as in Anderson et al. 2010), 14.1 m s^{−1} and 6.9 m s^{−1} for Coralie (ID 2) and HARPS (ID 3) Hellier et al. (2011)’s data, respectively, 20 m s^{−1} for Albrecht et al. (2012)’s (ID 4) and 17.8 m s^{−1} for Knutson et al. (2014)’s (ID 5) RVs. Then, we refined all the jitter values by minimizing the χ^{2} as:
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 WASP19Ab system, the datasets come from various instruments with different precisions and, thus, are nonhomogeneous. 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 χ^{2} 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 midtransit and midoccultation 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 midtimes, we followed the approach of Gimenez & GarciaPelayo (1983) and Csizmadia et al. (2019).
Fig. 4 Starplanet 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. 
Parameters of the WASP19Ab system.
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 socalled semigrid analysis. This procedure also helps to get an initial range for the fitted value of . The limiting values were and + 0.45°d^{−1}. We also fixed the epoch to 2 454 775.3372 BJD_{TDB}, while all other parameters were free in the DEMCMC analysis. The result is presented in Fig. 5, showing the assumed value of the periastron precession rate and the corresponding χ^{2} 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 WASP19Ab, 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 semigrid method (Sect. 5.3), here we fit along with the other parameters.
During the analysis, we noticed that different DEMCMC 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}:
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 midtimes 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 semigrid method approach; see Fig. 5).
The highest peak of (solution B) in the posterior distribution results around approximately +0.065°d^{−1} (summarized in Table 5). The positive solutions C, D, and E can be discarded because of their high χ^{2} 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 RossiterMcLaughlin 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^{−1} = −229.68″d^{−1}, 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^{−1}) 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 longterm (over a 10day 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 χ^{2} are shown in Table 4, along with their errorbars. With 99 RV points (outside the transit, 203 in total), 162 midtransit and 11 midoccultation points, and 19 free parameters, the reduced χ^{2} is 1.02. The observed value of the total contribution from the DEMCMC 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 midtimes with such a fit. Figure 9 shows the phasefolded 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 nonrotating (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 uptodate values of the tidal Love number for planets and satellites in the Solar System and exoplanets are listed in Table 6. The value obtained for the fluid Love number of WASP19Ab 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 highmass 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 secondorder 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 smallshaped 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 midtransit times, which compromise the reliability of the results, increasing the χ^{2} 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 midtimes 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 semigrid 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 DEMCMC 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 midtimes, calculated with the sidereal period. Transit midtimes are plotted as blue circles, occultation midtimes as red circles. The errorbars include the jitters, as in Table 4. 
Fig. 9 Data with errorbars, phasefolded 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 errorbars include the jitters, as in Table 4. In both plots, the intransit points are not shown as they were not used in the fit. 
Fluid and tidal Love numbers measured for Solar System planets and minor bodies and exoplanets.
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 longterm acceleration of the system was fitted by adding a γ acceleration to the RV Eq. (7) as:
We ran the DEMCMC algorithm in the case of a circular and an eccentric orbit. The posterior distribution of has median of (7.37 ± 5.02) 10^{−7}m s^{−2} for the circular case and (7.43 ± 4.89) 10^{−7} m s^{−2} 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^{−7} m s^{−2}. Our results increase the confidence level of a nondetection of longterm 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 WASP12b (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 semiamplitude K in Eq. (9) as a timevarying parameter:
Deriving Kepler’s laws with respect to time, we can see that a change in the RV semiamplitude results in an opposite change in the semimajor axis.
For the circular orbit case, we obtained while for the eccentric case (−2.60 ± 2.73) × 10^{−8} m s^{−2}. As Rosário et al. (2022), we did not find significant evidence of orbital decay. Both in the case of a longterm 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 wellestablished 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 WASP19Ab, by using archival midtransit 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 longterm 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 WASP19Ab 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 WASP19Ab is on the same order of magnitude as the one of other Jupiterlike planets (WASP4b, WASP18Ab, WASP103b and WASP121b). A low value suggests the presence of a dense core in the planet, similar to the case of Saturn, Neptune, and WASP121b. Note that the Love numbers of WASP12b and WASP4b (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, longperiod 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 WASP19A 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/141 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 WASP19Ab 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.C0849. 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:
and
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 secondorder fluid Love number. Here, r is the radius vector:
where a is the semimajor 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 asynchronism 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:
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 viceversa (see Kaula (1964) and reference therein).
The radius vector, r_{ϵ}, is defined as:
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:
for the star and
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 RossiterMcLaughlin measurements of WASP19A 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:
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 WASP19Ab (see Sect. 3) with the Love number of WASP18Ab (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 WASP19 (see Sect. 3.1) and the Love number of WASP18Ab Csizmadia et al. (2019). 
Appendix B Planetary eccentricity driven by WASP19B
Here, we show that a companion star on a very eccentric orbit is able to cause an eccentric orbit of the planet in the WASP19 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 coplanar orbits for the planet and star B, and we also assume that their semimajor 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.
where f (υ_{B}) is
and, A_{L},
Here υ_{B} is the true anomaly of star B and it holds the timedependence, 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 semimajor axis that yields the orbital period, P_{B}, via Kepler’s third law. In case of e_{B} = 0.995 a_{B} ≈ l0^{4} 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
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
Transits found in the literature: cycle number, referred to in Hebb et al. (2010)  midtransit point in BID_{TDB}, its error, details on the observations, and reference.
Occultations found in the literature.
Appendix D Tables of RVs
Radial velocity of WASP19Ab system.
References
 Abe, L., Gonçalves, I., Agabi, A., et al. 2013, A&A, 553, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Albrecht, S., Winn, J. N., Johnson, J. A., et al. 2012, ApJ, 757, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Anderson, D. R., Gillon, M., Maxted, P. F. L., et al. 2010, A&A, 513, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Anderson, D. R., Smith, A. M. S., Madhusudhan, N., et al. 2013, MNRAS, 430, 3422 [NASA ADS] [CrossRef] [Google Scholar]
 Arras, P., Burkart, J., Quataert, E., & Weinberg, N. N. 2012, MNRAS, 422, 1761 [Google Scholar]
 Baluev, R. V. 2009, MNRAS, 393, 969 [Google Scholar]
 Banerdt, W. B., Smrekar, S. E., Banfield, D., et al. 2020, Nat. Geosci., 13, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Barros, S. C. C., Akinsanmi, B., Boué, G., et al. 2022, A&A, 657, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Batygin, K., Bodenheimer, P., & Laughlin, G. 2009, ApJ, 704, L49 [NASA ADS] [CrossRef] [Google Scholar]
 Baumeister, P., Padovan, S., Tosi, N., et al. 2020, ApJ, 889, 42 [Google Scholar]
 Bean, J. L., Désert, J.M., Seifahrt, A., et al. 2013, ApJ, 771, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, J. C., & Batygin, K. 2013, ApJ, 778, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Borkovits, T., Csizmadia, S., ForgácsDajka, E., & Hegedüs, T. 2011, A&A, 528, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boss, A. P. 1997, Science, 276, 1836 [Google Scholar]
 Bouma, L. G., Winn, J. N., Baxter, C., et al. 2019, AJ, 157, 217 [NASA ADS] [CrossRef] [Google Scholar]
 Bouma, L. G., Winn, J. N., Howard, A. W., et al. 2020, ApJ, 893, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Buhler, P. B., Knutson, H. A., Batygin, K., et al. 2016, ApJ, 821, 26 [CrossRef] [Google Scholar]
 Burton, J. R., Watson, C. A., Littlefair, S. P., et al. 2012, ApJS, 201, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Campo, C. J., Harrington, J., Hardy, R. A., et al. 2011, ApJ, 727, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Claret, A. 2004, A&A, 424, 919 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 CortésZuleta, P., Rojo, P., Wang, S., et al. 2020, A&A, 636, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Csizmadia, S. 2020, MNRAS, 496, 4442 [Google Scholar]
 Csizmadia, S., Hellard, H., & Smith, A. M. S. 2019, A&A, 623, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Damasso, M., Bonomo, A. S., AstudilloDefru, N., et al. 2018, A&A, 615, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dorn, C., Khan, A., Heng, K., et al. 2015, A&A, 577, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dragomir, D., Kane, S. R., Pilyavsky, G., et al. 2011, AJ, 142, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Durante, D., Parisi, M., Serra, D., et al. 2020, Geophys. Res. Lett., 47, e86572 [NASA ADS] [CrossRef] [Google Scholar]
 Eastman, J., Siverd, R., & Gaudi, B. S. 2010, PASP, 122, 935 [Google Scholar]
 Ebr, J., Janeček, P., Prouza, M., et al. 2014, Rev. Mex. Astron. Astrofis. Conf. Ser., 45, 114 [NASA ADS] [Google Scholar]
 Efroimsky, M., & Makarov, V. V. 2013, ApJ, 764, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Einstein, A. 1915, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin), 831 [Google Scholar]
 Espinoza, N., Rackham, B. V., Jordán, A., et al. 2019, MNRAS, 482, 2065 [Google Scholar]
 Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gelman, A., & Rubin, D. B. 1992, Statist Sci., 7, 457 [NASA ADS] [Google Scholar]
 Gibson, N. P., Aigrain, S., Pollacco, D. L., et al. 2010, MNRAS, 404, L114 [NASA ADS] [Google Scholar]
 Gimenez, A., & GarciaPelayo, J. M. 1983, Ap&SS, 92, 203 [NASA ADS] [CrossRef] [Google Scholar]
 Hara, N. C., Boué, G., Laskar, J., & Correia, A. C. M. 2017, MNRAS, 464, 1220 [Google Scholar]
 Hardy, R. A., Harrington, J., Hardin, M. R., et al. 2017, ApJ, 836, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Harre, J.V., & Smith, A. M. S. 2023, Universe, 9, 12, 506 [NASA ADS] [CrossRef] [Google Scholar]
 Hebb, L., CollierCameron, A., Triaud, A. H. M. J., et al. 2010, ApJ, 708, 224 [NASA ADS] [CrossRef] [Google Scholar]
 Hellard, H., Csizmadia, S., Padovan, S., et al. 2019, ApJ, 878, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Hellard, H., Csizmadia, S., Padovan, S., Sohl, F., & Rauer, H. 2020, ApJ, 889, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Hellier, C., Anderson, D. R., CollierCameron, A., et al. 2011, ApJ, 730, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Holman, M. J., & Murray, N. W. 2005, Science, 307, 1288 [Google Scholar]
 Huitson, C. M., Sing, D. K., Pont, F., et al. 2013, MNRAS, 434, 3252 [NASA ADS] [CrossRef] [Google Scholar]
 Iess, L., Jacobson, R. A., Ducci, M., et al. 2012, Science, 337, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Ivshina, E. S., & Winn, J. N. 2022, ApJS, 259, 62 [CrossRef] [Google Scholar]
 Jackson, B., Greenberg, R., & Barnes, R. 2008, ApJ, 678, 1396 [Google Scholar]
 Jacobson, R. A. 2009, AJ, 137, 4322 [Google Scholar]
 Jacobson, R. A., & Lainey, V. 2014, Planet. Space Sci., 102, 35 [CrossRef] [Google Scholar]
 Kass, R. E., & Raftery, A. E. 1995, J. Am. Statist. Assoc., 90, 773 [CrossRef] [Google Scholar]
 Kaula, W. M. 1964, Rev. Geophys. Space Phys., 2, 661 [Google Scholar]
 Knutson, H. A., Fulton, B. J., Montet, B. T., et al. 2014, ApJ, 785, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Konopliv, A. S., & Yoder, C. F. 1996, Geophys. Res. Lett., 23, 1857 [CrossRef] [Google Scholar]
 Konopliv, A. S., Park, R. S., & Ermakov, A. I. 2020, Icarus, 335, 113386 [NASA ADS] [CrossRef] [Google Scholar]
 Kopal, Z. 1959, Close Binary Systems (The International Astrophysics Series, London: Chapman & Hall) [Google Scholar]
 Kopal, Z. 1978, Dynamics of Close Binary Systems (Astrophysics and Space Science Library, Dordrecht: Reidel) [Google Scholar]
 Lainey, V., Jacobson, R. A., Tajeddine, R., et al. 2017, Icarus, 281, 286 [Google Scholar]
 Lendl, M., Gillon, M., Queloz, D., et al. 2013, A&A, 552, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Love, A. E. H. 1911, Some Problems of Geodynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Mancini, L., Ciceri, S., Chen, G., et al. 2013, MNRAS, 436, 2 [Google Scholar]
 Mandell, A. M., Haynes, K., Sinukoff, E., et al. 2013, ApJ, 779, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20 [NASA ADS] [Google Scholar]
 Munk, W. H., & MacDonald, G. J. F. 1960, The Rotation of the Earth; a Geophysical Discussion (Cambridge: Cambridge University Press) [Google Scholar]
 Patra, K. C., Winn, J. N., Holman, M. J., et al. 2017, AJ, 154, 4 [Google Scholar]
 Patra, K. C., Winn, J. N., Holman, M. J., et al. 2020, AJ, 159, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9 [Google Scholar]
 Petrucci, R., Jofré, E., Gómez Maqueo Chew, Y., et al. 2020, MNRAS, 491, 1243 [NASA ADS] [Google Scholar]
 Poddaný, S., Brát, L., & Pejcha, O. 2010, New A, 15, 297 [CrossRef] [Google Scholar]
 Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Ragozzine, D., & Wolf, A. S. 2009, ApJ, 698, 1778 [Google Scholar]
 Rosário, N. M., Barros, S. C. C., Demangeon, O. D. S., & Santos, N. C. 2022, A&A, 668, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Russell, H. N. 1928, MNRAS, 88, 641 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, J. 2018, in Handbook of Exoplanets, eds. H. J. Deeg, & J. A. Belmonte (Berlin: Springer), 163 [Google Scholar]
 Sedaghati, E., Boffin, H. M. J., Csizmadia, S., et al. 2015, A&A, 576, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sedaghati, E., Boffin, H. M. J., MacDonald, R. J., et al. 2017, Nature, 549, 238 [Google Scholar]
 Sedaghati, E., MacDonald, R. J., CasasayasBarris, N., et al. 2021, MNRAS, 505, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Sterne, T. E. 1939, MNRAS, 99, 451 [NASA ADS] [Google Scholar]
 Taubner, R. S., Leitner, J. J., Firneis, M. G., & Hitzenberger, R. 2016, Origins Life Evol. Biosph., 46, 283 [NASA ADS] [CrossRef] [Google Scholar]
 Ter Braak, C. J. F. 2006, Statist. Comput., 16, 239 [Google Scholar]
 Torres, G., Fischer, D. A., Sozzetti, A., et al. 2012, ApJ, 757, 161 [Google Scholar]
 TregloanReed, J., Southworth, J., & Tappert, C. 2013, MNRAS, 428, 3671 [Google Scholar]
 Turner, J. D., RiddenHarper, A., & Jayawardhana, R. 2021, AJ, 161, 72 [NASA ADS] [CrossRef] [Google Scholar]
 Turner, J. D., Flagg, L., RiddenHarper, A., & Jayawardhana, R. 2022, AJ, 163, 281 [NASA ADS] [CrossRef] [Google Scholar]
 Valencia, D., Sasselov, D. D., & O'Connell, R. J. 2007, ApJ, 665, 1413 [NASA ADS] [CrossRef] [Google Scholar]
 Van Hoolst, T., Noack, L., & Rivoldini, A. 2019, Adv. Phys. X, 4, 1630316 [NASA ADS] [Google Scholar]
 Wagner, F. W., Sohl, F., Hussmann, H., Grott, M., & Rauer, H. 2011, Icarus, 214, 366 [NASA ADS] [CrossRef] [Google Scholar]
 Williams, J. G., Konopliv, A. S., Boggs, D. H., et al. 2014, J. Geophys. Res. Planets, 119, 1546 [NASA ADS] [CrossRef] [Google Scholar]
 Wong, I., Knutson, H. A., Kataria, T., et al. 2016, ApJ, 823, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Wong, I., Benneke, B., Shporer, A., et al. 2020, AJ, 159, 104 [Google Scholar]
 Wong, I., Shporer, A., Vissapragada, S., et al. 2022, AJ, 163, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Yee, S. W., Winn, J. N., Knutson, H. A., et al. 2020, ApJ, 888, L5 [Google Scholar]
 Zacharias, N., Finch, C. T., Girard, T. M., et al. 2012, VizieR Online Data Catalog: I/322A [Google Scholar]
 Zhou, G., Bayliss, D. D. R., KedzioraChudczer, L., et al. 2014, MNRAS, 445, 2746 [Google Scholar]
All Tables
Expected gravitational, rotational and tidal contribution to the periastron precession in the extreme Love number cases (masspoint and homogeneous planet) and the corresponding apsidal motion period, assuming synchronous rotation for the planet.
Fluid and tidal Love numbers measured for Solar System planets and minor bodies and exoplanets.
Transits found in the literature: cycle number, referred to in Hebb et al. (2010)  midtransit point in BID_{TDB}, its error, details on the observations, and reference.
All Figures
Fig. 1 Radial velocity observations of WASP19Ab, 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 l_{1} 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ésZuleta 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 Starplanet 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 semigrid 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 DEMCMC 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 midtimes, calculated with the sidereal period. Transit midtimes are plotted as blue circles, occultation midtimes as red circles. The errorbars include the jitters, as in Table 4. 

In the text 
Fig. 9 Data with errorbars, phasefolded 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 errorbars include the jitters, as in Table 4. In both plots, the intransit 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 WASP19 (see Sect. 3.1) and the Love number of WASP18Ab Csizmadia et al. (2019). 

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.