Dynamical masses of two infant giant planets


 Current theories of planetary evolution predict that infant giant planets have large radii and very low densities before they slowly contract to reach their final size after about several hundred million years 1, 2. These theoretical expectations remain untested to date, despite the increasing number of exoplanetary discoveries, as the detection and characterisation of very young planets is extremely challenging due to the intense stellar activity of their host stars 3, 4. However, the recent discoveries of young planetary transiting systems allow to place initial constraints on evolutionary models5–9. With an estimated age of 20 million years, V1298 Tau is one of the youngest solar-type stars known to host transiting planets: it harbours a multiple system composed of two Neptune-sized, one Saturn-sized, and one Jupiter-sized planets 10, 11. Here we report the dynamical masses of two of the four planets. We find that planet b, with an orbital period of 24 days, has a mass of 0.60 Jupiter masses and a density similar to the giant planets of the Solar System and other known giant exoplanets with significantly older ages 12, 13. Planet e, with an orbital period of 40 days, has a mass of 1.21 Jupiter masses and a density larger than most giant exoplanets. This is unexpected for planets at such a young age and suggests that some giant planets might evolve and contract faster than anticipated, thus challenging current models of planetary evolution.

and suggests that some giant planets might evolve and contract faster than anticipated, thus challenging current models of planetary evolution.
V1298 Tau is a very young K1 dwarf star with a mass of 1.170±0.060 M , a radius of 1.278±0.070 R , an effective temperature of 5050 ± 100 K, and solar metallicity (see Table 1). It is the physical companion of the G2 star HD 284154. Based on their position in color-magnitude diagrams, isochrone fitting, rotation period, lithium abundance, and X-ray emission, the pair belongs to the Group 29 stellar association 14 and has an age of 20 ± 10 Myr. and rapid rotation of the star. However, the different wavelength coverage of each instrument and the varying response of V1298 Tau's activity in the blue and red wavelengths -due to the different contrast of active regions at different wavelengths -introduced an additional complexity to the analysis of the RVs 24 . We approached this issue by allowing different amplitudes in the stellar model for the instruments with different spectral ranges.
Young stars show large variations in flux and in RV due to magnetic activity. In V1298 Tau these variations are much larger than those expected for any of the planetary-induced signals. These variations have also a higher frequency than those induced by the planetary orbits due to the relatively fast stellar rotation. To better monitor the changes caused by stellar activity, we performed contemporaneous V -band photometry with a cadence of 1 observation every 8 hours using the Las Cumbres Observatory (LCOGT) network 25 . We obtained 250 measurements with a typical precision of 10 parts per thousand (ppt). The data showed variations of up to 60 ppt, almost twice as large as those found in the K2 data. This is not surprising for a spot-dominated photosphere, as the Kepler passband is more extended towards red wavelengths, where the spot contrast is smaller.
During the 2019-2020 campaign, our photometric data showed a rotation-induced signal at a period of 2.9104 ± 0.0019 days, slightly different from the period measured during the 2015 K2 campaign (2.869 ± 0.013 days). This difference could be attributed to differential rotation.
To mitigate the effects of activity on the RV data of V1298 Tau, we opted for a global model combining the K2 photometry, RV, and one contemporaneous activity proxy 26 . For this particular purpose we used the LCOGT photometry as our proxy, as we could obtain a better observing cadence than for any of the spectroscopic indicators. The model uses the K2 photometry to constrain the planetary orbital periods, phases and radii, the LCOGT photometry and the RVs to constrain the timescales and amplitudes of the activity variations during our observing campaign, and the RVs to measure the masses of the planets. To model the activity components, we relied on a Gaussian processes (GP) regression 27,28 . Given the complexity of the activity variations, we tested several GP kernels. More details can be found in the supplementary material. To take advantage of the sampling of the photometric follow-up, we shared several hyperparameters between the RV and LCOGT time series. We added four Keplerian components to the RV model, representing the planets. The analysis considered a zero-point value and a noise term (jitter) for each dataset as free parameters to be optimised simultaneously, with the exception of the K2 data. For the K2 data we opted to manually include the white-noise component given in the discovery paper 29 . To sample the posterior distribution and obtain the statistical significance of the model we relied on a Markov chain Monte Carlo (MCMC) procedure, using emcee 30 .
We obtained significant measurements of the RV semi-amplitudes induced by planet b of 38 ± 12 m s −1 , and by planet e of 65 ± 16 m s −1 . For planets c and d, we can only set upper limits of 26 m s −1 and 24 m s −1 with a 99.7% confidence, respectively. Planet e was originally detected with a single transit, insufficient for the accurate measurement of its orbital period 11 . Using the RV data, we obtained the detection at a period of 40.2 ± 1.0 d, on the short end of the range expected from the transit duration 29 . If this period is correct, the K2 mission missed transits right before and after its campaign by just a few days. We measured the orbital eccentricities for planets b and e to be of 0.13 ± 0.07 and 0.12 ± 0.07 respectively. Figure 1 shows the phase folded curves of the RV signals attributed to planets b and e. The analysis of the RV and photometry during the 2019-2020 campaign yielded a stellar rotation period of 2.9104 ± 0.0019 days, with a semi-amplitude in RV of about 250 m s −1 , and of 5% of the flux in the V -band photometry. Our analysis of the K2 data yielded orbital periods, times of transit, and relative radii consistent with the discovery paper 10 .  Figure 2 shows the position of the planets orbiting V1298 Tau in a mass-radius diagram compared to the known population of exoplanets. Similar to the case of AU Mic b 9 -the only other exoplanet of similar age with a mass measurement -the mass-radius relation of the planets orbiting V1298 Tau resembles that of the planets of our Solar System and of the general population of known transiting exoplanets.
However, in contrast to the case of AU Mic b, the planets V1298 Tau b and e seem incompatible with the expected population derived from these models of evolution of planetary systems. Figure 3 shows the planets orbiting V1298 Tau with the expected population of exoplanets orbiting a 1 M star at the ages of 20 Myr and 5 Gyr (simulation NG76 35 using the Bern model 1 ) and with the massradius tracks coming from other different models 32, 33 . These planets should not reach this massradius configuration until hundreds of Myr later 2 . Considering their densities, it is not expected that the planets orbiting V1298 Tau will contract significantly in the future due to evaporation 36 .
Our result suggests that some giant planets might reach a mass-radius configuration compatible with the known mature population of exoplanets during their first 20 Myr of age.
An alternative explanation to the characteristics of V1298 Tau b and e could be offered by an extreme enrichment in heavy elements compared to giant planets of the Solar System and the general population of transiting exoplanets 37 . A fraction of heavy elements of 40-60% and 60-80% of the total mass of planets b and e would partially reconcile our results with the core-accretion evolutionary models 1, 32, 33 . These enrichment would correspond to 3-20 times the fraction of heavy elements of Jupiter for V1298 Tau b, and 5-25 times for V1298 Tau e. Such high metallicities are not expected for planets orbiting stars with solar metallicity. Figure 3 shows the planets orbiting   Competing Interests The authors declare that they have no competing financial interests.    It is located at the 3.5 m Zeiss telescope at the Centro Astronómico Hispano Alemán (Almería, Spain). We extracted the spectra with the CARACAL pipeline, based on flat-relative optimal extraction 51 . The wavelength calibration that was performed by combining hollow cathode lamps Effective temperature, surface gravity and metallicity of V1298 Tau We derived the stellar parameters and metallicity of V1298 Tau using the high-resolution HARPS-N spectra. The extracted, blaze-corrected 2D spectra were corrected for barycentric velocity (varying from +10.8 to −30.4 km s −1 ) and for RV (varying from +14.4 to +15.7 km s −1 ) and normalised to unity order by order with a third-order polynomial using our own IDL-based automated code 58 . All orders of all spectra were co-added and merged with a wavelength step of 0.01Å per pixel. The resulting 1D spectrum depicted in Figure ED 2 shows a signal-to-noise ratio of ∼ 107, 222, 284, 412 and 391 at 4200, 4800, 5400, 6000, and 6600Å, respectively.

HARPS-N and CARMENES RVs
Using the same automated code, we normalised and combined the HARPS-N spectra from our Third, we implemented a Bayesian python code that compares the observed spectrum with a synthetic spectrum in the spectral range 5350-5850Å (see Figure ED 1 Effective temperature, surface gravity and metallicity of HD 284154 HD 284154 is resolved as a double-lined spectroscopic binary (see Figure ED 3), which has been identified as a wide binary of V1298 Tau 57 . Using the FIES spectrum, we estimated a RV difference of δRV = 43.6 ± 1.0 km s −1 between the two stellar components, an identical stellar rotation of V rot = FIES spectrum cross-correlated with a mask of isolated and relatively strong Fe lines. Assuming both binary components have the same mass (see below), the center-of-mass RV is +14.8 ± 0.7 km s −1 , perfectly compatible with the center-of-mass RV of V1298 Tau.
We applied method #3 to derive the stellar parameters and metallicity of both stellar components of HD 284154. We assumed both stars have the same luminosity and origin and therefore the same stellar mass and metallicity. We then generated a grid of synthetic spectra in the parameter range T eff / log g/A(Fe) of 5500 − 6000/3.5 Masses, radii and luminosities We determined the bolometric luminosity of both V1298 Tau and HD 284154 by transforming the observed magnitudes into bolometric magnitudes using Gaia distances and colour-bolometric corrections 73 . We confirmed that the obtained values (shown in Table 1)  Uncertainties in the mass determination account for the temperature and luminosity uncertainties and also for the dispersion of the results from the different models including models with slightly different metallic composition. We obtain mass and radius estimates of 1.170 ± 0.060 M and  This suggests that the Group 29 and, hence, the V1298 Tau system, is older than the USco association (5-11 Myr 10 ) and has a similar age to Beta Pic (20 ± 10 Myr 79 , compatible with a stronger activity than the stars in the Pleiades 86, 87 , and has a UV excess (NUV−J = 8.15 ± 0.05 mag), characteristic of stars younger than ∼100 Myr 88 . In addition, employing the same Bayesian inference method used in previous section 78 , we also obtained estimates for the ages of V1298 Tau and HD 28415 of 9 ± 2 Myr and 13 ± 4 Myr, respectively. Using all previous results, we can constraints the age of V1298 Tau to be 20 ± 10 Myr.

Modelling
We fit simultaneously the K2 photometry, LCO photometry and RV time-series. We model the activity signals in RV and the LCO photometry using Gaussian Processes (GP) with celerite 89 .
We used a variation of the rotation kernel described in the original celerite article: where A represents the covariance amplitude, P rot is the rotation period, L 1 and L 2 represent the timescale of the coherence of the periodicity at the rotation period and its first harmonic, ∆ represents the scaling in amplitude of the variability at the first harmonic of the rotation period, and C the balance between the periodic and the non-periodic components. The equation also includes a term of uncorrelated noise (σ), independent for every instrument, added quadratically to the diagonal of the covariance matrix to account for all un-modelled noise components, such as uncorrected activity or instrumental instabilities. δ τ is the Kronecker delta function, and τ represents an inter- To sample the posterior distribution of the model we relied on emcee 30 , combining a differential evolution algorithm 92, 93 and a differential evolution algorithm with a snooker move 94 , with a probability of 80% and 20%, respectively, to be used in each iteration. We initialised a number of independent chains equal to four times the number of free parameters, and ran the sampler for 1 000 000 iterations. Every 10 000 iterations we computed the auto-correlation time of the posterior. We stopped earlier if the number of steps was larger than 50 times the auto-correlation timescale, and the timescale changed less than 1% from the previous measurement. The burn-in period is later defined as 20 times the auto-correlation timescale. Then we estimate the Bayesian evidence of the model (LnZ) using the posterior distribution 95 .
We significantly detected the signals corresponding to planets b and e and derived upper limits for the amplitudes of the signals corresponding to planets c and d. This was our most significant model, with a measured LnZ of -4477. To confirm our results we repeated the analysis described above using the combination of two simple harmonic oscillators (SHO) to model the RV and LCO variations. We obtained a similar result, with larger amplitude for planet b and smaller uncertainties, but less significant (LnZ of -4547). We also attempted to confirm the results using the Quasi-Periodic GP Kernel 28 to model the activity variations in the RV and LCO data, implemented using George 96 . Previous studies have found it effective to study young stars 97 . In this case we obtain lower amplitudes for the signals attributed to planets b and e, and higher for planet c. This was the least favoured of the models we tested (LnZ of -4553). For the most favoured model (PQP2) we tested the difference between having 4 planetary components in the RV, 2 planetary components (b and e) and no planetary components. We found that a model with 2 Keplerian components in the RV is 30 times more likely than a model with no planetary signals. The model with 4 Keplerian components is less significant than the model with 2 Keplerian components, which is not surprising considering we could not detect the RV signals of planets c and d.
As the results coming from the different GP models are slightly different, we performed simulations to test the accuracy of the amplitude measurements in this particular case. To do that we subtracted the detected planetary signals from the RV to create an "activity-only" dataset. Following the same procedure as with our original RV dataset, we tested that all the models recovered  To further test our results we opted for a different approach based on the correlation of the RVs with the photometric data. In spot-dominated stars, the activity-induced RV variations are correlated with the gradient of the flux 26 . This correlation can be used to detrend the data from stellar activity.
As we do not have simultaneous, but contemporanous, photometry, we calculated the gradient from the model of the photometry. Figure ED 10 shows how the RV data (from HARPS-N and CARMENES) compare to the gradient of the V-band flux. We modelled the rotation using a third order polynomial against the derivative of the flux. A first attempt left some residual power at the first harmonic of the rotation period, which led us to include a sinusoidal component at that period.
Our activity model is then defined as: where T0 was parametrised as JD 0 +P Rot · φ, with JD 0 = 2458791.627.
Using this model we recovered a very similar solution as with the mixture of two SHO kernels, although with much larger residuals. We detected the presence of the planets V1298 Tau b and e, and measured upper limits for the amplitudes of planets c and d. We measure the amplitude of planet e to be much larger than what was found with the GP models, which might be caused by the Keplerian model absorbing some unmodelled activity.
Table ED 1 shows the parameters used in the fit, the datasets involved in fitting every every parameter, the priors and the results obtained for the different models tested. Figure ED 7 shows the RV data with the best fit model for the planetary parameters, and a weighted average model for the activity component. Figure ED 8 shows the best fit to the contemporary photometry and Figure   ED 9 shows the best fit to the K2 observations.

Lessons learned and limitations
We found that the signal phase folded to the rotation period shows clearly the two modes of oscillation that our favoured GP Kernel describes. The amplitude of the rotation signal is 8 times larger than the amplitude measured for the signal related to V1298 Tau b, and 5 times larger than the signal related to planet e. Figure ED 11 shows the huge RV amplitude difference between the stellar variability and the planetary detected signals. In young exoplanets field the stellar activity signals engulf those signals related to the planets and therefore similarly large observational efforts with precise RV measurements will be required.
We found that not all GP Kernels behaved the same at all timescales in our dataset. The classic QP Kernel handles short-period signals quite well. However for longer period signals it seems to absorb a significant part of the Keplerian components, causing a clear underestimation of the measured amplitudes. The mixture of SHO kernels had the opposite behaviour. It underfits the activity component, leaving larger residuals and causing an overestimation of (some of) the Keplerian amplitudes. We found that our Kernel of choice (PQP2) provides better description of the activity variations of V1298 Tau and a more accurate determination of the Keplerian amplitudes.
It is important to remain cautious about the mass determined for the planet V1298 Tau e. The original detection did not constrain the orbital period, which is derived purely from the RV information. We