© The Authors 2023

1 Introduction

Despite the discovery of thousands of exoplanets, the mechanisms of planet formation and evolution have not been well tested with observations, mainly because of the lack of a complete photometric and spectroscopic characterization of young planets. The determination of their physical parameters is a great challenge because the stellar activity of host young stars is high. The activity is produced by a rapid rotation and strong magnetic activity, which, in most cases, produce photometric and radial velocity variations that are several times greater than the signals assigned to the planets. This also hampers the search for new young planet candidates.

One of the critical questions of the last years in exoplanet research is the existence of a bimodal distribution in the size of the small planets, that is, the so-called radius valley or Fulton gap (Fulton et al. 2017). Several evolution models have been proposed to explain this gap by considering that a single population of rocky planets (with a similar composition to that of Earth) is formed with a gaseous envelope that is sculpted by evolution through rapid loss of its atmosphere (about several tens of million years) by photoevaporation (Owen & Wu 2017) or slowly (~1Gyr) by core-powered mass loss (Ginzburg et al. 2018). However, the recent study by Luque & Pallé (2022) shows that in M dwarfs, the planets around the valley are distributed based on their composition (i.e. densities) into three classes: rocky, water-rich (i.e. planets that are 50% rocky and 50% water ice by mass), and gas-rich planets (i.e. could be rocky planets with massive H/He envelopes or water-rich planets with less massive envelopes). According to the models of Venturini et al. (2020), the formation of water and gaseous planets occurs beyond the ice line, and the planets later migrate inwards. The radius valley is explained in these models by the initial differences in the core sizes below or beyond ice lines and the subsequent loss of the envelopes of gas planets that migrated outside-in. To constrain the planet formation and evolution models, it is therefore critical to measure the density of planets with sizes between those of Earth and Neptune at different ages below 1 Gyr to determine whether most of the super-Earths are rocky planets from their origin, or if they are born as small gaseous planets that lose their atmospheres.

Recently, several planets with sizes between those of Earth and Neptune have been found orbiting stellar members of the Hyades (Mann et al. 2016a, 2018; Ciardi et al. 2018), the Praesepe (Obermeier et al. 2016; Mann et al. 2017; Rizzuto et al. 2018) and δ Lyr (Bouma et al. 2022b) open clusters, in the Ursa Major (Mann et al. 2020), β Pic (Plavchan et al. 2020), Pisces-Eridanus (Newton et al. 2021), Melange-1 (Tofflemire et al. 2021), and AB Doradus (Zhou et al. 2022) young moving groups, in the Cep-Her complex (Bouma et al. 2022a), and around other young field stars (David et al. 2018a,b; Sun et al. 2019; Zhou et al. 2021; Ment et al. 2021; Hedges et al. 2021; Kossakowski et al. 2021; Barragán et al. 2022; Vach et al. 2022) during the K2 and TESS missions. Only six of these young planetary systems have measured masses and accordingly, densities, derived from dedicated radial velocity (RV) campaigns: AUMicb and c (~20Myr; Klein et al. 2021, Klein et al. 2022; Cale et al. 2021; Zicher et al. 2022), TOI-1807b (~300 Myr; Nardiello et al. 2022), TOI-560b and c (~490 Myr; Barragán et al. 2022; El Mufti et al. 2023), K2-25 b (~725 Myr; Stefansson et al. 2020), K2-100b (~750 Myr; Barragán et al. 2019), and TOI-1201b (600–800 Myr; Kossakowski et al. 2021).

HD 63433 (TOI-1726, V377 Gem) is a bright (V = 6.9 mag) solar-type star member of the Ursa Major moving group (Montes et al. 2001) with an estimated age of 414 ± 23 Myr (Jones et al. 2015). Its two transiting exoplanets have orbital periods and radii of 7.11 days and 2.15 ± 0.10 R_⊕ for the inner planet, and 20.55 days and 2.67 ± 0.12 R_⊕ for the outer planet (Mann et al. 2020). Furthermore, by modelling the Rossiter-McLaughlin (RM) effect (Rossiter 1924; McLaughlin 1924), HD 63433 b and c have been measured a prograde orbits with a sky-projected obliquity of $\lambda ={1.0}_{-43.0}^{+41.0\circ }$ (Mann et al. 2020) and $\lambda =-{1.0}_{-32.0}^{+35.0\circ }$ (Dai et al. 2020), respectively. The atmospheric escape has also been studied for both planets, with the detection of Ly α absorption in HD 63433 c but not in HD 63433 b, whereas no helium absorption (10 833 Å) was detected in either planet (Zhang et al. 2022). These results seem to indicate that the two planets have different atmospheric compositions (Zhang et al. 2022).

In this paper, we present the mass determination of HD 63433 c and an upper limit for the mass of HD 63433 b. The paper is organised as follows. In Sect. 2, we describe the photometry of the Transiting Exoplanet Survey Satellite (TESS) and ground based follow-up observations of the system. In Sect. 3, we revise the physical properties of the star. In Sects. 4 and 5, we perform a transit and RV analysis of the planets, respectively. A joint-fit modelling of the photometric and RV time-series is carried out in Sect. 6. In Sect. 7, we discuss the composition of the planet and the main implications for theoretical models of exoplanets. We summarize our main results in Sect. 8.

Fig. 1 TESS TPF plots for HD 63433 (white crosses). The red squares indicate the best optimal photometric aperture used to obtain the SAP flux. G-band magnitudes from Gaia DR3 are shown with different sizes of red circles for all nearby stars with respect to HD 63433 up to 8 magnitudes fainter.

2 Observations

2.1 TESS photometry

HD 63433 was observed by TESS in five sectors with a 2-min short cadence. The sector used for the discovery (Mann et al. 2020), sector 20, ran from 24 December 2019 until 21 January 2020 as part of cycle 2 of the TESS primary mission. The other four sectors (44–47) were contiguous and provided a light curve between 12 October 2021 and 28 January 2022 during cycle 4 of the first extended mission. At the time of this publication, the target is not expected to be observed again with TESS.

All sectors were processed by the Science Processing Operations Center (SPOC; Jenkins et al. 2016) photometry and transit search pipeline at the NASA Ames Research Center. The light curves and target pixel files (TPFs) were downloaded from the Mikulski Archive for Space Telescopes¹ (MAST), which provides the simple aperture photometry (SAP) and the pre-search data conditioning SAP flux (PDCSAP). To confirm that the best photometric aperture for the TESS light curve does not include any additional bright sources down to 8 mag fainter and contaminates the light curve and the transit depth, we plot in Fig. 1 the TPF using tpfplotter² (Aller et al. 2020), which produces the best SAP flux. We also overplot the Gaia Data Release 3 (DR3) catalogue (Gaia Collaboration 2023) on top of the TPF. We searched for possible sources of contamination and verified that no source lay close to the chosen photometric aperture. For the rest of our analysis, we used the PDCSAP flux, which was corrected for instrumental errors and crowding. However, we tested that our analysis and results using SAP flux light curves are similar. The PDCSAP flux light curve for all sectors of HD 63433 is illustrated in Fig. 2 along with the best-fit model (see Sect. 6 for details). The light curve has a dispersion of σ_TESS ~ 4.3 parts per thousand (ppt), an average error bar of ~0.2 ppt, and a modulation with peak-to-peak variations up to 20 ppt whose amplitude changes significantly from one period to the next. Moreover, the light curve presents a few flares with amplitudes smaller than 1.5 ppt that were removed with the following procedure: first, we created a smoothed model using a Savitzky–Golay filter. After removing the model from the light curve, we iteratively removed all values higher than 3 times the root mean square of the residuals. We verified that most of the flares were removed through this process.

Fig. 2 Light curves of HD 63433 for the five TESS sectors. The blue dots correspond to the PDCSAP flux data. The black line and the grey shaded region indicate the best model and its 1σ uncertainty, respectively. The vertical lines show the times of the planetary transits for HD 63433 b (orange) and HD63433 c (purple).

2.2 LCO photometry

To monitor the photometric stellar activity contemporaneously with the RV, we observed HD 63433 using the 40 cm telescopes of Las Cumbres Observatory Global Telescope (Brown et al. 2013) in the V-band at the Teide, McDonald, and Haleakala observatories between 15 September 2020 and 28 December 2020, and between 14 September 2021 and 30 March 2022. We obtained 98 and 176 observing epochs, 82 and 141 of which were good, respectively, and we acquired 30 and 50 individual exposures of 1/2 s per epoch in the 40 cm telescopes. The 40 cm telescopes are equipped with a Зk×2k SBIG CCD camera with a pixel scale of 0.571 arcsec, providing a field of view of 29.2×19.5 arcmin. Weather conditions at observatories were mostly clear, and the average seeing varies from 1.0” to 3.0”. Raw data were processed using the BANZAI pipeline (McCully et al. 2018), which includes bad pixel, bias, dark, and flat-field corrections for each individual night. We performed aperture photometry in the V–band for HD 63433 and two reference stars of the field and obtained the relative differential photometry between the target and the brightest reference. We adopted an aperture of 10 pixels (~6”), which minimizes the dispersion of the differential light curve. The light curve has a dispersion of σ_Lco ~ 46 mmag, and the mean value of the uncertainties is ~l7mmag.

2.3 CARMENES

Between 19 September 2020 and 23 February 2022, we collected 157 high-resolution spectra, divided into two campaigns, with the Calar Alto high-Resolution search for M dwarfs with Exoearths with Near-infrared and optical Echelle Spectrographs (CARMENES) mounted on the 3.5 m telescope at Calar Alto Observatory, Almería (Spain), under the observing programs H20-3.5-027, 21B-3.5-015, and 22A-3.5-009. To properly model the stellar activity, we designed an observational strategy to obtain three to five spectra per stellar rotation period (~6.4 days; see Sect. 3.2 for details). Finally, we obtained about four spectra per stellar rotation on average.

The CARMENES spectrograph has two channels (Quirrenbach et al. 2014, Quirrenbach et al. 2018). The visible (VIS) channel covers the spectral range 520-960 nm, and the near-infrared (NIR) channel covers the spectral range 960–1710 nm with a spectral resolution of ℛ = 94 000 and ℛ = 80 400. One spectrum of each arm was ruled out because the drift correction was missing. Moreover, another six spectra were discarded in each arm due to their low signal-to-noise ratio (S /N < 50). The remaining observations were taken with exposure times of 150 s, resulting in S/Ns per pixel in the range of 68–250 at 745 nm. We used the VIS and NIR channel observations to derive RV measurements. The CARMENES performance, data reduction, and wavelength calibration were made using CARACAL and are described in Trifonov et al. (2018) and Kaminski et al. (2018). Relative RV values and activity indexes such as the Ha index, the Can IR triplet (IRT), and the NaiD values were obtained using the serval³ pipeline (Zechmeister et al. 2018). Furthermore, the CRX and dLW activity indicators were also computed, where the first indicator takes into account the dependence of the RV on the wavelength, and the second evaluates the variations in the widths of the lines with respect to the reference lines. In addition, this software also produces a high S/N template spectrum by co-adding all the observed spectra after correcting them for the wavelength shifts. The RV measurements were corrected for barycentric motion, secular acceleration, instrumental drift, and nightly zero-points. We searched for possible outliers in the datasets and discarded 4 points in the NIR time-series. This resulted in 150 and 146 RV data points in the VIS and NIR channels, respectively. The typical dispersion of the RV measurements is σ_{carmenes vis} ~ 19.8 ms⁻¹ for the VIS range and σ_{carmenes nir} ~ 25.0 m s⁻¹ in the NIR range. The uncertainties are in the range of 2.6–7.5 m s⁻¹ with a mean value of 3.9 m s⁻¹ in the VIS arm and in the range of 8–38 m s⁻¹ with a mean value of 14 m s⁻¹ in the NIR. The RV curve for both datasets is illustrated in Fig. 3 with its best-fit model (see Sect. 6 for details). We searched for RV measurements affected by flares by measuring the relative intensity of a set of spectral lines that are usually associated with stellar activity (Hα, Can IRT, Na I, and K I). We compared them with each other to look for significant variations, but none of the spectra lines seem to be affected by flares. In addition, because a significant part of our RV data is contemporaneous with the TESS photometry, we visually inspected whether any of the flares found in the TESS light curve coincided with the RV data. This was not the case, therefore we kept all RV data. Tables C.1 and C.2 give the time stamps of the spectra in BJD_TDB and the relative RVs measured with serval along with their 1σ error bars.

Fig. 3 CARMENES RV data for HD 63433 (blue and orange dots for the visible and infrared, respectively). Top panel: combined model (black line) with its 1σ level of confidence (grey shadow), and the dashed red line depicts the Keplerian model for two planets. Middle panel: Keplerian model alone (dashed red line) and CARMENES VIS and NIR data after subtracting the best activity model. Bottom panel: residuals for the best-fit.

3 Stellar properties

3.1 Physical parameters of HD 63433

HD 63433 is a G5V star (Gray et al. 2003) located at a distance of 22.34 ± 0.04 pc. The distance was determined from the Gaia DR3 parallax (Gaia Collaboration 2023). HD 63433 was first identified as a member of the Ursa Major moving group by Gaidos (1998) and was later confirmed by Gagné et al. (2018). By measuring the lithium abundance and rotation period of HD 63433, Mann et al. (2020) confirmed the young age of the star and the membership to the group. We adopt the recent age determination of 414 ± 23Myr for the HD 63433 system that was provided by Jones et al. (2015) for the Ursa Major moving group.

Using the stellar template, which combines all the 150 spectra, we computed the stellar atmospheric parameters of HD 63433 with the STEPARSYN code⁴ (Tabernero et al. 2022). This code implements the spectral synthesis method with the emcee⁵ Python package (Foreman-Mackey et al. 2013) to retrieve the stellar atmospheric parameters. We employed a grid of synthetic spectra computed with the Turbospectrum (Plez 2012) code, the MARCS stellar atmospheric models (Gustafsson et al. 2008), and the atomic and molecular data of the Gaia-ESO line list (Heiter et al. 2021). We employed a set of Fe I,II features that are well suited for analysing FGKM stars (Tabernero et al. 2022). We retrieved the following parameters: T_eff = 5553 ± 56 K, log g = 4.56 ± 0.08 dex, [Fe/H] = –0.07 ± 0.03 dex, and υ sin i = 7.76 ± 0.10 km s⁻¹.

We estimated the luminosity of HD 63433 by integrating the observed fluxes from the UV-optical to mid-infrared using VOSA (Bayo et al. 2008), including the Galaxy Evolution Explorer (GALEX; Bianchi et al. 2017), Tycho (Høg et al. 2000), Gaia (Gaia Collaboration 2023), the Sloan Digital Sky Survey (SDSS; Abdurro’uf et al. 2022), Two Micron All Star (2MASS; Skrutskie et al. 2006), AKARI (Murakami et al. 2007), and Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) passbands. We used BT-Settl (CIFIST) models (Baraffe et al. 2015) to reproduce the spectral energy distribution (SED) of the star and to extrapolate to bluer and longer wavelengths. We obtained a luminosity of 0.748 ± 0.010 L_⊙ for HD 63433, which is very similar to previous results by Mann et al. (2020). To estimate the radius and mass of the star, we followed the procedure presented in Schweitzer et al. (2019). Based on the estimated effective temperature and luminosity and using the Stefan-Boltzman relation, we derived a radius of 0.934 ± 0.029 R_⊙. Based on this radius and using the empirical mass-radius relations forsolar type stars from Eker et al. (2018), we determined a mass of 0.956 ± 0.022 M_⊙, assuming that the star is on the main sequence, as expected for its age. A summary with the main stellar parameters of HD 63433 can be found in Table 1.

3.2 Stellar rotation and activity

As a consequence of the youth and fast rotation of HD 63433, the star shows a remarkable stellar activity that is clearly visible in photometry as well as in RV due to the rotation of the star. HD 63433 has a stellar rotation period of 6.45 ± 0.05 days (Mann et al. 2020), as estimated from the combination of the Lomb-Scargle periodogram and the autocorrelation function of the TESS light curve in sector 20. We searched for periodic signals between 2 and 30 days using the larger temporal coverage of sectors 44–47 (more than 100 days of quasi-continuous observations), our dedicated photometric monitoring with LCO, and spectroscopic campaigns with the CARMENES instrument, using the RV and the stellar activity indicators that are generated by serval for the CARMENES VIS and NIR datasets. GLS periodograms (Zechmeister & Kürster 2009) were computed for these data and are shown in Fig. 4. The top two panels show the GLS periodograms of TESS for all the sectors and of each one of them, respectively. The third panel shows the LCO peri-odogram, and the following panels present the RV periodograms and different activity indicators for both the VIS and NIR. The periodograms of the TESS, dLW (NIR), Hα, and Ca II indexes show a period between 6.4 and 6.8 days as the most significant signal, in agreement with the published rotation period (Mann et al. 2020). On the other hand, the RV periodogram shows that the strongest signal is at ~3.2 days, that is, at half the rotation period. This signal is also very intense and comparable with the 6.4- to 6.8-day signal in the TESS periodogram. It is easily recognizable in Fig. 2. The signal of one-third of the rotation period (centred on 2.1 days) is also detected in the dLW (VIS) panel, in the RV, and the dLW (NIR) periodogram, but with a lower significance. Finally, the LCO, CRX, and Na I indicators do not show any significant signal. The photometric precision obtained with the LCO data in the V band (~ 17–46 mmag) was not sufficient to detect the photometric variations of the star, which show a peak to peak of ~20 mmag in TESS data, partially because there were no stars of similar brightness in the field of view and because the chromaticity of these activity signals is low. For this reason, and because no significant signals are actually observed in the periodogram, we decided not to use these data in the subsequent analysis. We conclude that the signal of the stellar rotation period is present (in its fundamental frequency and/or its harmonics) in most of our data. However, it is remarkable that the most significant peaks in the TESS periodograms (probably the best indicator, considering the cadence and baseline of the data) have a similar intensity and a FWHM of ~0.4 days. The second panel shows that the peak values change slightly from one sector to the next and between the fundamental signal and its first harmonic. This may be an indication of variations in the rotation period determination that may be caused by differential rotation or by variations in the activity regions on timescales of a few rotation periods. Lastly, to study the correlation between the CARMENES data products, we computed the Pearson r coefficient. We found correlations of p > 0.8 between RV VIS and RV NIR on the one hand and between the Ca II indices on the other hand. We found no evidence of a correlation for the other products.

In Sects. 4–6 we use Gaussian process regression (GP; Rasmussen & Williams 2006), which is implemented in the celerite (Foreman-Mackey et al. 2017) and George (Ambikasaran et al. 2015) packages, to derive the hyperparameters of stellar activity, including the parameter associated with the period. Using the data from TESS, CARMENES VIS and NIR, we measured a P_rot = 6.40 ± 0.02 days (see Sect. 6 and Table 3), consistent with the values analysed in the peri-odograms and with previous results. For all the considerations given above, we adopted a final rotation period of 6.4 ± 0.4 days, where we used the FWHM of the main peak of the TESS periodogram as the error (Fig. 4, top panel).

Fig. 4 GLS periodograms for the photometric light curves, RV, and spectral activity indicators from TESS, LCO, and CARMENES VIS and NIR data. The first panel shows the periodogram of the combination of all TESS sectors. A GLS periodogram for each sector (sectors 20, 44, 45, 46, and 47 depicted in grey, green, cyan, brown, and black, respectively) is shown in second panel. The third panel shows the periodogram of the LCO dataset as a solid black line. In the remaining panels, the periodograms calculated from CARMENES VIS data are displayed in blue, and those calculated for the NIR arm are plotted in red. In all panels, the solid purple and vertical orange lines indicate the orbital periods of planets b and c, respectively. The shaded yellow bands indicate the rotation period (6.4 days), half (3.2 days), and one-third of the rotation period (2.1 days). The dashed horizontal black lines correspond to FAP levels of 10%, 1%, and 0.1% (from bottom to top). The first two panels only include the 0.1% FAP multiplied by a factor 100 for clarity.

4 Transit analysis

Mann et al. (2020) identified four transits of HD 63433 b and two transits of HD 63433 c in TESS sector 20. We employed the new TESS sectors (44, 45, 46, and 47) to update the planetary transit parameters and to search for additional new planets, taking advantage of more than 100 days of continuous observations. This allowed us to identify planets with longer orbital periods. However, as shown in Figs. 2 and 3, young stars are dominated by intense periodic stellar activity with timescales of several days. This represents the main difficulties in finding young planets. Therefore, it is necessary to find the appropriate model to capture this type of activity. We chose GPs as an effective model that is flexible enough to model the variations of the amplitudes and the quasi-periodic behaviour that this stellar activity shows. The activity may show changes between two periods and even within the same cycle (Angus et al. 2018). In particular, to model the photometric stellar activity, following a similar kernel as in Suárez Mascareño et al. (2021), we used the double simple harmonic oscillator (dSHO), built as the sum of two SHO kernels, defined as $$\begin{array}{r}\hfill k_{\text{dSHO}}\left(\tau \right)=k_{\text{SHO}}\left(\tau ;{\eta }_{{\sigma }_1},{\eta }_{L_1},{\eta }_{\text{p}}\right)+k_{\text{SHO}}\left(\tau ;{\eta }_{{\sigma }_2},{\eta }_{L_2},{\eta }_{\text{p}}/2\right)\\ \hfill ={\eta }_{{\sigma }_1}^2{\text{e}}^{-\frac{\tau }{{\eta }_L{}_{{}_1}}}\left[\cos \left({\eta }_1\frac{2\pi \tau }{\eta \text{p}}\right)+{\eta }_1\frac{\eta \text{p}}{2\pi {\eta }_L{}_{{}_1}}\sin \left({\eta }_1\frac{2\pi \tau }{\eta \text{p}}\right)\right]\\ \hfill +{\eta }_{{\sigma }_2}^2{\text{e}}^{-\frac{\tau }{{\eta }_L{}_{{}_2}}}\left[\cos \left({\eta }_2\frac{4\pi \tau }{\eta \text{p}}\right)+{\eta }_2\frac{\eta \text{p}}{4\pi {\eta }_L{}_{{}_2}}\sin \left({\eta }_2\frac{4\pi \tau }{\eta \text{p}}\right)\right],\end{array}$$(1)

where τ ≡ |t_i – t_j| is the time difference between two data points, η ≡ |1 – (2πη_L|η_P)⁻²\)^1/2 and ${\eta }_{{\sigma }_i}$, ${\eta }_{L_i}$, and η_P are the hyperparameters that represent the amplitude of the covariance, the decay timescale, and the period of the fundamental signal, respectively. The induced stellar activity can act on different timescales, from hours or days (oscillations and granulation), and from days to months (rotation) or years (magnetic cycles). The hyperparameters of the GPs allowed us to adapt the model to the scales that we wished to fit. In our case, these were the variations due to stellar rotation. Therefore, we interpret the period of the process physically as the rotation period of the star, the amplitude of the covariance as the amplitude variations produced by the rotation of a particular configuration of the active regions, and the decay timescale as the evolution of these active regions of the star. The definition of this kernel is valid as long as η_P < 2πη_L, which it is a reasonable assumption in young stars because it is observed that a clearly periodic behaviour dominates in the activity. This kernel has been widely used in the literature (David et al. 2018a; Mann et al. 2020; Newton et al. 2021; Tofflemire et al. 2021) to model the light curves of young stars, for which the rotation period and its first harmonic is easily identifiable. This is consistent with our periodogram (TESS panel in Fig. 4), and when there is a sufficient cadence of data because it allows great flexibility. An instrumental offset (γ_TESS) was also considered for the TESS dataset, and a jitter (σ_{jit, TESS}) term was added in quadrature to the error bars. To explore the parameter space in our analysis, we employed an affine-invariant ensemble sampler (Goodman & Weare 2010) for the Markov chain Monte Carlo (MCMC) implemented in the emcee code. The parameter space here was also explored with another sampler algorithm, dynesty⁶ (Speagle 2020), which is based on nested sampling (Skilling 2004). The next sections show that the results are consistent.

4.1 Transit search

We searched for transits in the TESS light curves using the box least-squares periodogram (BLS; Kóvacs et al. 2002; Hartman & Bakos 2016) in each individual sector, but also in the combination of all of them. We modelled the stellar activity and subtracted it as mentioned above, and then applied the BLS to search for transits. To do this, we forced (setting normal priors) the amplitude covariance hyperparameters (${\eta }_{{\sigma }_1}$, ${\eta }_{{\sigma }_2}$) to values close to the standard deviation of the light curve (~4.3ppt) and fixed the σ_jit,TESS parameter to ~0.8 ppt, which corresponds to the depth of the largest transit of the two known planets (HD 63433 c). In this way, we modelled the stellar activity smoothly and prevented the GPs from removing or modelling possible transits.

After the stellar activity was subtracted, BLS found HD 63443 b as the most significant signal with 18 transits (vertical orange bands in Fig. 2). Then, we masked out this signal and iteratively applied the BLS algorithm again to search for additional features. It found 8 transits of HD 63433 c as the second strongest signal (vertical purple bands in Fig. 2). After this, we found no additional transit. Furthermore, we inspected the light curves by eye for isolated transits and found no variations compatible with transits.

4.2 Transit parameters

To improve the ephemerides of the two known planets through a larger baseline, we proceeded to simultaneously fit our photometric model as a combination of stellar activity and two planetary transit signals. To model the planetary transits, we used PyTгansit⁷ (Parviainen 2015), and the stellar activity was modelled in the same way as in Sect. 4.1, now leaving σ_jit,TESS as a free parameter. Our transit-only fit contains the following planetary parameters: the time-of-transit centre (T_c), the planetary radius (R_p), the orbital period of the planet (P), and the impact parameter (b). Models with non-circular orbits were also included. The eccentricity (e) and argument of perias-tron (ω) in them were parametrized as in Anderson et al. (2011, $\sqrt{e}\sin \omega$, $\sqrt{e}\cos \omega$. The input for stellar parameters was the stellar mass (M_⋆), the stellar radius (R_⋆), and the linear and quadratic limb-darkening coefficients (q₁,q₂). In the last case, we used the parametrization of Kipping (2013), where the initial values were previously calculated with the Python Limb Darkening Toolkit⁸ (PyLDTk; Parviainen & Aigrain 2015). We chose to sample the hyperparameters on a natural logarithmic scale for faster convergence. In addition, we partially constrained some of them assuming the following hypotheses: As mentioned before, the standard deviations of the process were favoured by setting normal priors to the scales of the standard deviation of the activity, the decay timescale was limited to values clearly greater than the rotation period and lower than the baseline of the data sample, and finally, the period of the process was allowed to sample between the main signals of activity (6.4 days, 3.2 days, and 2.1 days). The priors and posterior results are presented in Table A.1, assuming circular and eccentric orbits for the planets. Our posterior parameters are consistent with those obtained by Mann et al. (2020).

Fig. 5 TTVs for planet b (top panel) and planet c (bottom panel) represented as black dots with their 1σ uncertainty.

4.3 Transit-timing variations

Because we have a long temporal coverage (~750 days) and 18 and 8 transits for planet b and c, respectively, the transit-timing variations (TTVs) can be a tool for measuring the planetary masses. To do this, we proceeded in the same way as in the previous section, but used T_c for each individual transit. Figure 5 shows that the TTVs of planet b and c present variations consistent with zero within the error bars. Therefore, we conclude that no TTVs are detected in the system to measure planetary masses.

5 RV analysis

With the new ephemeris of the planets derived in the previous section, we checked whether some RV epochs fall during the transits because these RVs may be altered due to the RM effect. Three RV measurements fall during the transit of HD 63433 b (BJD = 2459513.5797, 2459634.4229) and HD 63433 c (BJD = 2459131.6311). Because the expected variation in amplitude due to RM effect of HD 63433 b and c can be a few ms⁻¹ (Dai et al. 2020; Zhang et al. 2022), we decided to remove these RVs from our data. Overall, our final dataset contains 147 RV points of CARMENES VIS and 143 of CARMENES NIR.

5.1 Periodogram analysis

In Fig. 6, we show a more detailed study of the RV periodograms calculated in Sect. 3.2. The CARMENES VIS and NIR datasets are highlighted in blue and red, respectively. The orbital periods of the transiting planets b and c are marked with vertical purple and orange lines, respectively. In the top panel in Fig. 6, we plot the RV data, where two peaks with a high significance (FAP<0.1%) are seen at ~ 3.2 and ~ 2.1 days. They correspond to the harmonics of the stellar rotation period (vertical yellow band) estimated in Sect. 3.2. As previously commented, the signal of the planets that orbit or transit young stars is usually hidden by the strong stellar activity, which is the dominating signal in periodograms. As a first approximation, we modelled the most significant signals with a sinusoidal function before subtraction. The second panel shows the residuals after the signal related to half of the stellar rotation period was subtracted. This is the most significant signal in both periodograms. After removal of this activity signal, the 2.1d signal becomes the most significant in both datasets, while the signals of the planets remain hidden. The data from both arms seem to be dominated by stellar activity so far, with a strong correlation between the two datasets (Sect. 3.2). In the third panel, we repeat the process and subtract the 2.1d signal. In the periodogram of the CARMENES VIS residuals, only one signal remains with an FAP<0.1%: the signal at 6.4 days. This signal is at the rotation period of the star. However, in the CARMENES NIR dataset, all the peaks have a significance lower than 0.1% FAP. In the fourth panel, we subtracted the three most significant signals associated with stellar activity only in the CARMENES VIS data. The signal with an FAP~1% close to the orbital period of planet c does not appear in the periodogram of CARMENES NIR because the larger uncertainties in the NIR measurements make it more difficult to detect signals with smaller amplitudes. In the fifth panel, we subtracted the stellar activity for each dataset by modelling in another way, using GP (see Sect. 5.2). As with the previous method, the signal of the outer planet (FAP~10–1%) is visible, but not that of the inner one. The orbital period of HD 63433 b is ~7.1d and the rotational period is ~6.4 days, that is, within a 10% margin. This quasi-coupling of the signals, together with the great difference in amplitude between them, makes it especially difficult to detect the planetary signal for the inner planet. Finally, in the sixth panel, we show the window functions of the VIS and NIR datasets. In summary, from the analysis of the periodograms of the VIS and NIR RV, we conclude that only the activity signal in both datasets and only the outer planet signal in CARMENES VIS is detected with sufficient significance.

Fig. 6 GLS periodograms for the CARMENES VIS (blue) and NIR (red) RV datasets. In all panels, the solid vertical purple and orange lines indicate the orbital periods of planets b and c, respectively. The shaded vertical yellow bands indicate the rotation period (6.4 days) and the half (3.2 days) and one-third (2.1 days) harmonics. The dashed, dotted-dashed, and dotted horizontal black lines correspond to FAP levels of 10%, 1%, and 0.1% from bottom to top, respectively.

5.2 RV modeling

After studying the signals present in the periodograms, we tested a battery of RV-only fits to determine the best activity model and the best Keplerian model of the planets. We considered two possible datasets: only the VIS dataset, and the combination of VIS and NIR measurements. To create a model of stellar activity, we used GP as in Sect. 4.2 and the quasi-periodic (QP) kernel defined in Aigrain et al. (2012) as $$k_{\text{QP}}\left(\tau \right)={\eta }_{\sigma }^2\exp \left[-\frac{{\tau }^2}{2{\eta }_L^2}-\frac{{\sin }^2\left(\frac{\pi \tau }{\eta \text{p}}\right)}{2{\eta }_{\omega }^2}\right],$$(2)

where τ, η_σ, η_L, and η_P are the same as in Eq. (1), and η_ω controls the weight given to the periodic part of the kernel. The QP kernel has one free parameter less than the dSHO, and there is no coupling between the periodic part and non-periodic part, as in the case of dSHO, which gives it less flexibility when the data sampling is smaller, as is the case with RV data. Once again, we imposed normal priors on the covariance amplitudes associated with the standard deviation of the data because the goal was to model stellar activity with a smooth function on period scales. As the rotation period signal and its harmonics are present in RV data, we also tested a combination of sinusoidal functions to fit the stellar activity. To do this, we considered several cases. The first case was an activity-only model, for which we assumed that planetary signals are not detected in the CARMENES RV data. We modelled the stellar activity with GPs or sinusoidal functions. The second case considered the existence of one transiting planetary signal (the inner or the outer planet) modelled as a Keplerian circular orbit with the initial planetary parameters obtained from a transit-only fit (Sect. 4.2), including the stellar activity as in the first case. The last case included a Keplerian model of both planets and the different models for stellar activity. The Keplerian model parameters were the T_c, P, and the RV amplitude of the planet (K) and were modelled with the RadVel⁹ (Fulton et al. 2018) package. In the photometric analysis, all models of the RV fit included an instrumental offset (γ_CARMENES) and a jitter term added in quadrature to the error-bar measurements (σ_jit,CARMENES).

To evaluate which model was better for each dataset, we considered the Bayesian log-evidence (ln 𝒵) criterion defined in Trotta (2008) and calculated from Díaz et al. (2016), where the model with a higher log-evidence is favoured and was used as reference (Δ ln 𝒵≡ ln 𝒵_model – ln 𝒵_{model ref}). When the difference between two models is |Δ ln 𝒵| < 1, the two models are considered indistinguishable. When 1 < |Δ ln 𝒵| < 2.5, the evidence in favour of one of the model is weak, and when the difference is 2.5 < |Δ ln 𝒵| < 5, the evidence is moderate. This criterion allows us to compare models with different numbers of parameters as long as the input observational dataset is the same. The results of the different models considered in this work are compiled in Table 2. For the CARMENES VIS dataset, the model with the best log-evidence is always the model that includes the a single planetary model, planet c (marked in bold in the table). This yields weak evidence compared to other planetary models. In contrast, including a model with planet b alone does not give a significant improvement. Moreover, better ln 𝒵 are clearly obtained when the stellar activity is modelled by GPs instead of sinusoidals. The amplitudes obtained for the planets using GPs or sinusoidal functions are comparable within 1σ, although those obtained for planet c using GPs have a slightly lower value.

In the second configuration, using both CARMENES VIS and NIR data, we only ran the activity model with GP, which has the best ln 𝒵 by far. Although the error bars in CARMENES NIR data are significantly higher, the stellar activity appears to be well detected in the periodograms. Therefore, a combined model of both datasets may help to better constrain stellar activity and, consequently, to better define planetary signals. Since there is a high correlation between VIS and NIR datasets, we have shared the hyper-parameters of η_L, η_ω, and η_P between both datasets in our model. At the bottom of Table 2, we see that the best model includes planet c, together with the activity, in the same way as with CARMENES VIS data alone. In this case, the amplitudes are similar to the previous case, but the errors are slightly lower. After this analysis, we decided to use the results that include CARMENES VIS and NIR data with two planetary models for two reasons. Although we do not see any signal that we can associate with the planets in the periodogram of CARMENES NIR, the stellar activity seems to be well detected. A better characterization of stellar activity (the hyperparameters of the model that combines the two datasets are slightly better constrained than when CARMENES VIS is used alone) appear to constrain the planetary parameters better. Finally, we also decided to adopt the two-planet models as the final model because we have a priori information of the existence of the two transiting planets, and with this model, we can place a constraint on the mass of planet b. The planetary parameters for this best RV-only fit are catalogued in Table B.1 for circular and eccentric orbits.

Fig. 7 Recovered RV amplitudes vs. orbital periods for planet b. The red dots and blue triangles show the amplitude we recovered for planets b and c, respectively. The dotted horizontal and vertical black lines indicate the injected amplitude and the orbital period for planet b, respectively.

5.3 Injection-recovery tests of planet b

As mentioned in Sect. 5.1, the proximity of the signals of the stellar rotation and the orbital period of the inner planet hamper a correct measurement of the RV amplitude of planet b. In order to test the dependence of the amplitude of planet b on the separation between the orbital and rotation period, we decided to carry out an injection-recovery test. First, we subtracted the signal of the inner planet obtained from the best RV model that combines GP and two Keplerians (K^b ~ 1.77 ms⁻¹). Then, we injected the signal of planet b with the same amplitude but slightly different orbital periods, and finally, we determined the amplitude of the inner planet using the same model as in Sect. 5.2. In Fig. 7, we represent the recovered amplitude for different orbital periods. For values close to the rotation period and its first harmonic, the recovered amplitude and the error bars can be twice or three times as high as the value of the injected amplitude. This confirms the dependence of the orbital and rotation periods. In addition, it is noteworthy that even for the orbital periods farthest from the rotation period, the uncertainty of the amplitude is never lower than ~0.8 m s⁻¹, indicating that with this dataset and model, amplitudes smaller than 2.4 m s⁻¹ could not be measured with a 3σ significance. The figure also shows that the amplitude measure of the outer planet depends on the orbital period of the inner planet. The amplitude obtained for planet c is displayed as blue triangles (and intentionally shifted to the left of each measure of planet b) to show that it is stable and independent of the measure of planet b.

6 Joint-fit analysis

After analysing the best photometry and RV models, we carried out a global modelling of the data. To do this, we used a transit model that included two planets (whose phase-folded transits are depicted in Fig. 8), a Keplerian model with two planets (whose phase-folded RV are displayed in Fig. 9), a stellar activity model that used a GP with a dSHO kernel to model photometry (Fig. 2), and a QP kernel to model the RV datasets (Fig. 3). Although we would have preferred to use the same GP for the joint fit, the large number of data in TESS (>80000) forced us to use a kernel implemented in celerite, which has a significantly lower computational cost. In the planetary models, the parameters of T_c and P share normal priors, while uniform priors were set for the rest of planetary parameters. For the activity models, we chose normal priors for the parameters of the covariance amplitudes, centred on the standard deviation of each dataset (consistent with the posteriors obtained in the independent analysis). The period of the process was shared by the two activity models, and the hyperparameters of η_L and η_ω were shared between RV datasets. All the prior and posterior results of the joint fit and their derived parameters are shown in Table 3.

7 Discussion

7.1 Densities and compositions

From the posterior parameters obtained in the joint fit, we derive an upper limit of $M_{\text{p}}^{\text{b}}<21.8M_{\oplus }$ on the mass of HD 63433 b at the level of 3σ and present a 4σ detection for HD 63433 c, yielding a mass of $M_{\text{p}}^{\text{c}}=15.5\pm 3.9M_{\oplus }$. These values, together with the radii of planets b and c updated in this work, give bulk densities of ρ^b < 13.0 g cm⁻³ and ρ^c =4.6 ± 1.3 g cm⁻³, respectively. The equilibrium temperature (T_eq) for both planets lies in the range of 769–967 K for HD 63433 b and 540–679 K for HD 63433 c, assuming planetary albedos (A) between 0 and 0.6. Our results imply that both planets are warm sub-Neptunes. The density of HD 63433 c is slightly lower than that of Earth.

Figure 10 shows a mass-radius diagram that presents all known small exoplanets (taken from the Extrasolar Planets Encyclopaedia¹⁰) with radius precisions better than 8% (measured through the transit method) and mass precisions better than 20% (measured through RV). The uncertainties of HD 63433 b and c are represented as a coloured-shaded region, according to the 1, 2, and 3σ levels of confidence. To place the planets in context, we also plot the young planets (≤900 Myr) as coloured dots according to their age. Only eight young planets have a measured radius and mass (AUMicc, TOI-1807b, TOI-560b and c, K2-25b, K2-100b, TOI-1201b, and now we add HD 63433 c). They represent less than 10% of the total population of exo-planet shown in the diagram. The left panel of Fig. 10 shows that the planets are not randomly distributed. Most of them are concentrated between the density lines of 1–10 g cm⁻³ with two main groups: those with masses lower than 8–10 M_⊕ with densities of 3–10 g cm⁻³, and those with higher masses, with typically a larger radius and lower densities. Young planets also follow this distribution so far, although only one planet, TOI-1807b, currently has the characteristic radius and mass of super-Earths. HD 63433 c also matches the planets in the field. Its density is only slightly lower than that of Earth, but is consistent within the error bars.

Radius and mass measurements can help to determine the internal composition of the planets. We plot models of different compositions by mass without an atmosphere (in the middle panel) and with atmospheres (on the right) from Zeng et al. (2019) in Fig. 10. In the first case, where the internal composition of the planets does not have a significant envelope, our measurements allow us to rule out that HD 63433 b contains more than 50% of iron by mass, and it may be compatible with planets with a higher silicate or water content. However, the composition of HD 63433 c does not contain a significant mass fraction of iron and is most probably formed by silicates and/or water. On the other hand, models with envelopes can only explain their positions in the mass-radius diagram if their atmospheres contain 2% of H₂ by mass at most for any composition mixture.

Fig. 8 Top panels: phase-folded TESS light curves of the transits of planets b and c (blue points), binned data (white dots), and our best transit-fit model (black line). Bottom panels: residuals for the best fit.

Fig. 9 Top panels: RV-folded CARMENES data for planets b and c (blue and orange dots for VIS and NIR, respectively), binned data (white dots), and our best Keplerian-fit model (black line). Bottom panels: residuals for the best fit.

7.2 Formation history

Formation models of planets with orbital periods shorter than 100 days suggest that the most massive planets (with a core mass greater than 5–10 M_⊕) should form in outer regions (beyond the ice line) where the mean size of ice pebbles is greater than the size of the silicates (Morbidelli et al. 2015; Venturini et al. 2020). Higher masses in this region also favour gas accretion, which notably increases the size of the planets. Then, interactions with the disc gas are expected to trigger the migration of several of these planets into the inner regions, causing them to mix with rocky planets that have already formed there. This mechanism is expected to occur on timescales shorter than 10 Myr, which is the typical lifetime of protoplanetary mass discs (Williams & Cieza 2011). The discovery of a Neptune-size planet in the Upper Scorpius region (5–10 Myr) confirms this (David et al. 2016; Mann et al. 2016b). At this stage, evolutionary mechanisms such as photoevaporation (Owen & Wu 2017), which takes place on timescales of several tens of million years, or core-powered mass-loss (Ginzburg et al. 2018; Gupta & Schlichting 2020), which acts on timescales of ~ 1 Gyr, may play an important role in removing part of these primordial atmospheres of short-period planets.

Young, multi-planet systems offer a unique opportunity to study the formation and the effects of evolution on planets. HD 63433 has an estimated intermediate age of ~400 Myr and therefore constitutes a good example on which these mechanisms can be tested. The mass measurement of HD 63433 c and the upper limit set on the mass of HD 63433 b presented in this paper allow us to determine the composition of these planets and discuss them in the context of their evolution. In this sense, planet HD 63433 c is mostly composed of rocks and water, with little or no gas envelope, while HD 63433 b may have a similar composition or a larger part of gas. Recently, Zhang et al.(2022) have shown that the atmosphere of both planets does not show signs of an intense evaporation. These results suggest that if HD 63433 c had a significant gas envelope in the past, as expected for such a massive planet according to formation models, it has already lost most of its gas content at the age of system. This favours rapid evolutionary mechanisms, such as photo-evaporation, which can explain the mass loss of the planet atmospheres on timescales shorter than a few hundred million years.

All the young planets shown in Fig. 10, except for TOI-1807 (catalogued as an ultra-short period planet), are on the upper side of the radius valley (R > 1.9 R_⊕), which corresponds to Neptune/sub-Neptune planets whose composition could be similar to HD 63433 c. This is probably an observational bias, as mentioned earlier, because of the high stellar activity of young stars that makes the detection and measurements of the physical properties of small planets around these stars difficult. Therefore, young planets are typically found with larger radii and higher masses because they are easier to detect. A recent work (Luque & Pallé 2022) has shown that the population of small planets orbiting M dwarfs is composed of three types of planets with different compositions: rocky planets, water worlds (50% water, 50% silicates), and small gas planets. Although HD 63433 c (and possibly HD 63433 b also) and other young small planets, such as TOI-1201 b, TOI-560b and c, seem to follow the sequence of water planets, while more massive young planets, such as K2-100 b or AU Mic b and c, seem to be small gas planets, more similar to Neptune and Uranus, the number is still low. It is necessary to obtain more dynamical mass determinations of young planets to compare the atmospheric compositions of young and old populations.

Fig. 10 Radius-mass diagram of HD 63433 b and c, together with all exoplanets (grey dots) with a precision better than 8% in radius (by transits) and 20% in mass (by RV). Young exoplanets with measured masses are plotted as coloured dots according to their ages. The uncertainties on the HD 63433 b and c planets are shown as coloured shaded regions with 1, 2, and 3σ levels of confidence. The dashed grey lines are constant density lines. Coloured lines indicate different composition models without gas (middle panel) and with a gas envelope (right panel) from Zeng et al. (2019). Earth and Neptune are also depicted as reference. B22 and M22 refer to Barragán et al. (2022) and El Mufti et al. (2023), respectively.

Fig. 11 Comparison between photometric and RV stellar activity. Top panel: TESS light curve of HD 63433 contemporaneously with the RV of CARMENES data. The blue dots are the PDCSAP flux, and the black line indicates the best activity model. Bottom panel: CARMENES RV VIS data shown with blue dots. The black line shows the GP activity model with its 1σ confidence level (grey shadow). The solid red line depicts the RV model predicted from the TESS light curve using FF′ method.

7.3 Dependence on flux and RV

Aigrain et al. (2012) presented a method with analytical formulas to predict the stellar activity in RV from the observed photometric fluxes. The high activity of young stars makes them an ideal laboratory for testing this FF′ method. Because part of our observations with CARMENES were intentionally obtained simultaneously with TESS, we can observationally verify the behaviour of stellar activity in photometry and RV. The upper panel of Fig. 11 shows the TESS data contemporaneously with the CARMENES RVs and the activity model obtained in the transit-only fit. The lower panel of Fig. 11 displays the CARMENES VIS data together with the activity model (solid black line) obtained in the RV-only fit for the same time interval as in the upper panel. The stellar activity in TESS presents a double modulation, in the same way as the RV. This agrees with the analysis of the RV periodograms, whose the maximum peak is at 3.2d, in agreement with the first harmonic of the rotation period. In addition, in the RV observations where the cadence is better (between 2459550 and 2459570 BJD), the model presents up to three peaks per rotation period. This also agrees with the analysis of the periodograms of the RV, and some indicators also show a significant peak at 2.1 days. Using the analytical formulas in Aigrain et al. (2012) based on a simple spot model, we calculated the ΔRV_rot from the stellar activity modeled in TESS. In the same panel, we represent the result of the FF′ method (normalized by a scale factor) with a solid red line together with the observed RV data. Both RV models, the one calculated from GPs and the other calculated with the FF′ method from photometry, behave similarly throughout the time series. Although the FF′ model never reproduces three peaks per rotation period, both models reproduce the amplitude variations between one period and another and between the two main peaks within each period in the same way. That is, the FF′ method predicts similar stellar activity as the GP model based on photometry and with the ΔRV_rot term alone, suggesting that most of the stellar activity affecting the photometry also affects the RV in a similar way. Therefore, obtaining contemporaneous photometry with high cadence as in TESS together with RV data maybe be a valuable resource for characterizing the source of stellar activity and correct it with the appropriate models to determine the masses of young or low-mass planets.

8 Conclusions

We have analysed precise photometric light curves of five sectors of the TESS mission and obtained ~150 precise RV measurements with CARMENES VIS and NIR of the HD 63433 star. We performed a joint-fit analysis of these data and present an update of the transit parameters and a characterization of the dynamic masses of the HD 63433 planetary system. The inner planet, HD 63433 b, has an orbital period of P^b = 7.108 days, a radius of $R_{\text{p}}^{\text{b}}=2.140\pm 0.087R_{\oplus }$, and an upper limit at 3σ level for the amplitude of $K_{\text{p}}^{\text{b}}<7.41\text{\hspace{0.17em}m\hspace{0.17em}}{\text{s}}^{-1}$. From these values, we derive a mass of $M_{\text{p}}^{\text{b}}<21.76M_{\oplus }$ and a planetary density of ρ^b < 13.00 ρ_⊕ On the other hand, for planet HD 63433 c, we obtain a period of P^c = 20.544 days, with a radius of $R_{\text{p}}^{\text{c}}=2.692\pm 0.108R_{\oplus }$ and an amplitude of $K_{\text{p}}^{\text{c}}=3.74\pm 0.93\text{m\hspace{0.17em}}{\text{s}}^{-1}$, translating into a mass of $M_{\text{p}}^{\text{c}}=15.54\pm 3.86M_{\oplus }$ and a density ρ^c = 4.63 ± 1.30 ρ_⊕.

According to theoretical models, planet HD 63433 c is composed mostly of rocks and water and has little or no gas envelope, while HD 63433 b may have a similar composition, but its gas content remains unconstrained due to the lack of a mass measurement. These results, together with the lack of signs of intense evaporation in both planets (Zhang et al. 2022), suggest that if HD 63433 c had a significant gas envelope in the past, it has already lost most of it at about 400 Myr. This favours rapid evolutionary mechanism of mass loss, such as photoevaporation.

We acknowledge that HD 63433 has been followed up also with the HARPS-N spectrograph (Damasso 2023). Our team and the GAPS team have coordinated the submission of two studies, which were intentionally carried out in an independent way.

Acknowledgements

M.M.D., N.L., V.J.S.B., M.R.Z.O., H.T., A.S.M., E.P., and D.M. acknowledge support from the Agencia Estatal de Investigación del Ministerio de Ciencia e Innovación (AEI-MCINN) under grant PID2019-109522GB-C5[1, 3, 4]. CARMENES is an instrument for the Centro Astronómico Hispano-Alemán de Calar Alto (CAHA, Almería, Spain). CARMENES is funded by the German Max-Planck-Gesellschaft (MPG), the Spanish Consejo Superior de Investigaciones Científicas (CSIC), the European Union through FEDER/ERF FICTS-2011-02 funds, and the members of the CARMENES Consortium (Max-Planck-Institut für Astronomie, Instituto de Astrofísica de Andalucía, Landessternwarte Königstuhl, Institut de Ciències de l’Espai, Institut für Astrophysik Göttingen, Universidad Complutense de Madrid, Thüringer Landessternwarte Tautenburg, Instituto de Astrofísica de Canarias, Hamburger Sternwarte, Centro de Astrobiología and Centro Astronómico Hispano-Alemán), with additional contributions by the Spanish Ministry of Economy, the German Science Foundation through the Major Research Instrumentation Programme and DFG Research Unit FOR2544 “Blue Planets around Red Stars”, the Klaus Tschira Stiftung, the states of Baden-Württemberg and Niedersachsen, and by the Junta de Andalucía. Based on observations collected at the Centro Astronómico Hispano-Alemán (CAHA) at Calar Alto, operated jointly by Junta de Andalucía and Consejo Superior de Investigaciones Científicas (IAA-CSIC) under programmes H20-3.5-020 and F21-3.5-010. Part of this work was supported by the German Deutsche Forschungsgemeinschaft, DFG project number Ts 17/2-1. We acknowledge the use of public TESS data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products. This paper includes data collected with the TESS mission, obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the TESS mission is provided by the NASA Explorer Program. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. This research has made use of the Simbad and Vizier databases, operated at the centre de Données Astronomiques de Strasbourg (CDS), and of NASA’s Astrophysics Data System Bibliographic Services (ADS).

Appendix A Transit-only fit

Appendix B RV-only fit

Appendix C RV data

References

All Tables

All Figures

	Fig. 1 TESS TPF plots for HD 63433 (white crosses). The red squares indicate the best optimal photometric aperture used to obtain the SAP flux. G-band magnitudes from Gaia DR3 are shown with different sizes of red circles for all nearby stars with respect to HD 63433 up to 8 magnitudes fainter.
In the text

	Fig. 2 Light curves of HD 63433 for the five TESS sectors. The blue dots correspond to the PDCSAP flux data. The black line and the grey shaded region indicate the best model and its 1σ uncertainty, respectively. The vertical lines show the times of the planetary transits for HD 63433 b (orange) and HD63433 c (purple).
In the text

Fig. 3 CARMENES RV data for HD 63433 (blue and orange dots for the visible and infrared, respectively). Top panel: combined model (black line) with its 1σ level of confidence (grey shadow), and the dashed red line depicts the Keplerian model for two planets. Middle panel: Keplerian model alone (dashed red line) and CARMENES VIS and NIR data after subtracting the best activity model. Bottom panel: residuals for the best-fit.

In the text

Fig. 4 GLS periodograms for the photometric light curves, RV, and spectral activity indicators from TESS, LCO, and CARMENES VIS and NIR data. The first panel shows the periodogram of the combination of all TESS sectors. A GLS periodogram for each sector (sectors 20, 44, 45, 46, and 47 depicted in grey, green, cyan, brown, and black, respectively) is shown in second panel. The third panel shows the periodogram of the LCO dataset as a solid black line. In the remaining panels, the periodograms calculated from CARMENES VIS data are displayed in blue, and those calculated for the NIR arm are plotted in red. In all panels, the solid purple and vertical orange lines indicate the orbital periods of planets b and c, respectively. The shaded yellow bands indicate the rotation period (6.4 days), half (3.2 days), and one-third of the rotation period (2.1 days). The dashed horizontal black lines correspond to FAP levels of 10%, 1%, and 0.1% (from bottom to top). The first two panels only include the 0.1% FAP multiplied by a factor 100 for clarity.

In the text

	Fig. 5 TTVs for planet b (top panel) and planet c (bottom panel) represented as black dots with their 1σ uncertainty.
In the text

Fig. 6 GLS periodograms for the CARMENES VIS (blue) and NIR (red) RV datasets. In all panels, the solid vertical purple and orange lines indicate the orbital periods of planets b and c, respectively. The shaded vertical yellow bands indicate the rotation period (6.4 days) and the half (3.2 days) and one-third (2.1 days) harmonics. The dashed, dotted-dashed, and dotted horizontal black lines correspond to FAP levels of 10%, 1%, and 0.1% from bottom to top, respectively.

In the text

	Fig. 7 Recovered RV amplitudes vs. orbital periods for planet b. The red dots and blue triangles show the amplitude we recovered for planets b and c, respectively. The dotted horizontal and vertical black lines indicate the injected amplitude and the orbital period for planet b, respectively.
In the text

	Fig. 8 Top panels: phase-folded TESS light curves of the transits of planets b and c (blue points), binned data (white dots), and our best transit-fit model (black line). Bottom panels: residuals for the best fit.
In the text

	Fig. 9 Top panels: RV-folded CARMENES data for planets b and c (blue and orange dots for VIS and NIR, respectively), binned data (white dots), and our best Keplerian-fit model (black line). Bottom panels: residuals for the best fit.
In the text

Fig. 10 Radius-mass diagram of HD 63433 b and c, together with all exoplanets (grey dots) with a precision better than 8% in radius (by transits) and 20% in mass (by RV). Young exoplanets with measured masses are plotted as coloured dots according to their ages. The uncertainties on the HD 63433 b and c planets are shown as coloured shaded regions with 1, 2, and 3σ levels of confidence. The dashed grey lines are constant density lines. Coloured lines indicate different composition models without gas (middle panel) and with a gas envelope (right panel) from Zeng et al. (2019). Earth and Neptune are also depicted as reference. B22 and M22 refer to Barragán et al. (2022) and El Mufti et al. (2023), respectively.

In the text

Fig. 11 Comparison between photometric and RV stellar activity. Top panel: TESS light curve of HD 63433 contemporaneously with the RV of CARMENES data. The blue dots are the PDCSAP flux, and the black line indicates the best activity model. Bottom panel: CARMENES RV VIS data shown with blue dots. The black line shows the GP activity model with its 1σ confidence level (grey shadow). The solid red line depicts the RV model predicted from the TESS light curve using FF′ method.

In the text