A hot Jupiter orbiting a 2-Myr-old solar-mass T Tauri star

Hot Jupiters are giant Jupiter-like exoplanets that orbit 100x closer to their host stars than Jupiter does to the Sun. These planets presumably form in the outer part of the primordial disc from which both the central star and surrounding planets are born, then migrate inwards and yet avoid falling into their host star. It is however unclear whether this occurs early in the lives of hot Jupiters, when still embedded within protoplanetary discs, or later, once multiple planets are formed and interact. Although numerous hot Jupiters were detected around mature Sun-like stars, their existence has not yet been firmly demonstrated for young stars, whose magnetic activity is so intense that it overshadows the radial velocity signal that close-in giant planets can induce. Here we show that hot Jupiters around young stars can be revealed from extended sets of high-resolution spectra. Once filtered-out from the activity, radial velocities of V830 Tau derived from new data collected in late 2015 exhibit a sine wave of period 4.93 d and semi-amplitude 75 m/ s, detected with a false alarm probability<0.03%. We find that this signal is fully unrelated to the 2.741-d rotation period of V830 Tau and we attribute it to the presence of a 0.77 Jupiter mass planet orbiting at a distance of 0.057 au from the host star. Our result demonstrates that hot Jupiters can migrate inwards in<2 Myr, most likely as a result of planet-disc interactions, and thus yields strong support to the theory of giant planet migration in gaseous protoplanetary discs.

Very few exoplanets have yet been discovered around young forming Sun-like stars aged <10 Myr 7,8 , the so-called T Tauri stars (TTSs), either through the radial velocity (RV) variations or the photometric transits they induce in the light of their stars. Yet detections of young planets are key for our understanding of how planetary systems form; this is especially true of young hot Jupiters (hJs) thought to have a critical impact on the early architecture of these systems. The first claimed detection of a young hJ 4 orbiting a TTS was quickly refuted; the reported periodic RV fluctuations were finally attributed to activity 5 and to cool spots at the stellar surface 9 . The recent candidate detection of a transiting hJ around a TTS 6 is still pending confirmation.
TTSs are known to harbor spots and plages at their surfaces, generating RV fluctuations with semi-amplitudes of several km/s 10 , i.e., much larger than the perturbations expected from a putative planet, even for close-in massive hJs inducing typical RV signals of ~0.1 km/s. Detecting hJs around TTSs through velocimetry or photometry is thus quite challenging and mandatorily requires efficient tools for filtering-out the dominant jitter that activity induces in the spectra and light-curves of young stars. We recently proposed a new method to achieve this goal 11 , whose first applications to wTTSs proved promising though inconclusive 12 .
V830 Tau is a ~2-Myr-old solar-mass TTS 12 contracting towards the main sequence and currently spinning in 2.741 d 13 , i.e., ~10x faster than the Sun. Evolutionary models 14 suggest that it is fully or largely convective. Unlike 80% of the TTSs in the Taurus star-forming region 15 , V830 Tau exhibits no significant infrared excess, implying that most of its inner accretion disc already dissipated. This is consistent with its status as a non-accreting weak-line TTS (wTTS), and makes it an ideal target to look for the presence of hJs at an early stage of star and planet formation.
In late 2015, we collected 48 high-resolution spectra of V830 Tau (see Extended Table 1) as part of the MaTYSSE Large Program aimed at detecting hJs around wTTSs 11 . Applying leastsquares deconvolution 16 (LSD) to our spectra, we derived accurate average line profiles and their temporal modulation over ~15 rotation cycles. Longitudinal magnetic fields were also derived from our circularly-polarized data and the Zeeman signatures that fields generate in spectral lines 16 . Using tomographic techniques inspired from medical imaging, one can reconstruct distributions of spots and plages at the surfaces of rotating cool active stars from sets of densely-sampled line profiles covering several rotation cycles. This method, called Doppler imaging 17 (DI), can also probe the photospheric shear associated with surface differential rotation through the amount of twisting it generates in brightness maps 18,19 .
Our DI code was previously applied to a small set of 15 LSD profiles of V830 Tau, from which the distribution of surface features and the differential rotation pattern were recovered 12 ; it even suggested the potential presence of a hJ, though with a very low confidence level. The version presented herein implements a novel technique to filter-out lunar contamination (plaguing spectra collected in non-photometric conditions) yielding good results when phase coverage is dense. Applying the code to our new set of 48 LSD profiles of V830 Tau yields the map and fit shown in Fig 1. We again clearly detect differential rotation (see Extended Fig 1) and confirm that it is ~3x weaker than that of the Sun, reflecting that V830 Tau is largely or fully convective 20 .
From the brightness image reconstructed with DI, we derive the model RV curve that V830 Tau should exhibit if all profile perturbations were attributable to surface features and differential rotation. By subtracting these modeled RVs from the observed ones (both computed as the first moment of the LSD profiles), we obtain the activity-filtered RVs of V830 Tau whose amplitude is typically ~10x lower than the raw RVs (see Fig 2). A clear RV signal, with a period of 4.93±0.05 d and a semi-amplitude of 75±11 m/s, is detected in the activity-filtered RVs at a confidence level >99.97% (see Figs 2 & 3). Using a Bayesian approach 21 , we find that the 4.93-d peak is at least 10 5 times more likely than any of the other features in the periodogram (see Extended Fig 2a). The regular phase coverage of our data also allows us to check that the 4.93-d signal is present in smaller subsets (e.g., first and second half, even and odd points) though expectedly with a lower confidence level; we similarly checked that our detection holds when profiles affected by lunar contamination are excluded (see Fig 3 and Extended Fig 2a) and when differential rotation is neglected. Periodograms of the longitudinal fields and of the Hα emission, both reliable proxies for the activity jitter plaguing RV curves 22 , show no power at a period of 4.93 d (see Fig 3b & Extended Fig 2b), demonstrating that the signal we report is unrelated to activity.
We interpret this RV signal as caused by a giant planet of mass 0.77±0.15 M♃ in circular orbit around V830 Tau at a distance of 0.057±0.001 au (see Extended Table 2). Although the filtered RVs are marginally better fitted for an eccentric orbit (e=0.30±0.15, see Fig 2), we still favor a circular orbit given the large error bar on the eccentricity 23 . Removing the planet signal from the original data and repeating the activity-filtering process yields residual RVs with a rms dispersion of 48 m/s (i.e., consistent with the average noise level, see Extended Table 1) and no significant peak left in the periodogram. Using the alternative option of fitting both the brightness map and the planet parameters simultaneously 24 yields identical results and demonstrates that our optimal model including a planet is orders of magnitude more likely than a model with no planet (see Fig 4). Simulations in conditions identical to those of our observations yield results in close agreement with those of Figs 2 and 3 (see Extended Figs 3 & 4), further demonstrating that our filtering process induces no spurious peak and that the RV signal we detect is unrelated to activity.
A careful re-analysis of our original data 12 demonstrates that, despite being affected by intrinsic variability from the host star, they nonetheless confirm the presence of the planet signal detected in our new data. In particular, applying our filtering analysis to the main subset of our original data (featuring dense enough coverage for our technique to perform reliably) yields filtered RVs agreeing well with those derived from the new data; fitting both sets together further improves the confidence level at which the planet is detected (see Extended Fig 2c).
The detection we report among the small sample of wTTSs (~10) already studied with MaTYSSE suggests that close-in giant planets are potentially more frequent around TTSs than around mature low-mass stars, ~1% of which are known to host hJs 25,26 . Poisson statistics however indicate that there is still a ~10% chance that hJs are similarly frequent for both populations. A more quantitative conclusion will have to await for a thorough analysis of all MaTYSSE data, later-on complemented by the large TTS survey to be carried out with SPIRou, the new nIR cryogenic spectropolarimeter and velocimeter to be installed at CFHT in 2017.
The broad consensus is that hJs form beyond a few au from their host stars and migrate inwards to their close-in orbits. The delivery of hJs may result from planet-disc interactions 2 (disc migration) or from dynamical interactions with a planetary 3 or a stellar companion, followed by orbital circularization (high-eccentricity migration). These highly-debated migration channels were proposed as tentative explanations for hJs with orbits either well aligned or misaligned with the spin axis of their host stars 27 . Our detection of a hJ on a ~5-d circular (or moderately eccentric) orbit around a 2-Myr-old star is most naturally explained with hJ delivery by disc migration, producing hJs with nearly circular orbits, rather than by planet-planet scattering, generating highlyeccentric (e>0.9) hJs whose circularization timescales are at least 2-3 orders of magnitude longer 28 than the age of V830 Tau (for typical tidal dissipation factors in giant planets 29 ). Our result thus yields strong support to the theory of Type-II migration of giant planets in gaseous protoplanetary discs 2 and confirms that the architecture of planetary systems, on which hJs have a strong impact, is likely quite dynamic right from the very early ages of planetary formation.
Global models of planet formation and evolution show that giant planets can reach a mass and an orbital period similar to those of V830 Tau b in 2-3 Myr, whether formation occurs through core accretion or disk gravitational instability. Large uncertainties in these models make it difficult to accurately predict the occurrence rate of hJs and thus to associate V830 Tau b to either formation scenario. The ~300 G dipole of the star's magnetic field 12 is strong enough to have disrupted the central 0.06 au of the now-dissipated disc for accretion rates <2.10 -10 M⊙/yr; for rates of ~10 -9 M⊙/yr, more typical to those of the classical TTSs (cTTSs) that still feed from their discs, a >700 G dipole field is required, again compatible with observations of cTTSs similar to V830 Tau 30,12 . This shows that the field of V830 Tau may well have stopped the planet within the magnetospheric gap 1 at the end of its disc migration, and saved it from falling into the star.
Fédérale Toulouse Midi-Pyrénées for awarding a "Chaire d'Attractivité" to GH, in the framework of which this work was done. SA acknowledges financial support from CNPq, CAPES and Fapemig.
Author Contributions This work merged data collected with 2 different instruments and 3 different telescopes. JFD led the data processing, analysis and manuscript preparation, CM, LM, LY and EH participated in the data collection and data analysis, whereas CB contributed to the theoretical implications of the results. All coauthors, including GH, SA, MT, FM, JB, PP, RD and ACC, were involved in elaborating the observing proposals, in discussing the results at various stages of the analysis, and in providing contributions to older versions of the manuscript.

Author Information
The authors declare no competing financial interests. Correspondence should be addressed to JFD (jean-francois.donati@irap.omp.eu)   Tau (open  symbols and 1σ error bars) and model inferred with Doppler imaging (cyan line); the model slowly evolves with time as a result of differential rotation. Open circles, squares and triangles depict ESPaDOnS, NARVAL and ESPaDOnS / GRACES data. Middle curve: activity-filtered RVs and best sine fit to the data (cyan line). The period and semi-amplitude of the planet RV signal are equal to 4.93±0.05 d and 75±11 m/s (1σ error bars). Bottom curve: residual RVs once the planet signal is removed and activity is filtered out, with a final rms dispersion of 48 m/s. Colors code the rotation cycles. b, Activity-filtered RVs phase-folded on the planet orbital period of 4.93 d. Although the fit to the data is marginally better with an eccentric orbit (dashed line) than with a circular orbit (solid line), the significance of the derived eccentricity (0.30±0.15) is too low to be reliable 23 .  Variations in the reduced χ 2 of the DI fit to the LSD profiles for a given spottedness level and assuming the presence of a planet in circular orbit, for a range of orbital periods Porb and semi-amplitudes K of the RV signature. (This is a 2D cut from a 3D map, with the phase of the RV signal also included as a search parameter). The location of the minimum and local paraboloid curvature yield the optimal values of Porb and K and their respective 1σ error bars 24 , equal to 4.94±0.05 d and 82±10 m/s, in agreement with the results of our main filtering technique. The outer color contour traces the projected 99.99% confidence interval, corresponding to a χ 2 increase of 21.1 for a 3-parameter fit to the 2208 data points of the LSD profiles. With respect to our best model incorporating a planet, a model with no planet corresponds to a χ 2 increase of 82, implying that the planet is detected with a FAP <10 -15 , much smaller than the FAP derived from the periodogram of the 48 RV points (see Fig 3) thanks to the larger number of data points in the fitted LSD profiles.

Spectropolarimetry with ESPaDOnS / NARVAL and the MaTYSSE programme
ESPaDOnS 31 and NARVAL are twin spectropolarimeters respectively installed at the Cassegrain focus of the 3.6m Canada-France-Hawaii Telescope (CFHT) atop Maunakea (Hawaii), and of the 2m Bernard Lyot Telescope atop Pic du Midi (France). Both include a fiber-fed bench-mounted high-resolution spectrograph, yielding full coverage of the 370-1,000 nm wavelength range in a single exposure at a spectral resolving power of 65,000. ESPaDOnS can also be fed from the 8m Gemini-N Telescope next to CFHT, through a 300m fiber link called GRACES 32 , yielding spectra with either similar resolving power (in star-only mode) or half of it (in star+sky mode). In this run, we secured 16 spectra with ESPaDOnS from 2015 Nov 17 to Dec 02 in fair weather conditions, 16 spectra with NARVAL from 2015 Nov 10 to Dec 17 in moderate to good weather, and 16 spectra with ESPaDOnS / GRACES on 4 different nights towards the end of the run (with 2 / 14 spectra in star-only / star+sky mode respectively). The full journal of observations is given as Extended Table 1 All spectra were derived from raw frames with the reference pipeline implementing optimal extraction and RV correction from telluric lines 16 , yielding a typical RV precision of 30 m/s rms 33 .

Deriving mean line profiles with least-squares deconvolution (LSD)
LSD 16 is a multiline technique similar to cross-correlation, used to derive line profiles with enhanced S/N from thousands of spectral lines simultaneously. For this study, the line list we used for LSD is derived from spectrum synthesis through model atmospheres computed assuming local thermodynamic equilibrium 34 , for atmospheric parameters relevant for V830 Tau 12 (effective temperature of 4250 K, logarithmic gravity of 4.0 and solar metallicity). Resulting S/N in LSD profiles are in the range 950-1540 (see journal of observations, Extended Table 1), corresponding to average multiplex gains in S/N of ~10. LSD was also applied to circularly polarized spectra to retrieve average Zeeman signatures and longitudinal field estimates 16 .

Doppler imaging (DI) of stellar surfaces and the modeling of differential rotation
DI is a tomographic technique inspired from medical imaging, with which distributions of brightness features and magnetic fields at the surfaces of rotating stars can be reconstructed from time-series of high-resolution spectropolarimetric observations. DI is based on the fact that, thanks to the Doppler effect, line profiles of rotating stars can be interpreted as one-dimensional (1D) images of stellar surfaces, resolved in the Doppler direction but otherwise blurred. By coupling many such 1D images recorded at different rotation phases, one can reliably reconstruct the parent surface distribution that gives rise to the observed line profiles and rotational modulation. First introduced in the late 1980s 17 , DI has been extensively used to investigate with unprecedented accuracy surface features and magnetic fields in cool stars other than the Sun 35,36 . Technically speaking, DI follows the principles of maximum-entropy image reconstruction, and iteratively looks for the image with lowest information content that fits the data at a given χ 2 level. For more details on the imaging process, our previous MaTYSSE studies 11,12 give a detailed account of all modeling steps. In the case of the present paper, we carried out DI modeling using either unpolarized (Stokes I) spectra only, or both unpolarized and circularly polarized (Stokes V) spectra simultaneously, with identical results regarding filtering performances and RV results. By looking at how surface maps get twisted as a function of time, DI can also estimate the amount of latitudinal differential rotation shearing stellar photospheres 18,19 . In this study we assume a typical solar-like differential rotation law in which the surface rotation rate varies with latitude θ as sin 2 θ, and depends on 2 main parameters, the rotation rate at the equator Ωeq and the difference in rotation rate dΩ between the equator and the pole. Both parameters are derived by looking for the pair that minimizes the χ 2 of the fit to the data (at constant information content in the reconstructed image, see Extended Fig 1), whereas the corresponding error bars are computed from the curvature of the χ 2 paraboloid at its minimum 19 . Although helpful to achieve a more accurate description of the activity jitter and a cleaner filtering of raw RV curves (at periods Prot and Prot/2 in particular), differential rotation as weak as that of V830 Tau has little impact on the filtered RVs; similar conclusions regarding the planet signal are obtained when assuming that V830 Tau is rotating as a solid body. A similar DI-based technique can be used to diagnose the presence of hJs around active stars 24 , with the planet parameters replacing those describing differential rotation. This alternate method yields identical results to those presented here for our data set (see Fig 4).

Filtering LSD profiles from lunar contamination
Spectra recorded in non-photometric conditions near full-moon epochs are often contaminated by solar light reflected off the moon and diffused by clouds. To filter-out this pollution, whose location and width is well known at any given epoch but whose strength we want to determine, we implement a dual-step DI process. The first step consists in applying DI to the set of original LSD profiles, with scaled-up error bars for all pixels potentially affected by lunar contamination; the strength of the lunar contamination is then measured with a gaussian fit to the residuals, and subtracted from the polluted LSD profiles. In the second step, conventional DI is applied to the set of filtered LSD profiles with original error bars. This dual-step process is found to be very efficient when applied to densely-sampled data sets like the one presented here, in which profiles at similar phases but different rotation cycles provide a strong constraint on the strength of lunar pollution. In our data, 6 LSD profiles suffer from a strong pollution (rotation cycles 6.0 to 7.6, see Fig 1), whereas 9 others are affected at a much weaker level. If excluding the 6 strongly (resp all 15) moon-polluted profiles from our data set, the RV signal from V830 Tau b is still clearly detected, though with a lower confidence rate of 99.9% (resp 99%) reflecting the poorer temporal coverage and the degraded window function (see Fig 3a and Extended Fig 2a). This check shows that our decontamination process is successful at restoring the original profile distortions and at retaining the RV content, provided the data set is dense enough.

Revisiting the original data set from 2014 Dec & 2015 Jan
A careful re-analysis of our original data (consisting of two subsets shifted in time by 17 d 12 , see Extended Table 3) indicates that variability occurred at the surface of the star between the two subsets. While DI succeeds at adjusting the main subset (9 evenly-spaced points secured in 2015 Jan) down to noise level, as for our new data, fitting both subsets together requires to lower S/Ns in the 2014 Dec subset by ~15% in order to reach unit reduced χ 2 . This is definite evidence that intrinsic variability (beyond pure differential rotation) occurred at the surface of V830 Tau between both subsets, with profile #6 of the 2014 Dec subset being the most affected. This variability reflects a modification in the brightness map, subtle enough to affect DI no more than moderately, yet large enough to significantly affect filtered RVs, which are quite sensitive to even small features in the brightness map. It illustrates how tricky activity filtering can get when dealing with intrinsic variability, and how critical dense and even phase coverage is to efficiently diagnose it. As a result of this variability, our filtering analysis can only be applied to the individual subsets of our original data, and in fact to no more than the 2015 Jan subset, the other being far too sparse and uneven for the technique to perform reliably. We find that the filtered RVs from the main subset (see Extended Table 3) agree with those of our new data; fitting them together improves the confidence level at which the planet is detected, but not the accuracy on the orbital period (see Extended Fig 2c). Enabling to shortcut the computation of filtered RVs, the DI-based method of adjusting the planet parameters simultaneously with the distribution of surface features 24 offers an alternative option for confirming that the RV signal from the detected planet is present in our original data. (This method however still suffers from DI's inability to describe intrinsic variability beyond differential rotation). Freezing the planet orbital period to the value found in our new analysis (4.93 d) and applying this technique to the full set of our original data with only profile #6 removed, we find a semi-amplitude of 67±18 m/s for the planet RV signature, which agrees well with the measurement derived from our main study. We also confirm that our original data are better explained by a model including a planet in a 4.93-d orbit than by one with no planet, with a false-alarm probability of ~0.01% (corresponding to a χ 2 increase of 19 for the 644 data points of the fitted LSD profiles).

Code availability
The Doppler imaging code used for this study is as yet undocumented and has thus not been released in the public domain.   Table 1 for the 2014 Dec and 2015 Jan data 12 . Our filtering technique is found to be only reliable for data sets with dense and regular phase coverage, hence the absence of filtered RVs for the 2014 Dec subset of 2x3 points, that does not satisfy this requirement. Filtered RVs from 2015 Jan 07-15 agree well with those derived from our new data (see Extended Data Table 1). Figure 1 | Estimating surface differential rotation parameter. Variations of the reduced χ 2 as a function of the surface differential rotation parameters Ωeq and dΩ, denoting respectively the rotation rate at the equator and the difference in rotation rate between the equator and the pole (and assuming a solar-like sine-square differential rotation law). The location of the minimum and the local paraboloid curvature yield the optimal parameters and their respective 1σ error bars 19  for simulated data computed using the brightness map, differential rotation and planet parameters inferred from the real data, and assuming the same coverage and similar S/N (equal for all LSD profiles). As for our observations, the planet signal is detected at a confidence level >99.9% in the filtered RVs despite being invisible in the raw RVs, and the planet parameters are well recovered. The periodogram of the raw RVs is very similar to that of Fig 3, featuring the main peaks (at Prot and Prot/2) and their aliases; residual RVs mostly reflect the noise in the data. b, Same as Extended Data Figure 3a but with no planet included in the simulation. No signal with a confidence level >90% is recovered in the filtered RVs, demonstrating that the filtering process is not generating spurious RV signals, in particular at a period of 4.93 d.