© ESO 2021

1 Introduction

In the last few years, initial observations and studies of young close-in exoplanets have been underway (e.g. Rizzuto et al. 2017; Berger et al. 2018a; Carleo et al. 2020). A crucial role in this field is played by space-based surveys for transiting planets, such as those performed with the NASA Kepler/K2 (Borucki et al. 2003; Howell et al. 2014) and the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015) missions. The latter is providing a significant number of close-in planet candidates around young stars (e.g. Nardiello et al. 2020; Battley et al. 2020; Nardiello 2020) thanks to systematic observations of about 85% of the sky. Early discoveries of hot Jupiters around very young stars (Donati et al. 2016; Yu et al. 2017; Rizzuto et al. 2020) have suggested a greater frequency of close-in massive planets with respect to older stars. The discrepancy between the two age regimes was interpreted with different scenarios, considering an intense atmospheric escape due to the strong XUV irradiation from the young host or planet engulfment in early evolutionary stages (Donati et al. 2017).

The apparent high frequency of young hot Jupiters has been questioned based on the recent results by Carleo et al. (2018), Donati et al. (2020) and Damasso et al. (2020), which have placed doubts on the detection of such type of companions. Moreover, the first results of TESS observations suggest that the typical radius of close-in planets with ages of <100 Myr is lower than one Jupiter radius (e.g. Benatti et al. 2019; Plavchan et al. 2020; Carleo et al. 2021). The ongoing contraction and the effect of the intense high-energy irradiation from their host star suggest the possibility that many of them could actually be highly inflated lower mass planets, as was first shown by Kraus et al. (2015) and Mann et al. (2016), and confirmed by theoretical predictions as in Linder et al. (2019). This working hypothesis implies:

• A dependency of the planetary mass-radius relation on the age, which is well constrained for the older population (e.g. Chen & Kipping 2017; Bashi et al. 2017). However, this aspect is still unexplored for young planets ranging from a few Myr up to hundreds of Myr.

• A key role for the photo-evaporation of the primordial hydrogen envelope in shaping the exoplanet population throughout the system lifetime. Super Earths and mini Neptunes (radii between 1 and 4 R_⊕ and masses lower than 20 M_⊕) at short orbital distances are more susceptible to lose part or all of their gaseous envelope (see e.g. Lopez et al. 2012), creating a gap between smaller super-Earth planets that have lost all the hydrogen/helium envelope, and the larger sub-Neptune planets that have retained some of it (e.g. Owen & Wu 2013).

• A difficult radial velocity (RV) characterisation since the host stars show an intense level of the stellar activity able to bury the small planetary signal (up to tens of times lower than the typical RV dispersion). The typically large rotation velocity of the star also degrades the quality of the spectrum in terms of RV precision and the use of very dense time sampling and sophisticated tools, necessary to model the stellar activity, is required (e.g. Damasso et al. 2020).

• A very low density and large scale height that make young planets, in principle, excellent targets for atmospheric characterisation; however, studies suggest the presence of high-altitude clouds, hazes and dust grains that increase atmospheric opacity, possibly leading to a null detection of atmospheric features (Wang & Dai 2019; Gao & Zhang 2020).

Benatti et al. (2019) and Newton et al. (2019) identified the first young planet detected by TESS in its first month of observations around the main component of the DS Tuc binary system in the Tuc-Hor association (G6 spectral type, visual magnitude = 8.2, estimated age ~40 Myr, TESS Object of Interest: TOI-200, Gaia DR3 source ID: 6387058411482257536). The planetary orbital period from the transit fit is 8.14 days and the planetary radius is about 0.5 R_J (5.6 R_⊕). The upper limit on the planetary mass is set to 1.3 M_J from Benatti et al. (2019) by using a few high-resolution archive spectra, and 0.09 M_J (28 M_⊕) from Newton et al. (2019) according to the mass-radius relation (Chen & Kipping 2017). Montet et al. (2020) and Zhou et al. (2020) found the planetary orbit well aligned with the equatorial plane of the host through the observation of the Rossiter-McLaughlin (RM) effect (Rossiter 1924; McLaughlin 1924) during the transit of DS Tuc A b.

Within the framework of coordinated observing programmes with the High Accuracy Radial velocity Planet Searcher (HARPS, Mayor et al. 2003) spectrograph at the 3.6 m-ESO telescope (La Silla Observatory, Chile), the Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO, Pepe et al. 2021) at VLT (Paranal Observatory, Chile), and XMM-Newton, we now aim to measure the mass of DS Tuc A b to confirm the orbital inclination, and to characterise its atmosphere. We show that the impact of the stellar activity on our data is particularly significant, but nevertheless, we were able to reduce the current value of the mass upper limit for the planet and to investigate the scenario of atmospheric escape by using mass and radius evolutionary models.

This paper is organised as follows. In Sect. 2, we describe the HARPS and ESPRESSO datasets. The RV modelling with the GP regression are provided in Sect. 3.2. The analysis of the RM effect from ESPRESSO data is reported in Sect. 4, while the results of the atmospheric characterisation are presented in Sect. 5. We present the results from XMM-Newton observations in Sect. 6 and describe our photo-evaporation models in Sect. 7. We provide a discussion in Sect. 8 and draw our conclusions in Sect. 9.

2 Datasets

In this section, we summarise the HARPS and ESPRESSO datasets collected for DS Tuc A.

2.1 HARPS–RV monitoring

We monitored DS Tuc A with HARPS at the 3.6 m telescope at the ESO-La Silla Observatory (Chile) within the framework of the HARPS time-sharing collaboration (programme ID 0103.C-0759(A), PI: Benatti). Thanks to this strategy, we benefited from the time shared among several observing programmes and gathered a time series composed by 76 spectra spanning from June to September 2019. As previously mentioned, DS Tuc A belongs to a visual binary system, with a projected separation from the stellar companion (K3 spectral type) of about 5′′. The presence of the companion does not affect our RV measurement, since the aperture of the HARPS fibre correspond to 1′′ on the sky. When possible, during the maximum visibility period of our target, we obtained two spectra in the same night with a few hours of separation. In doing so, we managed to collect nine couples of data points for a better knowledge of the intra-night RV scatter. Themean signal-to-noise ratio (S/N) of the spectra at 5500 Å is 108, with a typical exposure time of 900 s.

The HARPS spectra were reduced through the Data Reduction Software (DRS), which estimates the RVs by computing the Cross Correlation function (CCF) between the spectra and a digital mask depicting the features of a G2 star (seePepe et al. 2002 and references therein). In particular, since DS Tuc A is a moderately fast rotator (v sin i = 15.5 ± 1.5 km s⁻¹, Benatti et al. 2019) the CCF profile is broadened, so we adapted the width of the window for the CCF evaluation by using the offline version of the DRS. The RV dispersion of the whole time series and the mean RV uncertainty obtained with the HARPS DRS are 163 m s⁻¹ and 5.2 m s⁻¹, respectively.We also obtained an alternative processing of the RVs by using the pipeline TERRA (Anglada-Escudé & Butler 2012). In case of active stars, TERRA is more efficient in evaluating the RVs (see Carleo et al. 2020; Damasso et al. 2020)since it employs one of the observed spectra as a template of the line profiles, that is more suitable to describe potential distortions and asymmetries. For the TERRA dataset, the RV dispersion and the mean RV uncertainty are 138 and 5.4 m s⁻¹, respectively.The following analysis is based on the RV time series obtained with TERRA, reported in Table 1 with the flag ‘RV’ and displayed in Fig. 1.

Finally, we extracted the time series of several activity indicators such as the Ca II H+K lines and Hα by using the ACTIN code (Gomes da Silva et al. 2018), while the HARPS DRS provides the bisector span (BIS) of the CCF, which is useful for tracing the distortion of the line profiles produced by the stellar activity. These time series are reported in Table 1.

Fig. 1 HARPS RV time series of DS Tuc A as extracted from the TERRA pipeline. A bulk of data at JD = 8706 corresponds to the RVs collected during the planetary transit.

2.2 HARPS–transit

On the night of 2019-08-10, within the same observing programme described above, we observed the transit of DS Tuc A b in front of its star, collecting 30 consecutive data points. Due to the large airmass of the target (> 1.8) and the presenceof veils just before the event, we could not monitor the pre-transit phase. We collected 15 spectra during the transit and 15 spectra during the post-transit phase. To improve the temporal resolution, the exposure time was decreased to 600 s, with a mean S/N of 100, thanks to the good quality of the observing conditions during the event.

HARPS data were reduced as described in Sect. 2.1 and included in Table 1 with the flag ‘T’ (Transit). The RM effect detected with HARPS is depicted in Fig. 2 with red and green dots showing the comparison between the RV extraction obtained from Terra and DRS, respectively. Further comments on the figure are provided in the next section. The analysis of this transit is also reported in Montet et al. (2020), observed in that work with the PFS spectrograph, and Zhou et al. (2020), observed with the CHIRON spectrograph.

2.3 ESPRESSO

We observed a second transit of DS Tuc A b with ESPRESSO at VLT during the night of 2019-10-06 (programme 0104.C-0798, PI Borsa). We obtained 97 spectra with 300 s of exposure time that produced spectra with mean S/N of 111 at 5500 Å. The mean uncertainty of the radial velocity is 2.5 m s⁻¹. We collected 22 spectra before the transit event, 32 during the transit and 45 spectra were finally collected after the transit, for a total coverage of about ten hours. The object airmass at the beginning of the observation was 1.74, reaching the minimum value of 1.40 during the transit and 2.52 at the end of the night. This transit is also reported in Zhou et al. (2020), observed with the PFS spectrograph.

The ESPRESSO RVs were obtained from the ESPRESSO DRS (1.3.2), also based on the CCF method, by using the available G8 mask as a spectral template with the offline version¹ and a largerwidth of the window for the evaluation of the CCF. The ESPRESSO RVs are reported in Table 2 and are shown in Fig. 2 with blue dots as a function of time, compared with those obtained with HARPS on Aug 10. To obtain Fig. 2, we applied an offset of approximately − 25 and − 105 m s⁻¹ for the HARPS RVs obtained with TERRA and DRS, respectively, and approximately − 33 m s⁻¹ for the ESPRESSO RVs. The RM amplitude is is slightly lower than 50 m s⁻¹ and it is compatible for the two series of data (and consistent between the two different RV extractions for HARPS data). The out-of-transit series show a clear increasing trend which is in contrast with the expectation from Keplerian motion and is dominated by the stellar activity. Interestingly, for different transit events the slope of the HARPS out-of-transit series, as extracted from the official DRS (green dots), is the same as the one shown by ESPRESSO (blue dots). The two series deviate at about 1.5 hr after the egress, when the ESPRESSO time series shows an RV‘bump’ (investigated in Sect. 4). Instead, the out-of-transit series obtained with TERRA (red dots) shows a lower steepness after the fourth contact with respect to the DRS series. This demonstrates that the usage of TERRA provides RV estimations that are less affected by stellar activity since the disagreement between TERRA and DRS occurs when the RV are dominated by the activity (out of the transit) and not by the planet contribution (during the transit).

Finally, the comparison between the RM effect for the transit of Aug. 10 detected by HARPS and PFS (reported by Montet et al. 2020)shows that the out-of-transit modulation is actually due to the stellar activity, since the data from two different instruments with different RV extraction (TERRA vs. Iodine cell technique) are in full agreement. No comparison can be performed with the transits reported in Zhou et al. (2020), namely Aug. 10 with CHIRON and Aug. 18 and Oct. 6 with PFS since the data shown in their Fig. 3 are already corrected for the activity and have a shorter out-of-transit baseline.

Fig. 2 RM effect for the transits observed with HARPS (red dots, TERRA; green dots, DRS, systemic RV removed) and ESPRESSO (blue dots, DRS, systemic RV removed). Vertical lines indicate the first and fourth contact (dotted) and the centre of the transit (dashed), according to the ephemeris in Newton et al. (2019).

Fig. 3 GLS of the HARPS RV (black line) and BIS (red line) time series of DS Tuc A. The green vertical line indicates the locationof the planet orbital period, while the window function of the observations is reported in the inset plot. Dotted, dashed, and solid black horizontal lines indicate the FAP thresholds at 1, 0.1, and 0.01%, respectively, for the RV periodogram.

3 RV modelling and mass upper limit of DS Tuc A b

In this section, we present our RV analysis aiming to recover the planetary signal of DS Tuc A b by using the Gaussian process regression. Our results are subsequently compared with injection-retrieval simulations.

3.1 Frequency analysis

Figure 3 compares the Generalised Lomb-Scargle periodogram (GLS, Zechmeister et al. 2009) of the HARPS RVs (black line) and BIS (red line) time series. The stellar activity dominates the behaviour of the RVs since the two GLS periodograms show a similar shape with a significant periodicity around three days (False Alarm Probability lower than 0.01%, evaluated with 10 000 bootstrap random permutations), which is compatible with the rotation period of the star (2.85 days, see Benatti et al. 2019). The presence of additional peaks between 0.29 and 0.4 d⁻¹ is mainly due to the spectral window of the observations (reported in the inset plot). The expected location of the orbital period of DS Tuc A b is marked with a green vertical line. No evidence of power is present around that frequency value.

The GLS of the other activity indicators extracted from TERRA show the similar periodic content, but only the BIS shows a strong correlation with the RVs (Spearman coefficient close to 0.9). The decorrelation of the RVs with the values of the BIS produces a time series with dispersion of about 90 m s⁻¹. The GLS of the residuals does not show any periodicity at the planet period, only some additional power at around three days, suggesting that the RV contribution of the stellar activity must be modelled in a more sophisticated manner.

3.2 Gaussian process regression

3.2.1 Quasi-periodic kernel

We used a Gaussian process (GP) regression to model the activity-correlated signal in the RVs. The fitting analysis was carried out using the publicly available Monte Carlo nested sampler and Bayesian inference tool MultiNestv3.10 (e.g. Feroz et al. 2019), through the pyMultiNest wrapper (Buchner et al. 2014). The algorithm was setup to run with 1000 live points and a sampling efficiency of 0.3. We fitted the combined HARPS+ESPRESSO dataset, excluding the data taken during the transits. The zero-point corrections between the two datasets (γ_HARPS and γ_ESPRESSO in Table 3) are also fitted in the process. In what follows, we modelled the planetary signal with a sinusoid, assuming a circular orbit as in Benatti et al. (2019). To confirm this assumption, we performed a joint RV+TESS photometry Monte-Carlo analysis, using the same transit light curves analysed in Benatti et al. (2019). We found that the use of the RVs does not help in constraining the eccentricity, which turns out to be compatible with zero within ~1.5σ ($e_{\text{b}}={0.14}_{-0.10}^{+0.29}$). Hereafter, we present results from the analysis of only RVs, adopting the planet radius given by Benatti et al. (2019) and assuming, as the priors for the transit ephemeris, the more precise values of P_b and T_{b, conj.} obtained by Newton et al. (2019). Initially, theactivity term was modelled by adopting the quasi-periodic covariance matrix to model the correlated signal due to stellar activity, defined as follows: $$\begin{array}{lll}k(t,t^{\prime })\hfill & =\hfill & h^2\cdot \exp \left[-\frac{{(t-t^{\prime })}^2}{2{\lambda }^2}-\frac{{\text{sin}}^2(\pi (t-t^{\prime })/\theta )}{2w^2}\right]\hfill \\ \hfill & \hfill & +({\sigma }_{\text{RV}}^2(t)+{\sigma }_{\text{jit}}^2)\cdot {\delta }_{t,t^{\prime }}.\hfill \end{array}$$(1)

Here, t and t′ represent two different epochs of observations, σ_RV is the radial velocity uncertainty, and ${\delta }_{t,t^{\prime }}$ is the Kronecker delta. Our analysis takes into account other sources of uncorrelated noise – instrumental or astrophysical (or both) – by including a constant jitter term σ_jit, which is added in quadrature to the formal uncertainties σ_RV. The GP hyper-parameters are h, which denotes the scale amplitude of the correlated signal; θ, which represents the periodic time-scale of the modelled signal, and corresponds to the stellar rotation period; w, which describes the ‘weight’ of the rotation period harmonic content within a complete stellar rotation (i.e. a low value of w indicates that the periodic variations contain a significant contribution from the harmonics of the rotation periods); and λ, which represents the decay timescale of the correlations, and is related to the temporal evolution of the magnetically active regions responsible for the correlated signal observed in the RVs.

The GP regression and kernel implementation was performed using the publicly available PYTHON module GEORGE v0.2.1 (Ambikasaran et al. 2015). which was integrated within the nested sampling algorithm. A summary of priors and results is given in Table 3 and the posterior distributions for all the free parameters are shown in Fig. 4. The fitted planetary Doppler semi-amplitude K_b is consistent with zero, therefore, we can only provide an upper limit for the planet mass M_b] < 14.4 M_⊕ at 1σ level. Figure 5 shows the RV residuals (best-fit stellar signal and instrumental offsets subtracted from the original RVs) phase-folded to the orbital period of DS Tuc A b. According to our results, the evolutionary time scale λ of the quasi-periodic activity component is of the order of twice the rotation period. The quick evolution of the correlated activity signal is confirmed by analysing the Hα and the S-index activity diagnostics, calculated with TERRA from the HARPS spectra, with the GP quasi-periodic kernel. By adopting for the hyper-parameters the same priors used for the RVs (except for the amplitude h), we got ${\theta }_{\text{H}\alpha }={1.499}_{-0.008}^{+1.488}$ d and λ_Hα = 6.8 ± 1.2 d, ${\theta }_{\text{S-index}}={2.9}_{-1.4}^{+0.1}$ d and ${\lambda }_{\text{S-index}}={3.7}_{-1.0}^{+1.6}$ d. The posterior distribution of the rotation period for both activity indexes shows an over-density at the location of the first harmonic, as indicated by the asymmetric error bars, which is not found for the RVs. This evidence suggested a test with a different kernel for modelling the activity signal in the RVs, presented in the following subsection.

3.2.2 Simple harmonic oscillator kernel

For this test, we considered a mixture of two stochastically driven, damped simple harmonic oscillators (SHOs) with undamped periods of P_rot and P_rot/2. This can be done by combining two SHOTerm kernels included in the package celerite (Foreman-Mackey et al. 2017) as done, for instance, by David et al. (2019) to model the Kepler/K2 photometric variability of another very young star, V1298 Tau. The power spectral density corresponding to this kernel (also known as the ‘rotation kernel’) is: $$S(\omega )=\sqrt{\frac{2}{\pi }}\frac{S_1{\omega }_1^4}{{({\omega }^2-{\omega }_1^2)}^2+2{\omega }_1^2{\omega }^2}+\sqrt{\frac{2}{\pi }}\frac{S_2{\omega }_2^4}{{({\omega }^2-{\omega }_2^2)}^2+2{\omega }_2^2{\omega }^2/Q^2},$$(2)

where $$S_{\text{1}}=\frac{A}{{\omega }_{\text{1}}{Q}_{\text{1}}},$$(3) $$S_{\text{2}}=\frac{A}{{\omega }_{\text{2}}{Q}_{\text{2}}}\times f,$$(4) $${\omega }_{\text{1}}=\frac{4\pi Q_{\text{1}}}{P_{\text{rot}}\sqrt{4Q_{\text{1}}^2-1}},$$(5) $${\omega }_{\text{2}}=\frac{8\pi Q_{\text{2}}}{P_{\text{rot}}\sqrt{4Q_{\text{2}}^2-1}},$$(6) $$Q_{\text{1}}=\frac{1}{2}+Q_{\text{0}}+\Delta Q,$$(7) $$Q_{\text{2}}=\frac{1}{2}+Q_{\text{0}}.$$(8)

The parameters in Eqs. (3)–(8), where the subscripts 1 and 2 refer to the primary (P_rot) and secondary (P_rot/2) modes, are the inputs to the SHOTerm kernels. However, instead of using them directly, we adopt a different parametrisation using the following variables as free hyper-parameters in the MC analysis, from which S_i, Q_i, and ω_i are derived through Eqs. (2)–(7): the variability amplitude A (m² s⁻²), the stellar rotation period, P_rot (the periodicity of the primary term), the quality factor, Q₀, the difference ΔQ between thequality factors of the first and second terms, and the fractional amplitude, f, of the secondary mode relative to the primary². We adoptedpriors in the natural logarithm space for all the hyper-parameters, except for P_rot. Priors and results of the fit are shown in Table 3 and in the second panel of Fig. 4. Even in this case, we are able to provide only an upper limit for the planet mass (M_b] < 40.3 M_⊕ at 1σ level), higher than that found using a quasi-periodic kernel. The instrumental jitters are lower for both HARPS and ESPRESSO, denoting a better fit to the data than the quasi-periodic activity model. However, the model with the mixture of SHOs is statistically mildly favored over the quasi-periodic ($\Delta \ln \mathcal{Z}$ = 3.4), therefore, we do not have a strong argument to prefer this model over the other. Figure 6 shows a comparison between the best-fit activity model obtained with the quasi-periodic and rotation kernel from the HARPS+ESPRESSO joint analysis, for a subset of the HARPS RVs.

Fig. 4 Corner plots for GP quasi-periodic (upper panel) and SHO mixture (lower panel) models used to fit the HARPS+ESPRESSO radial velocities.

3.3 Testing the mass upper limit of DS Tuc A b through injection and recovery simulations

The main conclusions from the RV analysis are that we are able to provide only an upper limit for the mass of DS Tuc A b and this limit changes depending on how the correlated signal, due to the stellar activity, is modelled. It may be that the GPquasi-periodic kernel is absorbing part of the planetary signal, causing the mass upper limit to be nearly threetimes lower than that found using the rotation kernel. This result has important consequences on characterising the planetary structure and studying the evolution of the planet atmosphere within the first Myrs from the formation, based on models of photo-evaporation that require the current mass as input. Statistical considerations suggest that slight preference should be given to the results obtained with the adoption of the rotation kernel, therefore, to the higher mass limit that we found for DS Tuc A b, even if there is not a strong evidence in favor of it. Based on a similar statistical argument, the rotation kernel describes the time series of the Hα and CaII H&K activity diagnostics better than the quasi-periodic kernel. On the other hand, both kernels represent a priori a proper choice to treat the activity-related RV signals for a very active star such as DS Tuc A, and still very little is known in literature about their relative performances for this kind of stars. Since it is difficult to drawn general conclusions about this issue, tests should be done on a star-by-star basis. Therefore, we devised two different sets of injection-retrieval simulations, provided in the next subsections, to get further information on which kernel provides more reliable results for the specific case of DS Tuc A.

3.3.1 Injection and recovery test performed directly on the real data

First, we injected a planetary signal directly in the real data, namely using the same time stamps of the original RVs. This signal is represented by a sinusoid with period and phase fixed to the best-fit values found for DS Tuc A b from the analysis of the photometric light curve, and adopting a grid of values for the semi-amplitude K_b in the range 0.75–40 m s⁻¹. Then, each simulated dataset was fitted exactly within the same framework as for the real data, testing both the quasi-periodic and rotation kernels, and we compared the injected K_{b, inj} with the retrieved K_{b, retrieved}, which is expressed as the 50th percentile of the posterior, as usual. This test can help to understand which model recovers better the semi-amplitude of the planetary signal. Since the transit ephemeris are very precise, the injection of the Doppler planetary signal should only amplify what is already hidden in the data that we could not constrain by fitting the original RVs. The results of this set of simulations are shown in Table 4, and the posteriors of K_{b, retrieved} are shown in Fig. 7. There are some interesting features to point out here. Concerning the quasi-periodic model, K_{b, retrieved} is only slightly overestimated and accurate within less than 1 m s⁻¹ for K_{b, inj} ≥ 10 m s⁻¹, with a significance ranging from 2.4σ to nearly 8.2σ, denoting a successful recovery of the injected signal. The overestimate does not exceeds 3 m s⁻¹ for the lower injected signal K_{b, inj} = 0.75 m s⁻¹. For the case of the rotation/SHO mixture kernel, the picture is different. K_{b, retrieved} is surprisingly underestimated for K_{b, inj} >15 m s⁻¹, with the highest departure from the injected value occurring for K_{b, inj} = 40 m s⁻¹, which is supposed to be the easiest case for recovering the actual signal. While the overestimates observed for the quasi-periodic model is an expected result since the injected sinusoid adds up with the real signal due to planet b already present in the data, we interpret the underestimate observed for the case of the rotation kernel as a general failure of the model for values of K_{b, inj} sufficiently high. Further, in this case the error bars associated to K_{b, retrieved} are larger than for the quasi-periodic model, due to the density profile of the posteriors which follow an asymmetric distribution. For K_b,inj ≤ 10 m s⁻¹, the recovered semi-amplitude is overestimated by an amount larger than for the quasi-periodic model. We remark that the best-fit values of the GP hyper-parameters derived from each mock dataset perfectly agree with those determined for the original RVs, indicating that the recovered stellar activity signals were not altered by the injection of a sinusoid into the data.

The overall picture emerging from the first set of simulations lead us to conclude that, contrary to what suggested by the statistical evidence, the quasi-periodic model for the activity signal provides more reliable results and guarantees the best performancefor DS Tuc A. This in turn makes us more confident that the mass of DS Tuc A b should be not larger than 14.4 M_⊕ with 68% of confidence.

Fig. 5 Spectroscopic orbit of planet DS Tuc A b. The instrumental offsets and the best-fit stellar activity signal (GP quasi-periodic model) have been subtracted from the HARPS and ESPRESSO original RVs. The error bars include the uncorrelated jitter terms added in quadrature to the formal uncertainties. The curve in red represents the best-fit Doppler signal (circular obit; Kb = 3.2 m s−1).

Fig. 6 Summary of the GP best-fit model as obtained from the joint analysis of the HARPS+ESPRESSO RVs. Here, we show a portion of the HARPS time series, after removing the best-fit signal induced by planet b. The shaded area in gray represents the ±1σ confidence interval. Upper panel: best-fit model for the GP quasi-periodic kernel (orange line). Lower panel: same as above, but for the GP rotation kernel (blue line).

3.3.2 Injection and recovery test using mock RVs stemming from the best-fit solutions found for the real data

In order to test further this conclusion, we devised a second type of simulations. We generated ten mock RV time series using the best-fit solutions of the quasi-periodic and SHO mixture activity models, by randomly drawing the hyper-parameters from distributions based on the results in Table 3 and randomly sampling a list of predictions for the RV correlated signal at thereal observing epochs using the sample module in george and celerite. We then injected a planetary sinusoid with K_b = 40 m s⁻¹, P_b = 8.138268 d, and T_{b, conj} = 2 458 332.30997 BJD into the mock activity signals, including the best-fit RV instrumental zero-points γ_HARPS and γ_ESPRESSO in Table 3 (the same constant values are adopted for all the mock dataset). Finally, we added to the RVs calculated so far a term (positive or negative) to take the instrumental uncorrelated jitters into account. For each time stamp, this term is randomly drawn from a normal distribution $\mathcal{N}(0,{\sigma }_{\text{jitt}}^2)$ (for both HARPS and ESPRESSO), which corresponds to the operation σ_jitt × numpy.random.randn(N) using the syntax of python, where N is the number of observing epochs. For each mock dataset, the jitter terms σ_jitt are randomly drawn from the posterior distributions derived for the real data (Table 3). For the mock RVs we used the same error bars of the real data, and the instrumental zero-points and uncorrelated jitters were treated as free parameters in the Monte Carlo analysis, generally adopting the same framework for the analysis of the real RVs. We fitted each mock RV dataset using both the quasi-periodic and rotation GP models. This gives a total of four tests: (a) generated: quasi-periodic; fitted: quasi-periodic; (b) generated: quasi-periodic; fitted: rotation; (c) generated: rotation; fitted: quasi-periodic; (d) generated: rotation; fitted: rotation. The idea behind this test is to evaluate the general performance of the two models in retrieving the Doppler planetary signal with a cross check, taking advantage of the fact that we know the underlying (mock) stellar activity signal in each dataset.

For illustrative purposes, part of the ten posteriors of K_b obtained for each test is shown in Fig. 8. It can be seen that the injected planet is well recovered only for the case a (K_{b, retrieved} = 38.3 ± 2.9 m s⁻¹) with a median detection significance of 5.1σ over the ten mock dataset. For the case b, the semi-amplitude is systematically underestimated (K_{b, retrieved} = 32.2 ± 2.9 m s⁻¹) and the posteriors are wide, with a median detection significance of 2.1σ. These results are consistent with those in Table 4 for K_{b, inj} = 40 m s⁻¹ obtained using the real data. The fitting procedure equally fails in recovering the planetary signal when the mock dataset are generated using the SHO mixture kernel (K_{b, retrieved} = 29.9 ± 9.4 m s⁻¹ and K_{b, retrieved} = 33.1 ± 9.1 m s⁻¹, for case c and d, respectively). The result of this cross-check further support the conclusion that using a quasi-periodic kernel to model the stellar activity better describes the temporal variations observed in the real RV data.

Fig. 7 Posteriors of the Kb parameterretrieved from the simulations involving the original RV datasets, with the injected Kb, inj values indicated in the legend. The vertical dashed lines represent the median of the posteriors of the corresponding colour. The vertical black lines indicates the median of the posteriors for the real data (3.2 and 8.6 m s−1), as a reference.

4 RM modelling

The employment of the GP regression allowed us not only to model the stellar activity related to the time series, but also to make predictions on the activity behaviour during the transit on the night of Oct 6, observed with ESPRESSO, for which we have a dense monitoring session before and after the event. The activity model obtained with the SHO kernel was more suitable than the quasi-periodic kernel to describe the out-of-transit modulation. In particular, the quasi-periodic kernel failed in the modelling of the RV ‘bump’ occurring at about JD2458763.76 (see Fig. 9, upper panel), while this feature was better treated by the SHO. It is possible that the quasi-periodic kernel is more adapt to describe the global behaviour of the activity,while the SHO better manages events on a short-term basis. The SHO model is shown in the upper panel of Fig. 9 (orange line), where the residuals are also reported. The RV bump seems to be connected with a strong activity event like a stellar flare, as confirmed by the time series of the CaII H&K and Hα indices, obtained with the ACTIN code (and reported in Table 2). At the time corresponding to the RV bump, a sudden increase of the value of both indices can be observed (lower panel of Fig. 9), followed by a second peak and a slow decay, possibly down to the value before this event.

Since the HARPS monitoring missed the pre-transit-ingress phase, it was not possible to reconstruct a similar activity model for the transit of Aug. 10. For this reason we only considered the ESPRESSO time series to fit the RM effect, showing higher temporal resolution and lower RV uncertainties. However, we stress that HARPS data allowed us to confirm the amplitude of the RM effect detected with ESPRESSO.

We fitted a RM model (Ohta et al. 2005) to the activity-corrected RVs in a Bayesian framework by employing a differential evolution Markov chain Monte Carlo (DE-MCMC) technique (Ter Braak 2006; Eastman et al. 2013), running ten DE-MCMC chains of 100,000 steps and discarding the burn-in. The medians and the 15.86 and 84.14% quantiles of the posterior distributions were taken as the best values and 1σ uncertainties. The first six in-transit RV measurements were masked from the fit procedure, as they were clearly deviating from the other RVs and could have biased the retrieved value of the limb darkening parameter μ. We took into account the error given by the stellar activity correction performed by adding it in quadrature to the one of the RVs extracted by the DRS pipeline. Gaussian priors for i, a/Rs, and Rp/Rs were set to the values in Benatti et al. (2019), leaving, as free parameters, the linear limb darkening coefficient μ₁, v sin i and the projected spin-orbit angle λ_PSO (we use this notation to distinguish this parameter from the GP hyper-parameter λ used in Sect. 3.2). The v sin i was left free because different RM models are known to produce vsini measurements that are often in disagreement with each other and with estimates obtained from spectral line broadening (Brown et al. 2017), while for the linear limb darkening parameter μ, we left it free because it is sensitive to the stellar activity and often the theoretical values are in disagreement with the observed ones for active and spotted stars (e.g. Csizmadia et al. 2013).

The RVs and the best fit are shown in Fig. 10, with the best-fit values presented in Table 5. We note that μ takes a value which is quite distant from the theoretical one calculated for the star on the ESPRESSO wavelength range using LDTK (μ = 0.704 ± 0.002, Parviainen & Aigrain 2015), but still compatible within 3σ. When setting a prior to the LDTK estimated value, the returned fit is clearly not matching the observations even by eye, but we note that the resulting value of λ_PSO is perfectly compatible the one here reported.

The measurement of the projected spin-orbit obliquity, λ_PSO, is in well agreement with the one of Montet et al. (2020, 12 ± 13 degrees), thus confirming an almost (but not perfectly) aligned system. Zhou et al. (2020) reported a value of λ_PSO of 2.5 ± 1 deg, which is not compatible with our value. Indeed, the determination of λ_PSO through the RM effect is sensitive to stellar activity (Oshagh et al. 2018), so it is possible (and likely) that the differences are due to the activity correction performed. While in the RM fit procedure Montet et al. (2020) and Zhou et al. (2020) use different spot-modeling techniques, we let all the residual activity left from the Gaussian process model (which does not take into account the in-transit RV values) to be included in the limb-darkening estimation. A common analysis of all available datasets, which are taken with different instruments and have different quality, will help in check the stellar activity impact on the retrieved λ_PSO.

We also performed a line profile tomography aiming to reveal changes in the CCF morphology before, during, and after the transit of DS Tuc A b. As in the case of AU Mic (Fig. 4 in Palle et al. 2020), our Doppler tomography shows the very high activity level of the host (Fig. 11), which prevented us from obtaining a reference CCF for a proper subtraction and, in turn, to fit the planetary Doppler shadow. However, this feature can be noticed as a slight ripple occurring during the transit phase, compatible with an aligned planet. The flare event previously identified can be observed also in Fig. 11 as a clear increase of the distortion of the line profile (depicted in dark red) at phase ~ 0.015. We note that methods such as the reloaded RM (Cegla et al. 2016) which, in principle, are less affected by stellar activity in determining the projected spin-orbit angle, would also likely fail with such a level of activity.

Fig. 8 Posteriors of the Kb parameter retrieved from the simulations involving mock RV datasets build on the best-fit stellar activity signals based on a quasi-periodic and SHO mixture models (see the text for details). The simulated RV time series include an injected planetary signal of semi-amplitude Kb, inj = 40 m s−1. Only five out of a total of ten posteriors are shown for the sake of a better the readability First row: results from the fit of data generated using a quasi-periodic kernel. Second row: results from the fit of data generated using a SHO mixture kernel. In both rows, each colour identifies one mock dataset, therefore, the posterior distributions with the same colour can be compared. The vertical lines indicate the median of the posteriors.

Fig. 9 Radial velocities and activity indices variation during the transit of Oct. 6 of DS Tuc A b observed with ESPRESSO. Upper panel: RV (blue dots) with the GP model of the stellar activity (orange line) and the corresponding residuals. Lower panel: CaII H&K and Hα time series extracted from the same spectra.

Fig. 10 RM effect for DS Tuc A b from ESPRESSO data modelled by our fit (red line) with the corresponding residuals. Vertical dotted lines indicate the transit ingress and egress.

Fig. 11 CCF residuals during the transit of DS Tuc A b in the stellar restframe. Horizontal white lines indicate the transit ingress and egress.

5 Atmospheric characterisation

In the following sections, we provide our analysis of the ESPRESSO spectra taken during the transit of DS Tuc A b occurred on Oct. 6, 2019, aiming to detect traces of its atmospheric compounds.

5.1 Transmission spectrum

We extracted the planetary transmission spectrum from the ESPRESSO data following the classic procedure used in high-resolution transmission spectroscopy (e.g. Wyttenbach et al. 2015; Borsa et al. 2021). We shifted the spectra to the stellar restframe by using the Keplerian model of the system with the upper limits found in Table 3, then we normalised each spectrum. Telluric correction was performed by exploiting the scaling relation between airmass and telluric line strength (e.g. Snellen et al. 2008; Astudillo-Defru & Rojo 2013), rescaling all the spectra as if they were observed at the airmass of the transit center. Then we divided all the spectra by a master stellar spectrum created by averaging the out-of-transit spectra. To avoid being too much contaminated by stellar activity, we used only the out-of-transit data close to the transit, this means that we discarded the first 10 and the last 33 spectra of the ESPRESSO timeseries. The ‘wiggles’ pattern (e.g. Tabernero et al. 2021) is corrected in the same way presented in Borsa et al. (2021), namely by fitting a sinusoid with varying period, phase and amplitude for each spectrum, with the fit performed independently for each wavelength region where we looked for spectral features. At the end, we moved to the planetary restframe by shifting all the residual spectra for the theoretical planetary radial velocity. The planetary signal during the transit is always within the stellar line profile (Fig. 12). The transmission spectrum was then created by averaging all the full-in-transit residual spectra (see e.g. Fig. 13 for Hα, Hβ, Na D2, and Li lines).

Since this star has a very high level of activity, its line profile changes rapidly, thus the consequence is that the stellar activity deeply contaminates the transmission spectrum and prevents us from probing the planetary atmosphere as shown by Fig. 13. This was also noted for the case of AU Mic b (Palle et al. 2020), another planet transiting a young star. The movement of active regions on the star is sufficiently fast to give well-detectable changes in the line profiles during the time span of our observations. This is fully visible in Fig. 11, where these active regions are the blue and red stripes (flux decrease and increase with respect to the average out-of-transit CCF, respectively) moving along the residuals of the CCF, which can be considered as an average stellar line profile. The slope of these stripes in the velocity-phase 2D plot (~4–5 km s⁻¹ in ~4 h) is compatible with regions moving on the stellar surface when assuming a v sini ~ 20 km s⁻¹ and stellar rotation of 2.9 days (Table 3). In order to extract the transmission spectrum, we are forced to set some baseline out-of-transit, but we could not find any combination of spectra for which the out-of versus in-transit active regions configuration is the same. Interestingly, the effect of the contamination in the transmission spectrum looks wavelength dependent (Fig. 13), being more prominent toward blue wavelengths and going to almost disappear toward the red edge of the spectrograph. We expect that features like starspots are the main responsible for this effect since the temperature difference with the photosphere (500–2000 K) is much larger than the one predicted with the faculae (~ 100 K in the solar case).

We searched for possible absorption given by Hα (R_eff < 1.22 R_p), Hβ (R_eff < 1.3 R_p), Na D doublet (R_eff < 1.28 R_p), K (R_eff < 1.16 R_p), Mg (R_eff < 1.21 R_p), Li (R_eff < 1.12 R_p), without finding any trace of significant absorption beyond the variation caused by stellar activity. We estimated indicative upper limits to the effective planetary radius R_eff by using the standard deviation of the transmission spectrum in a range of ± 1 Å, centered at the rest frame wavelength of the element. We translated this value to R_eff by assuming $R_{\text{eff}}^2/R_{\text{p}}^2=(\delta +h)/\delta$, where δ is the transit depth and h the calculated standard deviation.

Fig. 12 2D tomography of Hα in the stellar restframe. The vertical dotted line shows the stellar restframe. The diagonal magenta dashed line shows the expected planetary restframe. The horizontal white lines mark the beginning and end of the transit.

5.2 Cross-correlation with theoretical templates

In order to try to overcome the stellar variability issue, we tried to look for elements or molecules that we expect to find possibly more in the planetary atmosphere than in the star. We generated model templates for V, Cr, AlO, VO, TiO (Plez linelists), and Y by using petitRADTRANS (Mollière et al. 2019), assuming an isothermal atmospheric profile with T = 1000 K, continuum pressure level of 10 mbar, and solar abundances. The atmospheric models were then translated into flux (R_p /R_s)², convolved at the ESPRESSO resolving power and continuum normalised. Cross-correlation between model templates and the single transmission spectra was then performed in the same way as in Borsa et al. (2021). Unfortunately, we find no evidence in the planetary atmosphere for the presence of any of the species we were looking for.

6 X-ray emission

In this section, we report our analysis of the XMM-Newton observations of the DS Tuc binary system.

6.1 XMM-Newton data

We observed DS Tuc with XMM-Newton (obsid:0863400901, PI S. Wolk), aiming to assess the spectrum and time variability of the host star in X-rays and, hence, to characterise the coronal activity level and the dose of high-energy radiation received by its young planet. A measure of the strong irradiation by the host star can give us useful constraints with which to model the evolution of the planetary atmosphere due to photo-evaporation processes.

The exposure time of the observation was about 40 ks and the prime instrument was EPIC with full frame window imaging mode and the Medium filter. The Observation Data Files were retrieved from the XMM-Newton archive and reduced with the Science Analysis System (SAS, ver.18.0.0³). We obtained FITS lists of photons detected with the MOS1, MOS2, and pn CCD cameras, calibrated in energy, arrival time and astrometry with the evselect SAS task. Since the two components of the DS Tuc system are separated by about 4.8′′, they are partially resolved in the MOS images (2′′ pixel resolution), but not with the pn CCD (about 4′′ pixel resolution). Figure 14 shows the MOS1 image of DS Tuc, with the system primary (A) being the source at the bottom, and the secondary (B) at the top. The two sources appear with similar intensities, but the tight separation implies a cross-contamination of order of 16% in estimating their fluxes.

In order to estimate the contamination and extract the correct flux for each component, we compared the events extracted from two different sets of regions. The first set is made of two semicircular regions with radius 30′′ centered half way between the two sources. These regions collect most of the counts of the two sources, but they suffer from a high cross-talk contamination. The second set is made of two smaller circles of 10′′ radius, displaced by ~9′′ from each emission peak along the direction connecting the pair. In this way, we minimise the cross-talk contamination still retaining enough counts for light curves and spectra.

The extracted light curves are shown in Fig. 15. These plots show that the star A is less variable and slightly brighter than the star B, with the latter characterised by at least a couple of flares of duration 7–10 ks that contaminate the light curve of A accumulated in the large semicircular region.

From the events in the 10′′ regions, we extracted spectra and relative instrument response files with SAS. The spectra for the A and B sources are shown in Fig. 16. We fitted the spectra with XSPEC, adopting a multi-temperature optically thin coronalemission model (APEC) multiplied by a global interstellar absorption component (TBABS). The results are reported in Table 6. In order to obtain a statistically acceptable fit, we need to describe thecorona of each source with at least three isothermal components: a cool temperature component at ~ 0.3 keV, another component at ~0.95 keV and a hot component at ~2−5.5 keV. From the estimates of the plasma emission measure, the 0.95 keV component is the most dominant, followed by the cool component at 0.3 keV, and then the hot component, which is less constrained.

Fig. 13 Transmission spectrum extracted in the zone of the Hβ (left panel), Na D2 (central left panel), Hα (central right panel), and Li (right panel) lines. The vertical dashed lines show the planetary restframe. Black points refer to 0.1 Angstrom binning. Different lines we searched for are affected by the contamination given by the stellar activity in different ways.

Fig. 14 XMM/Newton MOS1 image of DS Tuc. The two red points indicate the ICRS positions of the two stars from SIMBAD, the centroid of X-ray emission is offset because of the proper motions of the system. The contours show the emission of the two sources, marginally resolved. The large circle (in red) has 30′′ radius, while the smaller circles (in yellow) are 10′′ wide. They were used to extract the spectra and light curves of the A and B system components (see text).

6.2 Results

We evaluated the luminosity for both sources also in the ROSAT band (0.1–2.4 keV) for a direct comparison with the previous measurements. In this band, the luminosity for the primary (A) is 9.57 × 10²⁹ erg s⁻¹, and 7.27 × 10²⁹ erg s⁻¹ for the secondary (B). Since the X-ray luminosity evolution is usually calibrated in the ROSAT band (e.g. Penz et al. 2008), we considered this value for DS Tuc A to test photo-evaporation models of the planetary atmosphere (see the next section).

The gas absorption is better constrained from the spectrum of DS Tuc A, with a best-fit interstellar hydrogen column density N_H ~ 2.4 × 10²⁰ cm⁻², which is compatible with the low value of optical A_V = 0.15 from the Gaia DR1 data analysis by Gontcharov & Mosenkov (2018), typical for nearby stars within ~ 50 pc.

We verified that the total observed X-ray flux of the combined A+B source is about 70% the one derived from the detection reported in the ROSAT All-Sky Survey Point Source Catalogue (2RXS), assuming the same three-component thermal model best-fitting the XMM-Newton data. This result implies very little variability on a time scale of about 30 yr, given that DS Tuc A is an active star, and this is comparable to the short-term variability that we observed in the time frame of the XMM-Newton observation (Δt ~ 0.5 d). Since we have now resolved for the first time the two stellar components, we can recompute more precisely the ratios of X-ray to bolometric luminosities: assuming L_bol,A = 2.80 × 10³³ erg s⁻¹ and L_bol,B = 1.26 × 10³³ erg s⁻¹ (Newton et al. 2019), we derive logL_x∕L_bol = −3.55 and − 3.31 for the G6V primary and the K3V secondary, respectively. This result implies that DS Tuc A shows X-ray emission just below the threshold activity level for saturated coronal sources (Pizzolato et al. 2003), that solar-type stars maintain until their rotation period remains below ~2 days. In fact, P_rot,A = 2.85 ± 0.02 d (Benatti et al. 2019). Instead DS Tuc B is likely at the saturation level, possibly due to its lower mass, although we lack a measure of P_rot,B. Finally, we have checked that the X-ray luminosities of both components are consistent with the expectation for stars in open clusters of similar age, such as the 50 Myr old cluster of α Per (Pillitteri et al. 2013). These results allow us to model the evolution of the average high-energy irradiation received by the planet, and to estimate its atmospheric evaporation in the last few tens of millions of years, that is, since the dissipation of the parent circumstellar disc.

7 Atmospheric evolution induced by photo-evaporation

In Sect. 3.2, we obtained, for DS Tuc A b, an upper limit of 14.4 M_⊕ on the mass. In this section, we investigate the predicted atmospheric evolution induced by the stellar high energy radiation derived in Sect. 6.2 (hereafter XUV) for different values of the planetary mass, in order to evaluate the different scenarios.

Fig. 15 Mos1 light curves for both components of DS Tuc. Black line is the light curve from the semi-circular region of radius 30′′, red line is the light curve from the circular region of radius 10′′. Error bars on rates are at 68% confidence level. Upper panel: DS Tuc A. Lower panel: DS Tuc B.

7.1 Set-up of the models

Proceeding in a similar fashion as in Carleo et al. (2021) and Georgieva et al. (2021), we evaluated the mass loss rate of the planetary atmosphere using the hydrodynamic-based approximation developed by Kubyshkina et al. (2018). For the X-ray evolution, we use two different laws both based on cluster and field star data, namely, those proposed by Penz et al. (2008, hereafter P08a), and Jackson et al. (2012, hereafter J12). Both describe the evolution of the X-ray luminosity, L_x, with a saturation and a decay phase. The main difference is that J12 assume a constant X-ray emission up to an age of ~ 70 Myr, while P08a predict high-energy irradiation decreasing with a shallow power-law up to 600 Myr; for older ages (decay phase), the steeper decrease of L_x is also described with slightly different slopes.

In order to evaluate the stellar irradiation in extreme ultraviolet band we adopted the scaling law derived by Sanz-Forcada et al. (2011), providing a relation between EUV and X-ray luminosity, calibrated on the stars observed in both bands. We took into account the evolution of the planetary radius through the relation given in Lopez & Fortney (2014), as already adopted, for instance, by Poppenhaeger et al. (2021). This relation was developed for H-He dominated atmospheres, and provides the envelope radius R_env, as a function of the planetary mass, the bolometric flux received, the atmospheric mass fraction f_atm, and the age of the system, allowing us to account also for the gravitational shrinking. Considering an upper limit for the mass of DS Tuc A b of 15 M_⊕, we performed several simulations for different values of the mass in the range 5–15 M_⊕. For each value, we estimated the core radius, R_core, through the following analytic relation, derived by fitting the data in Fig. 1 of Mocquet et al. (2014) for planets with Earth-like composition of the core: $$R_{\text{core}}=R_{\oplus }{\left(\frac{M_{\text{P}}}{M_{\oplus }}\right)}^{\frac{1}{3.7}}{e}^{[w(\frac{M_{\text{P}}-M_{\oplus }}{M_{\text{P}}})]},$$(9)

where M_P is the mass of the planet, and w = −0.014 is the fitting parameter. Given R_core and the measured planetary radius R_P = 5.6 R_⊕, we estimated the initial R_env as the difference between R_P and R_c. Inverting the relation of Lopez & Fortney (2014), we calculated the initial atmospheric mass fraction for each assumed planetary mass. For each time step of the simulation, we updated f_atm and the planetary mass in response to the mass loss, obtaining a new value of R_env. The latter quantity added to the core radius provided the updated planetary radius. For each mass of our grid we let the system to evolve from the current age, supposed to be 40 Myr, until 5 Gyr (Fig. 17, upper panel). We stopped the simulations when the cumulative percentage of the mass loss reached the corresponding atmospheric mass fraction, that is when the planet has lost completely its atmosphere.

Fig. 16 X-ray spectra of DS Tuc A (upper panel) and DS Tuc B (lower panel). Each panel is subdivided in two plots: spectrum with the best fit model (a sum of 3 APEC components, each one indicated by the dotted lines) and the terms of χ2 for each bin in units of χ2 standard deviation.

7.2 Results

In almost all cases the planet loses its atmosphere within 1 Gyr. The two X-ray evolution laws give results that are fully compatible when lower masses are considered. The most discrepant case is for a planet of 15 M_⊕. Following the P08a evolution, such a planet is subject to a complete atmospheric loss in 600 Myr, while J12 foresees that at such an age the planet retains still an atmosphere of ~ 1% of the planetary mass, which is completely dissipated at ~3 Gyr.

In the particular case in which the mass of the planet is ~ 8  M_⊕, or lower, the atmosphere is completely lost in less than 1 Myr, which is the time step of our simulations. In this case, if a planet atmosphere would be detected, we were in the fortunate position to witness the final stages in the life of a planetary atmosphere. In other words, if DS Tuc A b still has an atmosphere that we can detect, we can conclude that its mass is very likely greater than 8 M_⊕, given the known age of the system. However, the scenario in which the planet has a mass up to 8 M_⊕ and no atmosphere still remains plausible. These results can be compared with those presented by Owen (2020), which excluded a planetary mass lower than ~4.5 M_⊕. Instead, a mass up to 7 M_⊕ requires a boil-off phase during planet formation, while masses between 7 and 10 M_⊕ would be consistent with a standard picture of core accretion.

We also investigated the separate contribution of the X-ray irradiation with respect to the X+EUV irradiation. Figure 17 (upper panel) shows the previous simulations, in the case of the P08a evolution (dashed lines), compared to the same model evaluated considering the X-ray flux only (dotted lines). Clearly, the EUV contribution to the evaporation is dominant, since the atmosphere of planets with masses > 10 M_⊕ would be retained up to 5 Gyr in absence of it. In the case of a 8 M_⊕ planet, complete evaporation would occur in 200–300 Myr.

According to our models, at the end of the evaporation phase the planetary radius will be about 2.0, 1.9, 1.8, and 1.7 R_⊕, namely coincident with the core radius, in the simulations with initial mass of 15, 12, 9, and 8 M_⊕, respectively. Figure 17 (lower panel) shows the radius evolution foreseen with the P08a evolution. Similarly, the final planetary mass is evaluated as 13.6, 11.0, 8.2, and 7.3 M_⊕, respectively.

In all of the cases explored, at the end of the evaporation phase the planet would lose approximately ~ 10% of its total mass, whereas its radius would be reduced by a factor of ~3. Our models foresee that DS Tuc A b would conclude in any case its evolution within the gap of the Fulton radius distribution (Fulton et al. 2017), where approximately half of all planets have lost completely their atmospheres (Modirrousta-Galian et al. 2020). The final density will be about 20 times higher than the initial one.

Fig. 17 Summary of our simulations on the atmospheric evolution of DS Tuc A b. Upper panel: cumulative percentage of the planetary mass loss evolution considering the X+EUV irradiation (dashed) or X-ray irradiation alone (dotted). We used the results obtained with the evolution by Penz et al. (2008). Horizontal solid lines represent the limits of the atmospheric mass concerning four representative cases with different starting masses, colour-coded as in the legend. Lower panel: radius evolution corresponding to the representative cases according to our models from the present age of the system to the end of the photo-evaporation phase.

8 Discussion

8.1 Young planet parameters in context: mass-radius

The ongoing characterisation of transiting exoplanets at young age is providing us the first clues of their population distribution,at least with respect to the mature systems. By using the available data, we aimed to update the emerging picture with a sample of well constrained adult systems. We selected those exoplanets having measured masses and radii with uncertainty better than 20 and 10%, respectively, by using the Exo-MerCat tool (Alei et al. 2020)⁴. We used these systems to plot a mass-radius diagram (grey dots) in Fig. 18. We positioned the young transiting planets with measured mass known to date as a function of their age (left panel) and planet insolation inEarth units (right panel). We also placed dashed lines indicating the traces with equal planet density. For planets younger than ~100 Myr we only have a mass upper limit, as in the case of DS Tuc A b, and for this reason, they are indicated in the plot with triangles pointing toward lower values of the planetary mass. Their actual location should be thus shifted more on the left. This condition is clearly related to the difficult RV modelling due to the contributionof the stellar activity. To built these plots, we considered those planets found by Kepler/K2 in young associations and the ones detected by TESS, besides DS Tuc A b. Among them, HIP67522 b (Sco-Cen, age ~ 17 Myr, Rizzuto et al. 2020), TOI-251 b (age ~18 Myr, Zhou et al. 2021), AU Mic b and c (age ~22 Myr, Plavchan et al. 2020; Martioli et al. 2021), TOI-837 b (age ~ 35 Myr, Bouma et al. 2020) and TOI-942 b and c (age ~50 Myr, Carleo et al. 2021; Zhou et al. 2021). From Fig. 18 (left panel), we notice that all the planets younger than 100 Myr show radii larger than ~4 R_⊕ and they are located in the upper portion of the distribution for intermediate age and mature planets with similar masses. DS Tuc A b is the low-mass planet with the lowest density in the diagram, except for the system of Kepler-51⁵. Currently, its density (<0.44 g cm⁻³) is similar to the Jupiter-like planet Kepler-289 d and, according to our simulations, its mass-radius evolution is expected to follow the black arrow down to the hatched rectangle. The expected density at this location ranges from 5.6 to 12.8 g cm⁻³ and according to the models by Zeng & Sasselov (2013), this should correspond to a composition between 100% MgSiO₃ and 50%Fe + 50% MgSiO₃. The most curious result is that the final radius of DS Tuc A b is expected to be within the Fulton gap.

The right panel of Fig. 18 shows that DS Tuc A b is the less dense planet in the insolation range between ~10 and ~500 S_⊕. Again, we must recall that some of those planets without a robust measured mass could be less dense as well. More in general, Fig. 18 shows that, except for a few outliers⁶, young planets seems to lie within or close to the distribution of relatively older planets, with DS Tuc A b located on the low-density edge of it.

Fig. 18 Planetary mass-radius diagram for a sample of characterised young transiting planets as a function of the stellar age (left) and planet insolation (right). Coloured circles indicate planet with measured mass, while triangles indicate planets withmass upper limit only. Grey dots represent a sample of well-characterised mature exoplanets for comparison. Dashed lines represent loci of equal density, as indicated at the top of the figure. Hatched rectangle indicates the locus of the final mass and radius of DS Tuc A b, according to our planetary evolution models.

Fig. 19 Planetary period-radius diagram for the sample of characterised young transiting planets as a function of the stellar age. Grey dots represent all the known planets with measured period and radius. Black bar line indicates the locus of the final radius of DS Tuc A b according to our planetary evolution models.

8.2 Young planet parameters in context: period-radius

Similarly to Fig. 18, we produced a period-radius diagram in Fig. 19. It shows a more clear separation between young and larger radii planets (red dots) with respect to intermediate-age and smaller radii planets (other colours), further noticing that almost all the young planets lay in the less populated regions of the plot. We must take into account that this distribution could be the result of observation biases. For instance, transiting planets with shorterradii can hardly be detectable when orbiting young stars affected by high photometric variability. This case includes a potential population of non-inflated young planets. As a consequence, these type of planets possibly have lower masses that are more difficult to measure with RVs, as clearly demonstrated in this work. To investigate this issue, a frequency and population synthesis study should be done, as the one presented by Mann et al. (2010) in the case of low-mass planets at wide separations. However, this topic is out of the purpose of the present study.

Another possible bias concerns the timescale of the low-mass planet formation if they follow the same formation path of the terrestrial planets in the Solar System. According to the Grand Tack theory, the latter are formed from the material left behind by the violent destruction of the first generation planets, due to the gravitational interaction produced by Jupiter (Batygin & Laughlin 2015). If the timescale of the second generation planet formation ranges between 30 and 100 Myr (Kleine et al. 2009), planets around younger stars are still not present in the observed distribution. Nevertheless, the diagram in Fig. 19 clearly shows the expected radius evolution toward lower values as the system ages, but it can also be used to evaluate possible orbital evolution that could modify the planet period. In the case of DS Tuc A b, we do not expect a further inward migration through the disc since no IR excess has been detected (Zuckerman et al. 2011) and, therefore, the disc is already dissipated. On a larger timescale we must consider dynamic interactions with additional bodies in the system able to interact with DS Tuc A b. From Benatti et al. (2019), we excluded planets with 2–3 Jupiter masses and period shorter than 10 days and 7 M_J within 100 days, from RV data, and 7–8 M_J at separation larger than 40 AU from direct imaging. DS Tuc B resides at 240 au apart. From the available information, we cannot confirm the presence of additional perturbing bodies in the system. Finally, star-planet tidal interactions can occur and produce modifications of the semi-major axis of the orbit. According to the models by Bolmont & Mathis (2016), tidal interactions between a solar-type star and a close-in planet with mass of 10 M_⊕ could shrink its orbit if the planet semi-major axis is ≲0.04 au after the disc dissipation. Given that the distance between DS Tuc A and its planetary companion is about 0.09 au, we do not expect a significant modification of the planetary orbit through this specific channel.

The expected location of DS Tuc A b in the diagram at the end of the photo-evaporation phase is therefore within the known distribution (indicated with a black bar line).

8.3 Challenging detection of the young planets’ atmosphere

As reported in Sect. 5, no evidence for atmospheric compounds in the atmosphere of DS Tuc A b is found in our high-resolution ESPRESSO data. The lack of Hα detection could be explained by a large fraction of ionised planetary escaping material produced by the interaction with the strong XUV stellar radiation (see e.g. Shaikhislamov et al. (2020) in the case of Lyα in π Men c). Alternatively, a suppression of the planetary mass-loss rate can be observed in presence of strong magnetic field both of the star and the planet (Owen & Adams 2014).

Recent studies seem to confirm that a flat spectrum is a characteristic of young planetary atmospheres, despite their large scale height. Libby-Roberts et al. (2020) and Thao et al. (2020) explained the featureless spectra of the intermediate-age planets Kepler-51, b and d, and K2-25 b, respectively, as due to the presence of photochemical hazes. This condition is predicted for Neptune-size companions by Crossfield & Kreidberg (2017), since the equilibrium temperature of those planets is less than 800 K.In the case of AU Mic b, Palle et al. (2020) indicated the stellar activity as the main responsible for the non-detection of atomic species. We must note that the latter case is more similar to the one we have analysed in this work. Because of its young age (about 20 Myr), a high level of stellar activity is present in AU Mic, possibly comparable to the one of DS Tuc A. Their equilibrium temperatures are 850 K (DS Tuc A b, Newton et al. 2019)and 600–800 K (AU Mic, Palle et al. 2020), which are close to the transition threshold set by Crossfield & Kreidberg (2017) from obscuring hazes to clear atmosphere. This fact complicates possible considerations unlike the case of Kepler-51 and K2-25. However, the possible presence of clouds or hazes with Rayleigh scattering masking the Hα detection could be further explored with low resolution transmission spectroscopy.

We can conclude that at very young age, the contribution of the activity from the star dominates the behaviour of planetary spectral features with respect to all the other components of the planetary atmosphere.

9 Conclusions

In this work, we performed a first attempt to characterise the young star-planet system of DS Tuc A by using different techniques. We monitored this target with HARPS and ESPRESSO aiming to search for a RV signal of the planetary companion and the detection of its atmospheric compounds. Despite the fact that the high level of the stellar activity in such a young host made our analysis very challenging, we were able to provide a robust mass upper limit of the planetary companion (14.4 M_⊕) that allowedus to confirm its low density. We confirm a certain degree of alignment of the system through the modelling of the RM effect. Stellar activity also impacted on the atmospheric detection, reinforcing the evidence that the atmosphere of young planets are not easily accessible. We also observed this target with XMM-Newton and we were able to estimate for the first time the flux, luminosity, and coronal plasma temperatures for each stellar component separately. The X-ray luminosity of DS Tuc A served as input for models of mass evolution of young planets in order to estimate mass evaporation as a function of age. From the photometric point of view, a deeper characterisation of DS Tuc A b and its host will be provided by the new available observations of the TESS cycle-2 (Sectors 27 and 28), and from a dedicated monitoring with CHEOPS (AO-1:ID-0020, PI: Desidera, Scandariato et al. in prep.). To overcome the effect of the stellar activity in the detection of atmospheric compounds, observation in the near-Infrared band is an option. Besides lowering the impact of the activity, this spectral range offers the possibility of recovering Helium (e.g. Allart et al. 2019; Guilluy et al. 2020) and several molecules. Spectrographs like CRIRES+ (Follert et al. 2014) at VLT and NIRPS (Bouchy et al. 2017) at ESO 3.6 m telescope could give an important contribution to this issue.

Acknowledgements

We thank the anonymous Referee for the useful comments and suggestions that improved the quality and the clarity of the paper. We thank F. Bouchy and X. Dumusque for managing the sharing of observing time between several programmes approved for HARPS in P103 and the observers of these programmes for carrying out observations of our target in their nights. We thank C. Lovis for the support in the reprocessing of the HARPS data. The authors acknowledge the support of the ARIEL ASI-INAF agreement n.2018-22-HH.0. We acknowledge financial support from the ASI-INAF agreement n.2018-16-HH.0. We acknowledge the support by INAF/Frontiera through the ‘Progetti Premiali’ funding scheme of the Italian Ministry of Education, University, and Research.

References

All Tables

All Figures

	Fig. 1 HARPS RV time series of DS Tuc A as extracted from the TERRA pipeline. A bulk of data at JD = 8706 corresponds to the RVs collected during the planetary transit.
In the text

	Fig. 2 RM effect for the transits observed with HARPS (red dots, TERRA; green dots, DRS, systemic RV removed) and ESPRESSO (blue dots, DRS, systemic RV removed). Vertical lines indicate the first and fourth contact (dotted) and the centre of the transit (dashed), according to the ephemeris in Newton et al. (2019).
In the text

	Fig. 3 GLS of the HARPS RV (black line) and BIS (red line) time series of DS Tuc A. The green vertical line indicates the locationof the planet orbital period, while the window function of the observations is reported in the inset plot. Dotted, dashed, and solid black horizontal lines indicate the FAP thresholds at 1, 0.1, and 0.01%, respectively, for the RV periodogram.
In the text

	Fig. 4 Corner plots for GP quasi-periodic (upper panel) and SHO mixture (lower panel) models used to fit the HARPS+ESPRESSO radial velocities.
In the text

	Fig. 5 Spectroscopic orbit of planet DS Tuc A b. The instrumental offsets and the best-fit stellar activity signal (GP quasi-periodic model) have been subtracted from the HARPS and ESPRESSO original RVs. The error bars include the uncorrelated jitter terms added in quadrature to the formal uncertainties. The curve in red represents the best-fit Doppler signal (circular obit; Kb = 3.2 m s−1).
In the text

Fig. 6 Summary of the GP best-fit model as obtained from the joint analysis of the HARPS+ESPRESSO RVs. Here, we show a portion of the HARPS time series, after removing the best-fit signal induced by planet b. The shaded area in gray represents the ±1σ confidence interval. Upper panel: best-fit model for the GP quasi-periodic kernel (orange line). Lower panel: same as above, but for the GP rotation kernel (blue line).

In the text

	Fig. 7 Posteriors of the Kb parameterretrieved from the simulations involving the original RV datasets, with the injected Kb, inj values indicated in the legend. The vertical dashed lines represent the median of the posteriors of the corresponding colour. The vertical black lines indicates the median of the posteriors for the real data (3.2 and 8.6 m s−1), as a reference.
In the text

Fig. 8 Posteriors of the Kb parameter retrieved from the simulations involving mock RV datasets build on the best-fit stellar activity signals based on a quasi-periodic and SHO mixture models (see the text for details). The simulated RV time series include an injected planetary signal of semi-amplitude Kb, inj = 40 m s−1. Only five out of a total of ten posteriors are shown for the sake of a better the readability First row: results from the fit of data generated using a quasi-periodic kernel. Second row: results from the fit of data generated using a SHO mixture kernel. In both rows, each colour identifies one mock dataset, therefore, the posterior distributions with the same colour can be compared. The vertical lines indicate the median of the posteriors.

In the text

	Fig. 9 Radial velocities and activity indices variation during the transit of Oct. 6 of DS Tuc A b observed with ESPRESSO. Upper panel: RV (blue dots) with the GP model of the stellar activity (orange line) and the corresponding residuals. Lower panel: CaII H&K and Hα time series extracted from the same spectra.
In the text

	Fig. 10 RM effect for DS Tuc A b from ESPRESSO data modelled by our fit (red line) with the corresponding residuals. Vertical dotted lines indicate the transit ingress and egress.
In the text

	Fig. 11 CCF residuals during the transit of DS Tuc A b in the stellar restframe. Horizontal white lines indicate the transit ingress and egress.
In the text

	Fig. 12 2D tomography of Hα in the stellar restframe. The vertical dotted line shows the stellar restframe. The diagonal magenta dashed line shows the expected planetary restframe. The horizontal white lines mark the beginning and end of the transit.
In the text

	Fig. 13 Transmission spectrum extracted in the zone of the Hβ (left panel), Na D2 (central left panel), Hα (central right panel), and Li (right panel) lines. The vertical dashed lines show the planetary restframe. Black points refer to 0.1 Angstrom binning. Different lines we searched for are affected by the contamination given by the stellar activity in different ways.
In the text

Fig. 14 XMM/Newton MOS1 image of DS Tuc. The two red points indicate the ICRS positions of the two stars from SIMBAD, the centroid of X-ray emission is offset because of the proper motions of the system. The contours show the emission of the two sources, marginally resolved. The large circle (in red) has 30′′ radius, while the smaller circles (in yellow) are 10′′ wide. They were used to extract the spectra and light curves of the A and B system components (see text).

In the text

	Fig. 15 Mos1 light curves for both components of DS Tuc. Black line is the light curve from the semi-circular region of radius 30′′, red line is the light curve from the circular region of radius 10′′. Error bars on rates are at 68% confidence level. Upper panel: DS Tuc A. Lower panel: DS Tuc B.
In the text

	Fig. 16 X-ray spectra of DS Tuc A (upper panel) and DS Tuc B (lower panel). Each panel is subdivided in two plots: spectrum with the best fit model (a sum of 3 APEC components, each one indicated by the dotted lines) and the terms of χ2 for each bin in units of χ2 standard deviation.
In the text

Fig. 17 Summary of our simulations on the atmospheric evolution of DS Tuc A b. Upper panel: cumulative percentage of the planetary mass loss evolution considering the X+EUV irradiation (dashed) or X-ray irradiation alone (dotted). We used the results obtained with the evolution by Penz et al. (2008). Horizontal solid lines represent the limits of the atmospheric mass concerning four representative cases with different starting masses, colour-coded as in the legend. Lower panel: radius evolution corresponding to the representative cases according to our models from the present age of the system to the end of the photo-evaporation phase.

In the text

Fig. 18 Planetary mass-radius diagram for a sample of characterised young transiting planets as a function of the stellar age (left) and planet insolation (right). Coloured circles indicate planet with measured mass, while triangles indicate planets withmass upper limit only. Grey dots represent a sample of well-characterised mature exoplanets for comparison. Dashed lines represent loci of equal density, as indicated at the top of the figure. Hatched rectangle indicates the locus of the final mass and radius of DS Tuc A b, according to our planetary evolution models.

In the text

	Fig. 19 Planetary period-radius diagram for the sample of characterised young transiting planets as a function of the stellar age. Grey dots represent all the known planets with measured period and radius. Black bar line indicates the locus of the final radius of DS Tuc A b according to our planetary evolution models.
In the text