# A hot Jupiter around the very active weak-line T Tauri star TAP 26

L. Yu, J. -F. Donati, E. Hébrard, C. Moutou, L Malo, K Grankin, G. Hussain, A. C. Cameron, A. A. Vidotto, C. Baruteau, et al.

► To cite this version:

L. Yu, J. -F. Donati, E. Hébrard, C. Moutou, L Malo, et al. A hot Jupiter around the very active weak-line T Tauri star TAP 26. Monthly Notices of the Royal Astronomical Society, 2017, 467 (2), pp.1342-1359. 10.1093/mnras/stx009. hal-01769730

HAL Id: hal-01769730

https://hal.science/hal-01769730

Submitted on 28 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A hot Jupiter around the very active weak-line T Tauri star TAP 26

L. Yu<sup>1,2*</sup>, J.-F. Donati<sup>1,2</sup>, E. M. Hébrard<sup>3</sup>, C. Moutou<sup>4</sup>, L. Malo<sup>5</sup>, K. Grankin<sup>6</sup>, G. Hussain<sup>7,1</sup>, A. Collier Cameron<sup>8</sup>, A. A. Vidotto<sup>9</sup>, C. Baruteau<sup>1,2</sup>, S.H.P. Alencar<sup>10</sup>, J. Bouvier<sup>11,12</sup>, P. Petit<sup>1,2</sup>, M. Takami<sup>13</sup>, G. Herczeg<sup>14</sup>, S. G. Gregory<sup>8</sup>, M. Jardine<sup>8</sup>, J. Morin<sup>15</sup>, F. Ménard<sup>11,12</sup> and the MaTYSSE collaboration

*All affiliations are listed at the end of the paper*

Accepted XXX. Received YYY; in original form ZZZ

## ABSTRACT

We report the results of an extended spectropolarimetric and photometric monitoring of the weak-line T Tauri star TAP 26, carried out within the MaTYSSE programme with the ESPaDOnS spectropolarimeter at the 3.6 m Canada-France-Hawaii Telescope. Applying Zeeman-Doppler Imaging to our observations, concentrating in 2015 November and 2016 January and spanning 72 d in total, 16 d in 2015 November and 13 d in 2016 January, we reconstruct surface brightness and magnetic field maps for both epochs and demonstrate that both distributions exhibit temporal evolution not explained by differential rotation alone. We report the detection of a hot Jupiter (hJ) around TAP 26 using three different methods, two using Zeeman-Doppler Imaging (ZDI) and one Gaussian-Process Regression (GPR), with a false-alarm probability smaller than 6 $10^{-4}$. However, as a result of the aliasing related to the observing window, the orbital period cannot be uniquely determined; the orbital period with highest likelihood is 10.79±0.14 d followed by 8.99±0.09 d. Assuming the most likely period, and that the planet orbits in the stellar equatorial plane, we obtain that the planet has a minimum mass *M* sin*i* of 1.66±0.31 M<sub>Jup</sub> and orbits at 0.0968±0.0032 au from its host star. This new detection suggests that disc type II migration is efficient at generating newborn hJs, and that hJs may be more frequent around young T Tauri stars than around mature stars (or that the MaTYSSE sample is biased towards hJ-hosting stars).

**Key words:** magnetic fields – planets and satellites: formation – stars: imaging – stars: rotation – stars: individual: TAP 26 – techniques: polarimetric

## 1 INTRODUCTION

Studying young forming stars stands as our best chance to progress in our understanding of the formation and early evolution of planetary systems. For instance, detecting hot Jupiters (hJs) around young stars (1-10 Myrs) and determining their orbital properties can enable us to clarify how they form and migrate, and to better characterise the physical processes (e.g. planet-disc interaction, planet-planet scattering, Baruteau et al. 2014, in-situ formation, Batygin et al. 2016) responsible for generating such planets.

However, young stars are enormously active, rendering planet signatures in their spectra and / or light-curves ex-

tremely difficult to detect in practice. Until very recently, most planets found so far around stars younger than 20 Myr were distant planets detected with imaging techniques (e.g. β Pic b, Lagrange et al. 2010, and LkCa 15, Sallum et al. 2015). Early claims of hJs orbiting around T Tauri stars (e.g. TW Hya, Setiawan et al. 2008) finally proved to be activity signatures mistakenly interpreted as radial velocity (RV) signals from close-in giant planets (Huélamo et al. 2008).

Following the recent discovery of newborn close-in giant planets (Mann et al. 2016; David et al. 2016; Donati et al. 2016) or planet candidates (van Eyken et al. 2012; Johns-Krull et al. 2016) around forming stars, time is ripe for a systematic exploration of hJs around T Tauri stars, and in particular the so called weak-line T Tauri stars (wTTSs), whose accretion disc has just dissipated. This is one of the

<footnote>* E-mail: louise.yu@irap.omp.eu</footnote>
<table><caption>**Table 1.** Journal of ESPaDOnS observations of TAP 26 collected in 2015 November (first 16 lines) and 2016 January (last 13 lines). Each observation consists of a sequence of 4 subexposures, each lasting 695 s. Columns 1 to 4 respectively list (i) the UT date of the observation, (ii) the corresponding UT time (at mid-exposure), (iii) the Barycentric Julian Date (BJD) in excess of 2,457,300, and (iv) the peak signal to noise ratio (per 2.6 km s⁻¹ velocity bin) of each observation. Column 5 lists the root-mean-square (rms) noise level (relative to the unpolarised continuum level *I*<sub>c</sub> and per 1.8 km s⁻¹ velocity bin) in the circular polarisation profiles produced by Least-Squares Deconvolution (LSD) and column 6 lists the signal to noise ratio in the unpolarised profiles produced by LSD, measured from the noise level in intervals of continuum of the LSD profiles. Column 7 indicates the rotational cycle associated with each exposure (using the ephemeris given by Eq. 1). Column 8 lists the raw RVs computed from the unpolarised spectra, column 9 the filtered RVs (see Sec. 5.1) and column 10 the 1σ error bar on both RV<sub>raw</sub> and RV<sub>filt</sub>. Column 11-13 list values for activity proxies mentioned in App. B: the line-of-sight-projected magnetic field averaged over the visible stellar hemisphere (also called longitudinal field) and the equivalent width of the Hα emission (counted from above the continuum level, expressed in km s⁻¹, and with a typical 1σ error bar of 3.0 km s⁻¹).</caption><thead><tr><th>Date</th><th>UT<br>(h:m:s)</th><th>BJD<br>(2,457,300+)</th><th>S/N</th><th>σ<sub>LSD</sub><br>(10<sup>-4</sup>)</th><th>S/N<sub>I</sub></th><th>Cycle</th><th>RV<sub>raw</sub><br>(km s<sup>-1</sup>)</th><th>RV<sub>filt</sub><br>(km s<sup>-1</sup>)</th><th>σ<sub>RV</sub><br>(km s<sup>-1</sup>)</th><th>B<sub>t</sub><br>(G)</th><th>σ<sub>B<sub>t</sub></sub><br>(G)</th><th>EW<sub>Hα</sub><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>18 Nov</td><td>09:36:28</td><td>44.90594</td><td>140</td><td>3.3</td><td>1867</td><td>0.148</td><td>1.049</td><td>0.141</td><td>0.075</td><td>99</td><td>45</td><td>39.3</td></tr><tr><td>22 Nov</td><td>12:11:18</td><td>49.01352</td><td>140</td><td>3.3</td><td>1835</td><td>5.905</td><td>-1.115</td><td>0.026</td><td>0.076</td><td>-72</td><td>47</td><td>37.6</td></tr><tr><td>23 Nov</td><td>11:20:34</td><td>49.97830</td><td>140</td><td>3.1</td><td>1862</td><td>7.258</td><td>0.677</td><td>-0.120</td><td>0.075</td><td>-20</td><td>46</td><td>36.2</td></tr><tr><td>24 Nov</td><td>11:20:25</td><td>50.97819</td><td>140</td><td>3.0</td><td>1890</td><td>8.659</td><td>0.915</td><td>-0.020</td><td>0.074</td><td>-143</td><td>45</td><td>43.1</td></tr><tr><td>25 Nov</td><td>07:41:04</td><td>51.82588</td><td>140</td><td>3.3</td><td>1804</td><td>9.847</td><td>-0.017</td><td>-0.149</td><td>0.078</td><td>-182</td><td>47</td><td>44.0</td></tr><tr><td>25 Nov</td><td>13:49:53</td><td>52.08201</td><td>140</td><td>3.2</td><td>1861</td><td>10.206</td><td>1.204</td><td>-0.077</td><td>0.075</td><td>-28</td><td>46</td><td>29.3</td></tr><tr><td>26 Nov</td><td>10:09:09</td><td>52.92871</td><td>150</td><td>3.0</td><td>1922</td><td>11.393</td><td>-0.791</td><td>-0.176</td><td>0.073</td><td>71</td><td>44</td><td>26.9</td></tr><tr><td>27 Nov</td><td>11:36:33</td><td>53.98941</td><td>120</td><td>3.9</td><td>1866</td><td>12.879</td><td>-0.590</td><td>-0.087</td><td>0.075</td><td>-44</td><td>46</td><td>52.7</td></tr><tr><td>28 Nov</td><td>11:25:28</td><td>54.98171</td><td>110</td><td>4.0</td><td>1849</td><td>14.270</td><td>0.491</td><td>-0.019</td><td>0.076</td><td>-59</td><td>46</td><td>37.4</td></tr><tr><td>29 Nov</td><td>08:19:32</td><td>55.85260</td><td>140</td><td>3.1</td><td>1894</td><td>15.491</td><td>0.224</td><td>-0.016</td><td>0.074</td><td>26</td><td>45</td><td>38.7</td></tr><tr><td>29 Nov</td><td>11:15:55</td><td>55.97508</td><td>140</td><td>3.3</td><td>1870</td><td>15.662</td><td>1.007</td><td>0.052</td><td>0.075</td><td>-129</td><td>46</td><td>42.1</td></tr><tr><td>30 Nov</td><td>07:30:58</td><td>56.81887</td><td>150</td><td>3.2</td><td>1863</td><td>16.845</td><td>0.508</td><td>0.184</td><td>0.075</td><td>-199</td><td>46</td><td>44.7</td></tr><tr><td>01 Dec</td><td>08:19:49</td><td>57.85279</td><td>140</td><td>3.2</td><td>1879</td><td>18.294</td><td>0.273</td><td>0.187</td><td>0.075</td><td>-107</td><td>45</td><td>47.2</td></tr><tr><td>01 Dec</td><td>11:18:25</td><td>57.97681</td><td>130</td><td>3.4</td><td>1909</td><td>18.468</td><td>0.158</td><td>0.084</td><td>0.074</td><td>40</td><td>45</td><td>44.1</td></tr><tr><td>02 Dec</td><td>07:48:41</td><td>58.83116</td><td>150</td><td>3.1</td><td>1887</td><td>19.665</td><td>1.068</td><td>0.097</td><td>0.074</td><td>-164</td><td>45</td><td>45.9</td></tr><tr><td>03 Dec</td><td>09:55:37</td><td>59.91929</td><td>150</td><td>3.0</td><td>1899</td><td>21.190</td><td>1.147</td><td>0.082</td><td>0.074</td><td>51</td><td>45</td><td>30.4</td></tr><tr><td>17 Jan</td><td>09:19:04</td><td>104.89186</td><td>130</td><td>3.5</td><td>1759</td><td>84.221</td><td>0.200</td><td>-0.070</td><td>0.080</td><td>-45</td><td>49</td><td>34.0</td></tr><tr><td>18 Jan</td><td>05:01:52</td><td>105.71318</td><td>140</td><td>3.2</td><td>1816</td><td>85.372</td><td>-0.500</td><td>-0.144</td><td>0.077</td><td>-15</td><td>47</td><td>24.5</td></tr><tr><td>19 Jan</td><td>05:02:31</td><td>106.71356</td><td>140</td><td>3.4</td><td>1772</td><td>87.774</td><td>0.594</td><td>-0.140</td><td>0.079</td><td>-36</td><td>48</td><td>57.9</td></tr><tr><td>20 Jan</td><td>07:55:33</td><td>107.83363</td><td>100</td><td>4.8</td><td>1708</td><td>88.344</td><td>-0.478</td><td>-0.078</td><td>0.082</td><td>-48</td><td>50</td><td>26.6</td></tr><tr><td>21 Jan</td><td>05:04:22</td><td>108.71467</td><td>140</td><td>3.4</td><td>1792</td><td>89.579</td><td>0.613</td><td>-0.067</td><td>0.078</td><td>71</td><td>48</td><td>37.6</td></tr><tr><td>22 Jan</td><td>05:04:03</td><td>109.71438</td><td>120</td><td>4.1</td><td>1738</td><td>90.980</td><td>-0.937</td><td>0.068</td><td>0.081</td><td>-201</td><td>49</td><td>44.0</td></tr><tr><td>23 Jan</td><td>06:06:31</td><td>110.75767</td><td>140</td><td>3.3</td><td>1802</td><td>92.442</td><td>0.376</td><td>0.190</td><td>0.078</td><td>1</td><td>47</td><td>38.9</td></tr><tr><td>24 Jan</td><td>05:05:28</td><td>111.71519</td><td>140</td><td>3.2</td><td>1780</td><td>93.784</td><td>0.944</td><td>0.102</td><td>0.079</td><td>-127</td><td>48</td><td>46.4</td></tr><tr><td>25 Jan</td><td>06:30:41</td><td>112.77428</td><td>140</td><td>3.3</td><td>1805</td><td>95.269</td><td>-0.014</td><td>0.169</td><td>0.078</td><td>27</td><td>47</td><td>37.0</td></tr><tr><td>26 Jan</td><td>06:03:54</td><td>113.75560</td><td>140</td><td>3.5</td><td>1767</td><td>96.644</td><td>0.778</td><td>0.100</td><td>0.079</td><td>-51</td><td>48</td><td>44.9</td></tr><tr><td>27 Jan</td><td>06:58:50</td><td>114.79365</td><td>140</td><td>3.4</td><td>1774</td><td>98.099</td><td>-1.185</td><td>-0.011</td><td>0.079</td><td>-2</td><td>48</td><td>39.7</td></tr><tr><td>28 Jan</td><td>06:59:12</td><td>115.79383</td><td>140</td><td>3.4</td><td>1737</td><td>99.501</td><td>0.548</td><td>-0.019</td><td>0.081</td><td>70</td><td>49</td><td>41.8</td></tr><tr><td>29 Jan</td><td>06:05:30</td><td>116.75644</td><td>130</td><td>3.5</td><td>1758</td><td>100.850</td><td>0.958</td><td>0.062</td><td>0.080</td><td>-71</td><td>49</td><td>41.7</td></tr></tbody></table>

main goals of the MaTYSSE (Magnetic Topologies of Young Stars and the Survival of close-in massive Exoplanets) large-programme allocated on the 3.6m Canada-France-Hawaii Telescope (CFHT), thanks to which the youngest hJ discovered so far was detected (Donati et al. 2016, 2017) and within which this study places.

In this paper, we present results for another wTTS, the young pre-main sequence (PMS) solar-mass star, TAP 26, (Feigelson et al. 1987; Grankin et al. 2008; Grankin 2013), located in the Taurus star-forming region. TAP 26 was observed in late 2015 and early 2016 with both the ESPaDOnS spectropolarimeter and the 1.25 m telescope at the Crimean Astrophysical Observatory (CrAO). After documenting our observations (Sec. 2), we derive the stellar parameters of TAP 26 (Sec. 3), before reconstructing the surface magnetic and brightness maps by applying Zeeman-Doppler Imaging (ZDI) to our data (Sec. 4). We finally detail in Sec. 5 our detection of a planet radial velocity (RV) signal in its spectrum, using three different methods. The first two methods

are based on ZDI following previous studies (Donati et al. 2015, 2016; Petit et al. 2015), and the third one exploits Gaussian-Process Regression (GPR, Haywood et al. 2014; Rajpaul et al. 2015, see Sec. 5).

## 2 OBSERVATIONS

TAP 26 was observed in November 2015 and January 2016 using the high-resolution spectropolarimeter ESPaDOnS at the 3.6-m CFHT at Mauna Kea (Hawaii). ESPaDOnS collects stellar spectra spanning the entire optical domain (from 370 to 1,000 nm) at a resolving power of 65,000 (i.e., resolved velocity element of 4.6 km s⁻¹) over the full wavelength range (Donati 2003). A total of 29 unpolarised (Stokes I) and circularly-polarised (Stokes V) spectra were collected over a timespan of 72 d, 16 spectra over 16 nights in 2015 Nov, and 13 spectra over 13 nights in 2016 Jan. The rate was of one spectrum per night, except at the beginning of the 2015
<table><caption>**Table 2.** Journal of contemporaneous CrAO multicolour photometric observations of TAP 26 collected in late 2015 and early 2016, respectively listing the UT date and Heliocentric Julian Date (HJD) of the observation, the measured V magnitude ($1\sigma$ error bar of 0.016 mag) and V-R$_J$ Johnson photometric colours, and the corresponding rotational phase (using again the ephemeris given by Eq. 1). The table is divided into three periods spanning 1.5-2.5 months each, the second one covering the 2015 Nov set of spectropolarimetric observations and the third one overlapping the 2016 Jan set of spectropolarimetric observations.</caption><thead><tr><th>Date</th><th>HJD<br>(2,457,200+)</th><th>V<br>(mag)</th><th>V-R<sub>J</sub></th><th>Cycle<br>(-120+)</th></tr></thead><tbody><tr><td>25 Aug</td><td>60.569</td><td>12.291</td><td>-</td><td>1.946</td></tr><tr><td>30 Aug</td><td>65.592</td><td>12.269</td><td>0.986</td><td>8.987</td></tr><tr><td>31 Aug</td><td>66.583</td><td>12.261</td><td>1.010</td><td>10.375</td></tr><tr><td>09 Sep</td><td>75.557</td><td>12.297</td><td>1.016</td><td>22.953</td></tr><tr><td>11 Sep</td><td>77.562</td><td>12.331</td><td>1.022</td><td>25.763</td></tr><tr><td>16 Sep</td><td>82.564</td><td>12.329</td><td>1.004</td><td>32.774</td></tr><tr><td>18 Sep</td><td>84.594</td><td>12.259</td><td>1.004</td><td>35.619</td></tr><tr><td>19 Sep</td><td>85.530</td><td>12.300</td><td>1.007</td><td>36.930</td></tr><tr><td>22 Sep</td><td>88.529</td><td>12.260</td><td>1.003</td><td>41.134</td></tr><tr><td>23 Sep</td><td>89.505</td><td>12.245</td><td>1.014</td><td>42.501</td></tr><tr><td>24 Sep</td><td>90.517</td><td>12.282</td><td>0.988</td><td>43.920</td></tr><tr><td>25 Sep</td><td>91.550</td><td>12.246</td><td>0.988</td><td>45.369</td></tr><tr><td>26 Sep</td><td>92.524</td><td>12.320</td><td>1.001</td><td>46.733</td></tr><tr><td>28 Sep</td><td>94.550</td><td>12.238</td><td>0.968</td><td>49.573</td></tr><tr><td>03 Oct</td><td>99.588</td><td>12.283</td><td>1.030</td><td>56.633</td></tr><tr><td>04 Oct</td><td>100.513</td><td>12.276</td><td>0.983</td><td>57.930</td></tr><tr><td>09 Oct</td><td>105.545</td><td>12.280</td><td>1.016</td><td>64.982</td></tr><tr><td>15 Oct</td><td>111.600</td><td>12.232</td><td>0.967</td><td>73.469</td></tr><tr><td>16 Oct</td><td>112.605</td><td>12.292</td><td>0.976</td><td>74.877</td></tr><tr><td>17 Oct</td><td>113.595</td><td>12.269</td><td>1.000</td><td>76.265</td></tr><tr><td>19 Oct</td><td>115.597</td><td>12.261</td><td>0.984</td><td>79.070</td></tr><tr><td>20 Oct</td><td>116.584</td><td>12.233</td><td>0.963</td><td>80.454</td></tr><tr><td>25 Oct</td><td>121.564</td><td>12.263</td><td>1.014</td><td>87.434</td></tr><tr><td>27 Oct</td><td>123.507</td><td>12.247</td><td>0.994</td><td>90.157</td></tr><tr><td>30 Oct</td><td>126.442</td><td>12.280</td><td>1.024</td><td>94.270</td></tr><tr><td>03 Nov</td><td>130.564</td><td>12.220</td><td>1.012</td><td>100.048</td></tr><tr><td>13 Nov</td><td>140.585</td><td>12.229</td><td>0.989</td><td>114.092</td></tr><tr><td>16 Dec</td><td>173.373</td><td>12.245</td><td>1.003</td><td>160.046</td></tr><tr><td>17 Dec</td><td>174.306</td><td>12.238</td><td>0.979</td><td>161.354</td></tr><tr><td>03 Jan</td><td>191.364</td><td>12.215</td><td>0.976</td><td>185.262</td></tr><tr><td>17 Jan</td><td>205.347</td><td>12.306</td><td>0.983</td><td>204.860</td></tr><tr><td>24 Jan</td><td>212.316</td><td>12.245</td><td>1.009</td><td>214.626</td></tr><tr><td>30 Jan</td><td>218.296</td><td>12.297</td><td>1.019</td><td>223.008</td></tr><tr><td>10 Feb</td><td>229.258</td><td>12.217</td><td>0.975</td><td>238.371</td></tr><tr><td>22 Feb</td><td>241.262</td><td>12.245</td><td>0.982</td><td>255.195</td></tr><tr><td>05 Mar</td><td>253.253</td><td>12.293</td><td>0.987</td><td>272.002</td></tr><tr><td>08 Mar</td><td>256.285</td><td>12.238</td><td>0.992</td><td>276.251</td></tr><tr><td>15 Mar</td><td>263.268</td><td>12.299</td><td>1.002</td><td>286.038</td></tr></tbody></table>

Nov session where a three-day gap following the first observation was compensated by pairs of observations on Nov 25, Nov 29 and Dec 01. However, given the 0.71 d rotation period of TAP 26, phase coverage is not optimal and the 2015 Nov data set presents gaps of 0.15-0.25 rotation cycle (see Table 1).

Each polarisation exposure sequence consists of 4 individual subexposures taken in different polarimeter configurations to allow the removal of all spurious polarisation signatures at first order. All raw frames are processed with the nominal reduction package LIBRE ESPRIT as described in the previous papers of the series (e.g., Donati et al. 2010, 2011, 2014), yielding a typical rms RV precision of 20-30 m s⁻¹ (Moutou et al. 2007; Donati et al. 2008). The peak signal-to-noise ratios (S/N, per 2.6 km s⁻¹ velocity bin) achieved on the collected spectra range between 100 and 150 (median 140), depending mostly on weather/seeing conditions. The full journal of observations is presented in Table 1.

Rotational cycles (noted *E* in the following equation) are computed from Barycentric Julian Dates (BJDs) according to the ephemeris:

$$BJD (d) = 2,457,344.8 + P_{\text{rot}}E \quad (1)$$

in which the photometrically-determined rotation period $P_{\text{rot}}$ (equal to 0.7135 d, Grankin 2013) is taken from the literature and the initial Julian date (2,457,344.8 d) is chosen arbitrarily.

Least-Squares Deconvolution (LSD, Donati et al. 1997) was applied to all spectra. The line list we employed for LSD is computed from an ATLAS9 LTE model atmosphere (Kurucz 1993) featuring $T_{\text{eff}}=4,500$ K and $\log g=4.5$, the most appropriate model for TAP 26 (see Sec. 3). Only moderate to strong atomic spectral lines are included in this list (see, e.g., Donati et al. 2010, for more details). Altogether, about 7,800 spectral features (with about 40% from Fe I) are used in this process. The Stokes *I* and Stokes *V* LSD profiles can be seen in Sec. 4. Significant distortions are visible in all Stokes *I* LSD profiles, indicating the presence of brightness inhomogeneities covering a large fraction of the surface of TAP 26 at the time of our observations. The noise level in Stokes *I* LSD profiles is measured from continuum intervals (see Table 1), and includes not only the noise from photon statistics, but also the (often dominant) noise introduced by LSD.

Among the 29 profiles we used, 11 were contaminated by solar light reflected off the Moon (5 in 2015 Nov, the Moon being at 9.5° from TAP 26 and at 99% illumination on 2015 Nov 26, and 6 in 2016 Jan, the Moon being at 12° from TAP 26 and at 85% illumination on 2016 Jan 19); we applied a two-step process involving tomographic imaging, described in Donati et al. (2016), to filter out this contamination from our Stokes *I* LSD profiles.

Regarding the Stokes *V* profiles, Zeeman signatures are detected in all observations, featuring amplitudes of typically 0.1%. Expressed in units of the unpolarised continuum level $I_c$, the average noise levels of the Stokes *V* LSD signatures (dominated here by photon statistics) range from 2.3 $10^{-4}$ to 3.9 $10^{-4}$ per 1.8 km s⁻¹ velocity bin - with a median value of 2.8 $10^{-4}$.

The emission core of the Ca II Infrared Triplet (IRT) lines exhibit an average equivalent width of ≈10 km s⁻¹, corresponding to the amount expected from chromospheric emission for such a wTTS. The He I $D_3$ line is relatively faint (average equivalent width of ≈5 km s⁻¹), demonstrating that accretion is no longer taking place at its surface, in agreement with previous studies (Donati et al. 2014, 2015). The $H\alpha$ line is also relatively weak by wTTS standards (Kenyon & Hartmann 1995), with an average equivalent width of 40 km s⁻¹, and is modulated with a period of 0.7132±0.0002 d (see App. B).

Contemporaneous VR$_J$ photometric observations were also collected from the Crimean Astrophysical Observatory
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Observed location of TAP 26 in the HR diagram. The red and blue open squares (with $1\sigma$ error bars) depict the location of TAP 26 using 2 different ways of estimating the inclination angle of the rotation axis - with the red one showing our best estimate used throughout the paper. The PMS evolutionary tracks for 0.9 $M_\odot$, 1.0 $M_\odot$, 1.1 $M_\odot$, 1.2 $M_\odot$ and 1.3 $M_\odot$, and corresponding isochrones for 10 Myr, 15 Myr and 20 Myr (Siess et al. 2000) assume solar metallicity and include convective overshooting. The green lines depict where models predict PMS stars' radiative core reaches a radius of 0.5 $R_\star$ and 0.6 $R_\star$.</figcaption></figure>

(CrAO) 1.25 m telescope between 2015 Aug and 2016 Mar. They indicate a brightness modulation with a period of 0.7138±0.0001 d of full amplitude 0.116 mag in V (see Table 2). By analogy with other wTTSs, these photometric variations can be safely attributed to the presence of brightness features at the surface of TAP 26 modulated by rotation. The small difference with the value found in Grankin (2013) suggests the presence of differential rotation in TAP 26 (see Sec. 4).

## 3 EVOLUTIONARY STATUS OF TAP 26

TAP 26 is a well-studied single wTTSs, close enough to T Tau, both spatially and in terms of velocity, to assume a distance of 147±3 pc (Loinard et al. 2007; Torres et al. 2009), with an error bar similar to that found on other regions of Taurus like L1495.

Applying the automatic spectral classification tool especially developed in the context of MaPP (Magnetic Planets and Planets) and MaTYSSE, following that of Valenti & Fischer (2005) and discussed in Donati et al. (2012), we find that the photospheric temperature and logarithmic gravity of TAP 26 are respectively equal to $T_{\text{eff}}=4,620\pm50$ K and $\log g=4.5\pm0.2$ (with g in cgs units). This is warmer than the temperature quoted in the literature (4340 K, Grankin 2013), which is derived from photometry and thus expected to be significantly less accurate than ours, derived from high-resolution spectroscopic data, enabling to find the actual temperature without the disturbance of circumstellar and interstellar reddening.

Long-term photometric monitoring of TAP 26 indicates that its maximum V magnitude is equal to 12.16 (Grankin et al. 2008). Following Donati et al. (2014, 2015), we assume

<table><caption>**Table 3.** Parameters for TAP 26, inferred from the photometric and spectroscopic measurements and the ZDI analysis (see Sec. 4). Respectively: distance to Earth *d*, mass $M_\star$, radius $R_\star$, effective temperature $T_{\text{eff}}$, decimal logarithm of surface gravity $\log g$, logarithmic luminosity $\log(L_\star/L_\odot)$, age, rotation period $P_{\text{rot}}$, inclination of the rotation axis to the line of sight *i*, line-of-sight-projected equatorial rotation velocity $v\sin i$, equatorial rotation rate $\Omega_{\text{eq}}$, difference $d\Omega$ between equatorial and polar rotation rates and mean radial velocity in the barycentric rest frame $v_{\text{rad}}$ (which was derived from our spectropolarimetric runs, see Sec. 4). T09 and G13 in the references respectively stand for Torres et al. (2009) and Grankin (2013).</caption><thead><tr><th>Parameter</th><th>Value</th><th>Reference</th></tr></thead><tbody><tr><td><i>d</i> (pc)</td><td>147±3</td><td>T09</td></tr><tr><td><i>M</i><sub>⋆</sub> (M<sub>☉</sub>)</td><td>1.04±0.10</td><td></td></tr><tr><td><i>R</i><sub>⋆</sub> (R<sub>☉</sub>)</td><td>1.17±0.17</td><td></td></tr><tr><td><i>T</i><sub>eff</sub> (K)</td><td>4,620±50</td><td></td></tr><tr><td>log <i>g</i></td><td>4.5</td><td></td></tr><tr><td>log(<i>L</i><sub>⋆</sub>/<i>L</i><sub>☉</sub>)</td><td>-0.25±0.10</td><td></td></tr><tr><td>Age (Myr)</td><td>≈17</td><td></td></tr><tr><td><i>P</i><sub>rot</sub> (d)</td><td>0.7135</td><td>G13</td></tr><tr><td><i>i</i> (°)</td><td>55±10</td><td></td></tr><tr><td>ν sin <i>i</i> (km s<sup>-1</sup>)</td><td>68.2±0.5</td><td></td></tr><tr><td>Ω<sub>eq</sub> (rad d<sup>-1</sup>)</td><td>8.8199±0.0003</td><td></td></tr><tr><td><i>dΩ</i> (rad d<sup>-1</sup>)</td><td>0.0492±0.0010</td><td></td></tr><tr><td>ν<sub>rad</sub> (km s<sup>-1</sup>)</td><td>16.25±0.20</td><td></td></tr></tbody></table>

a spot coverage<fnref n="1" /> of ≈25% at maximum brightness, typical for active stars (and caused by, e.g., the presence of high-latitude cool spots and / or of small spots evenly spread over the whole stellar surface), we derive an unspotted V magnitude of 11.86±0.20. From the difference between the B-V index expected at the temperature of TAP 26 (equal to 0.99±0.02, Pecaut & Mamajek 2013) and the averaged value measured for TAP 26 (equal to 1.13±0.05, see Kenyon & Hartmann 1995; Grankin et al. 2008), and given the very weak impact of starspot on B-V (Grankin et al. 2008), we derive that the amount of visual extinction $A_V$ that our target suffers is equal to 0.43±0.15 (within 1.5σ of the value of Herczeg & Hillenbrand 2014, despite the very different methods used to estimate this parameter). Using the visual bolometric correction expected for the adequate photospheric temperature (equal to -0.55±0.05, see Pecaut & Mamajek 2013) and the distance estimate assumed previously (147±3 pc), corresponding to a distance modulus of 5.84±0.04, we finally obtain a bolometric magnitude of 5.04±0.26, or equivalently a logarithmic luminosity relative to the Sun of -0.12±0.10. Coupling with the photospheric temperature obtained previously, we find a radius of 1.36±0.17 $R_\odot$ for our target star.

The rotation period of TAP 26 is well determined from long-term multi-colour photometric monitoring, with an average value over the full data set equal to 0.7135 d (Grankin 2013). Coupling this rotation period along with our measurements of the line-of-sight-projected equatorial rotation velocity $v\sin i$ of TAP 26 (equal to 68.2±0.5 km s⁻<fnref n="1" />, see Sec. 4), we can infer that $R_\star\sin i=0.96\pm0.05 R_\odot$, where $R_\star$ and $i$ de-

<footnote n="1">¹ Spot coverage: integral of the difference between local brightness and photosphere brightness over the surface of the star, in units of photosphere brightness.</footnote>
note the radius of the star and the inclination of its rotation axis to the line of sight. Comparing with the radius derived from the luminosity and photometric temperature, we derive that $i=45\pm8^\circ$.

Using ZDI, we actually infer from our data that $i=55\pm10^\circ$(see Sec. 4). The $1\sigma$ difference with the previous estimate can be simply interpreted as an over-estimate in spottedness at maximum brightness. Assuming now a spottedness of 12% at maximum brightness (instead of 25%) reconciles both approaches and yields a logarithmic luminosity of $-0.25\pm0.10$ and thus a radius of $1.17\pm0.17 R_\odot$, in good agreement with other studies ($1.18 R_\odot$ in Herczeg & Hillenbrand 2014).

Using the evolutionary models of Siess et al. (2000) (assuming solar metallicity and including convective overshooting), we obtain that TAP 26 is a $\approx 17$ Myr star (in good agreement with the estimate of Grankin 2013) and that its mass is $M_*=1.04\pm0.10 M_\odot$ (see Fig 1). The average equivalent width of the 670.7 nm Li line is equal to 0.045 nm, in good agreement with that measured for solar-mass PMS stars in the 10-15 Myr Sco-Cen association at the corresponding temperature (Pecaut & Mamajek 2016), which further confirms our age estimate and thus the evolutionary status of TAP 26.

Referring to Donati et al. (2015, 2016), TAP 26 closely resembles an evolved version of the 2 Myr star V830 Tau that would have contracted and spun up by 4x towards the zero-age main sequence, with the rotation period and radius of V830 Tau being respectively 2.741 d and $2.0\pm0.2 R_\odot$. The increase in rotation rate matches quite well the predicted decrease in the moment of inertia between both epochs according to evolutionary models of Siess et al. (2000). Given the prominent role of the disc in braking the rotation of the star and thus decreasing its angular momentum (Gallet & Bouvier 2015; Davies et al. 2014), this also suggests that TAP 26 dissipated its accretion disc very early, typically as early as, or earlier than V830 Tau. We also note that our target is located past the theoretical threshold at which stars start to be more than half radiative in radius, suggesting that the magnetic field of TAP 26 already started to evolve into a complex topology (Gregory et al. 2012).

The stellar parameters inferred and used in this study are summarised in Table 3.

# 4 TOMOGRAPHIC IMAGING

In order to model the activity jitter of TAP 26 (see Sec. 5), we applied ZDI (Semel 1989; Brown et al. 1991; Donati & Brown 1997) to our data. ZDI takes inspiration from medical tomography, which consists of constraining a 3D distribution using series of 2D projections as seen from various angles (Vogt et al. 1987). In our context, ZDI inverts simultaneous time series of 1D Stokes I and V LSD profiles into 2D brightness and magnetic field maps of the stellar surface (see Donati et al. 2014). The magnetic field is decomposed into its poloidal and toroidal components, both expressed as spherical harmonics expansions (Donati et al. 2006).

Synthetic LSD profiles are derived from brightness and magnetic maps by summing up the spectral contribution of all cells, taking into account the Doppler broadening caused by the rotation of the star, the Zeeman effect induced by magnetic fields and the continuum center-to-limb darkening.

Local Stokes I and V profiles are computed using Unno-Rachkovsky's analytical solution to the polarised radiative transfer equations in a Milne-Eddington model atmosphere (Landi degl'Innocenti & Landolfi 2004). The local profile used in this study has a central wavelength, a Doppler width and a Landé factor of typical values 670 nm, $1.8 \text{ km s}^{-1}$ and $1.2$ respectively, and an equivalent width of $4.6 \text{ km s}^{-1}$ corresponding to the LSD profiles of TAP 26. Technically, ZDI applies a conjugate gradient technique to iteratively reconstruct the brightness and magnetic surface maps with minimal information content (i.e. maximum Shannon entropy) that matches our observed LSD profiles at a given reduced chi-square ($\chi^2$, defined as $\chi^2$ divided by the number of data points$^2$) level. Concerning the brightness, we note that, unlike in Donati & Collier Cameron (1997) where we fit a spot filling factor with preset spot parameters, here we fit the local brightness $b_k$ of cell k, relative to the quiet photosphere ($0<b_k<1$ for dark spots and $b_k>1$ for bright plages), as described in Donati et al. (2014).

ZDI can also take into account and model latitudinal differential rotation, shearing the brightness distribution and magnetic topology at the surface of the star, and assuming a solar-like surface rotation rate, $\Omega(\theta)$, varying with latitude, $\theta$, as:

$$ \Omega(\theta) = \Omega_{eq} - d\Omega(\sin\theta)^2 \quad (2) $$

where $\Omega_{eq}$ is the equatorial rotation rate and $d\Omega$ is the difference between the equatorial and the polar rotation rates.

For a given set of parameters, ZDI looks for the map with minimal information content that matches the LSD profiles at $\chi^2=1$. As a by-product, we obtain the optimal stellar parameters for which the reconstructed images contain minimal information: $i=55\pm10^\circ$, $v\sin i=68.2\pm0.5 \text{ km s}^{-1}$ and $v_{rad}=17.0\pm0.2 \text{ km s}^{-1}$ (the RV the star would have if unspotted and planet-free). Regarding differential rotation, we obtain $\Omega_{eq}=8.8199\pm0.0003 \text{ rad d}^{-1}$ and $d\Omega=0.0492\pm0.0010 \text{ rad d}^{-1}$, as outlined in more detail in Sec. 4.2.

## 4.1 Brightness and magnetic imaging

Given the long time span between our two data sets (about 60 d, see Table 1), we start by reconstructing separate brightness and magnetic maps for each epoch (2015 Nov and 2016 Jan), before investigating the temporal variability between both in more detail.

The Stokes I and V LSD profiles, which are displayed in Fig. 2, were used simultaneously to reconstruct both surface brightness and magnetic field maps. The synthetic LSD profiles presented in the figure match the observed ones at $\chi^2=1$, or, equivalently, at a $\chi^2$ equal to 1484 for the 2015 Nov dataset and 1157 for the 2016 Jan dataset, and for both sets of Stokes I and V LSD profiles. The iterative reconstruction starts from unspotted magnetic maps corresponding to $\chi^2=13$ (2015 Nov) and 9 (2016 Jan), showing that the iterative algorithm of ZDI successfully manages to reproduce the

<footnote>$^2$ This follows the usual convention in regularised tomographic imaging techniques, where the number of model parameters is much smaller than the number of fitted data points and not taken into account in the expression of $\chi^2$ (Donati et al. 2017).</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Maximum entropy fit (thin red lines) to the observed (thick black lines) Stokes *I* (left) and *V* (right) LSD profiles. The 2015 Nov dataset is represented in the 1st and 3rd panels and the 2016 Jan dataset in the 2nd and 4th panels. The Stokes *I* LSD profiles before the removal of lunar pollution are coloured in cyan, and 3σ-error bars are displayed for the Stokes *V* profiles. The rotational cycles are written beside their corresponding profiles, in concordance with Table 1.</figcaption></figure>

data at noise level. In the particular case of Stokes *I* profiles, whose noise includes a significant level of systematics (see Sec. 2), we find that smaller error bars make ZDI unable to fit the data down to $\chi_r^2=1$; on the opposite, greater error bars result in a fit to the Stokes *I* profiles for which the raw radial velocities are not properly reproduced (see Sec. 5). This gives us confidence that the S/N values derived for the Stokes *I* LSD profiles (see Table 1) are accurate and reliable within 10%.

The reconstructed brightness maps for 2015 Nov and 2016 Jan are shown in Fig. 3, at an epoch corresponding to rotation cycle 10.0 (in the ephemeris of Eq. 1) for 2015 Nov, and 92.0 for 2016 Jan (see Table 1); the colour scale codes the logarithmic relative brightness compared to that of the photosphere. The surface spot coverage we derive is similar at both epochs, reaching 10% in the 2015 Nov map (5% / 5% of cool spots / hot plages respectively) and 12% in the 2016 Jan map (7% / 5% of cool spots / hot plages respectively). Both reconstructed maps share some similarities, such as a large cool polar cap resembling that reconstructed on other rapidly rotating wTTSs (e.g. Skelly et al. 2010; Donati et al. 2014), plus a number of smaller features located at lower latitudes (in particular the two equatorial spots located at phases 0.22 and 0.92 in 2015 Nov, 0.27 and 0.97 in 2016 Jan) interleaved with bright plages. We stress that ZDI is only sensitive to the medium and large bright-

ness features and misses small spots evenly distributed over the whole stellar surface, implying that the spottedness we recover for TAP26 is likely an underestimate. We observe a number of differences between both images potentially attributable to differential rotation and/or intrinsic variability (see Sec. 4.2); however, the limited phase coverage at both epochs makes the direct comparison of individual surface features between maps ambiguous and hazardous. We caution that the smallest-scale structures may reflect to some extent the limited phase coverage and be subject to phase ghosting (e.g. Stout-Batalha & Vogt 1999).

Using the brightness maps reconstructed with ZDI, we can predict photometric light curves at both epochs, which are found to compare well with our contemporaneous CrAO observations (see Fig. 4). Note the small but significant temporal evolution of the light-curve that we predict between the two epochs; this variability is however not obvious from the observed photometric data given their limited sampling and comparatively large error bars (rms 16 mmag).

The reconstructed magnetic topology is shown in Fig. 5. The large-scale field reconstructed for TAP 26 features a rms magnetic flux of 330 and 430 G in 2015 Nov and 2016 Jan respectively. The field is found to be mainly poloidal (70% of the reconstructed magnetic energy), though with a significant toroidal component (30% of the reconstructed
<figure><img src="image_4.png" /></figure>

**Figure 3.** Flattened polar view of the surface brightness maps for the 2015 Nov dataset (top panel) and 2016 Jan dataset (bottom panel). The equator and the 60°, 30° and -30° latitude parallels are depicted as solid and dashed black lines respectively. The colour scale indicates the logarithm of the relative brightness, with brown/blue areas representing cool spots/bright plages. Finally, the outer ticks mark the phases of observation.

<figure><img src="image_5.png" /></figure>

**Figure 4.** Photometry curves of the relative brightness as function of the rotation phase. The light curves synthesised from the reconstructed brightness maps for 2015 Nov and 2016 Jan are represented by a dashed red line and a dotted blue line respectively. The CrAO measurements are represented as dots with $1\sigma$ error bars, with the observations from 2015 Aug to 2015 Oct in black circles, the observations from 2015 Oct to 2015 Dec in red upward-pointing triangles and the observations from 2015 Dec to 2016 Mar in blue downward-pointing triangles.

magnetic energy). It is also largely axisymmetric (50% and 80% of the poloidal and toroidal field energy respectively).

The dipolar component of the large-scale field has a strength of 120±10 G at both epochs, corresponding to about 10% of the reconstructed poloidal field energy, and is tilted at 40±5° to the line of sight, i.e., midway to the equator, towards phase 0.73±0.03 and 0.85±0.03 in 2015 Nov and 2016 Jan respectively. The increase in the phase towards which the dipole is tilted suggests that intermediate to high latitudes (at which the dipole poles are anchored) are rotating more slowly than average by 0.19%, i.e., with a period of ≈0.7148 d; this is confirmed by the fact that the line-of-sight projected (longitudinal) magnetic field $B_\ell$ (proportional to the first moment of the Stokes V profiles, e.g., Donati et al. 1997, and most sensitive to the low-order components of the large-scale field) exhibit a recurrence timescale of 1.0014±0.003 $P_{\text{rot}}$ (see App. B), i.e., slightly longer than $P_{\text{rot}}$ by a similar amount. Higher order terms in the spherical harmonics expansion describing the field (in particular the quadrupolar and octupolar modes) get stronger between 2015 Nov and 2016 Dec, with total magnetic energies increasing from 85% to 93% of the poloidal field.

Finally, we show a large-scale extrapolation of the magnetic field (under the assumption of a potential field) in Fig. 6. Similarly to the brightness maps, the magnetic maps seem to point to a variation of the surface topology between late 2015 and early 2016, which is not explained by differential rotation alone, though the limited phase coverage calls for caution when comparing features between those maps.

The magnetic maps suggest that the magnetic topology at the rotation pole underwent a ≈0.1 phase shift between both dates.

## 4.2 Intrinsic variability and differential rotation

When applying ZDI to the whole dataset, i.e., modelling all Stokes I and V profiles with only one brightness map and one magnetic topology (see App. A), we obtain a minimum $\chi_r^2$ value of 1.4, even when taking into account differential rotation (starting from an initial value $\chi_r^2=20$). This indicates that intrinsic variability occurred during the 45 d gap (or 63 rotation cycles) separating both datasets.

Despite this variability, we attempted to retrieve differential rotation from the whole data set. The search for differential rotation parameters is done by minimising the value of $\chi_r^2$ at a fixed amount of information, in this present case using the Stokes I profiles and brightness map reconstruction only. From the curvature of the $\chi_r^2$ paraboloid around the minimum, one can infer error bars on differential rotation parameters (Donati et al. 2003). The spot coverage is fixed at 13% (chosen to be slightly higher than the values found in each reconstruction) and the values we found are $\Omega_{\text{eq}}=8.8199\pm0.0003 \text{ rad d}^{-1}$ and $d\Omega=0.0492\pm0.0010 \text{ rad d}^{-1}$, with a minimum $\chi_r^2$ of 1.4116. A map of $\Delta\chi^2$ is shown in Fig. 7, which presents a very clear paraboloid around the minimum we found, even if, due to our phase coverage, these precise values ask for further confirmation with the help of future data. This value of $d\Omega$ is close to the solar differential rotation (0.055rad d$^{-1}$). The case with no differential rotation yields $\chi_r^2=2.6907$. Normalising $\Delta\chi^2$ by the minimum $\chi^2$ achieved over the map (to scale up error bars as a way to account for the contribution from the reported intrinsic vari-
<figure><img src="image_6.png" /><figcaption>**Figure 5.** From left to right: radial, meridional and azimuthal component of the surface magnetic field (labelled in G), reconstructed with ZDI from the 2015 Nov dataset (top panels) and the 2016 Jan dataset (bottom panels).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Potential field extrapolations of the reconstructed magnetic topology as seen by an Earth-based observer, in 2015 Nov (left) and in 2016 Jan (right) both at phase 0.8. Open and closed field lines are shown in blue and white respectively, whereas colours at the stellar surface depict the local value of the radial field (in G, as shown in the left panels of Fig. 5). The source surface at which the field becomes radial is set to $4 R_\star$, slightly larger than the corotation radius of about $3 R_\star$ (at which the Keplerian period equals the stellar rotation period) and beyond which field lines are expected to quickly open under centrifugal forces.</figcaption></figure>
<figure><img src="image_8.png" /></figure>

**Figure 7.** Map of $\Delta\chi^2$ as a function of $\Omega_{\text{eq}}$ and $d\Omega$, derived from the modelling of our Stokes *I* LSD profiles of TAP 26 at constant information content. A well-defined paraboloid is observed with the outer colour contour corresponding to the 99.99% confidence level area (i.e., a $\chi^2$ increase of 18.4 for the 2581 Stokes *I* data points). The minimum value of $\chi_r^2$ is 1.4116. The minimum $\chi_r^2$ achieved is above unity due to intrinsic variability affecting the LSD profiles but not being taken into account within ZDI. The derived differential rotation parameters are $\Omega_{\text{eq}}=8.8199\pm0.0003 \text{ rad d}^{-1}$ and $d\Omega=0.0492\pm0.0010 \text{ rad d}^{-1}$.

ability) still yields a value in excess of 3300 and a negligible false alarm probability (FAP), unambiguously demonstrating that the star is not rotating as a solid body.

The differential rotation parameters we obtain imply a lap time of 128±3 d, with rotation periods of 0.71239±0.00003 d and 0.71638±0.00008 d for the equator and pole respectively, in good agreement with the range of rotation periods derived from photometry (ranging from 0.7135 to 0.7138, Grankin 2013). The 0.7132 d period found for the equivalent width of the $H\alpha$ line and the 0.7145 d period found for the longitudinal magnetic field $B_\ell$ (see App. B) are also consistent. We note that the rotation periods found with photometry, the longitudinal magnetic field and $H\alpha$ line correspond to latitudes ranging from 30° to 50°, indicating that an important amount of activity is concentrated at these mid-latitudes, with the dipole pole located in the upper part of this range, in good agreement with the ZDI reconstruction (see Sec. 4.1).

# 5 MODELLING THE PLANET SIGNAL

We describe below three different techniques aimed at characterising the RV curve of TAP 26. The first two are those used in Donati et al. (2016): filtering out the activity modelled with the help of ZDI, and the simultaneous fit of the planet parameters and the stellar activity. The third method follows the approach of Haywood et al. (2014) and Rajpaul et al. (2015) and uses Gaussian-Process Regression (GPR) to model the activity directly from the raw RVs. The results obtained from these three methods are outlined and discussed in the following sections.

## 5.1 Jitter activity filtering

The first technique consists of using the previously reconstructed maps to predict the pollution to the RV curve caused by activity (called activity jitter in the following) and subtract it from the raw RVs. From the observed Stokes *I* LSD profiles, we compute, at both epochs, the raw RVs RV<sub>raw</sub> (and error bars, see Table 1), as the first-order moment of the continuum-subtracted corresponding profiles (Donati et al. 2017). Likewise, the synthesised Stokes *I* LSD profiles derived from the brightness maps yield the synthesised activity jitter of the star (RV signal caused by the brightness distribution and stellar rotation). By subtracting the activity jitter from the raw RVs, we obtain filtered RVs RV<sub>filt</sub> (see Table 1). We observe that the jitter has a mean semi-amplitude of 1.81 km s<sup>-1</sup> in 2015 Nov and 1.21 km s<sup>-1</sup> in 2016 Jan, whereas the filtered RV curve features a signal with a semi-amplitude of ≈0.15 km s<sup>-1</sup> (Fig. 8), i.e., 8 to 12 times smaller than the activity signal we filtered out. We note the very significant evolution in the activity curve between 2015 Nov and 2016 Jan, demonstrating that the brightness distribution has evolved at the surface of TAP 26, through differential rotation and intrinsic variability (see Sec. 4).

With a rms dispersion of 109 m s<sup>-1</sup>, the filtered RVs clearly show the presence of a signal. Looking for a planet signature, we want to fit a sine curve (of semi-amplitude *K*, period $P_{\text{orb}}$, phase of inferior conjunction $\phi$, and offset RV<sub>0</sub>) to these filtered RVs, which corresponds to a circular orbit (see Fig. 9). The phase of inferior conjunction, i.e., corresponding to the epoch at which the planet is closest to us, is defined relatively to the reference date BJD<sub>c0</sub>=2,457,352.6485 (rotation cycle 11.0, approximately at the centre of the 2015 Nov observation run), such that the inferior conjunction occurs at BJD<sub>c</sub>=BJD<sub>c0</sub>+(φ-1)P<sub>orb</sub>. Due to the gap between both observing runs, several sine fits with different frequencies match the RV<sub>filt</sub> as local minima of $\chi_r^2$. The four best fits are shown in Fig. 9 and their characteristics are given in Table 4, with the value of the log likelihood as computed from the $\Delta\chi^2$ over these 29 RV data points. The residual RVs, derived from subtracting the best sine fit to the filtered RVs (shown in Fig. 9), feature a rms value of 51 m s<sup>-1</sup>.

Plotting Lomb-Scargle periodograms for the raw RVs, filtered RVs and residual RVs further demonstrates the presence of a periodic signal in the filtered RVs (Fig. 10). The above-mentioned dominant periods are seen as peaks in the periodogram; periodograms of partial data (only the 2015 Nov dataset, only the 2016 Jan dataset, odd points and even points) are also shown, yielding peaks at the same frequencies albeit with a lower power. We highlight the fact that the highest peaks in the raw RVs correspond to the activity jitter and are located at $P_{\text{rot}}/2$ and its aliases, whereas little power concentrates at $P_{\text{rot}}$ itself. A zoom-in of the filtered RV periodogram is also shown in Fig. 10 (bottom panel). The FAP is 0.06% for the highest peak ($P_{\text{orb}}$ = 13.41 d = 18.80 $P_{\text{rot}}$), and no significant period stands out in the residual RVs after filtering out both the activity jitter and the planet signal corresponding to the highest peak. We carried out simulations to ensure that the detected peaks are not generated by the filtering process, see details in App. C. Study of other activity proxies shows that the detected orbital periods are not present in the activity signal either (App. B).
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Top panels: RV (in the stellar rest frame) of TAP 26 as a function of rotation phase, as measured from our observations (open circles) and predicted by the tomographic maps (blue line). The synthesised raw RV curves exhibit only low-level temporal evolution resulting from differential rotation. Bottom panels: filtered RVs derived by subtracting the modelled activity jitter from the raw RVs, with a 10x zoom-in on the vertical axis.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Top: filtered RVs of TAP 26 and four sine curves representing the best fits. The thick green curve represents the case $P_{orb}/P_{rot}=18.80$, the thin magenta one $P_{orb}/P_{rot}=15.27$, the dash-and-dotted blue one $P_{orb}/P_{rot}=24.56$ and the dotted black one $P_{orb}/P_{rot}=12.76$. Bottom: residual RVs resulting from the subtraction of the best fit (green curve) from the filtered RVs. The residual RVs feature a rms value of 51 m s⁻¹.</figcaption></figure>

<table><caption>**Table 4.** Characteristics of the four best sine curve fits to the filtered RVs, and the case without planet. Respectively: semi-amplitude *K*, orbital period *P*<sub>orb</sub> in units of *P*<sub>rot</sub>, orbital period *P*<sub>orb</sub> in days, phase of inferior conjunction *φ* relative to cycle 11.0 (see ephemeris in Eq. 1), BJD of inferior conjunction, mean RV, corresponding *χ*<sup>2</sup>, difference in *χ*<sup>2</sup> with the best fit (Δ*χ*<sup>2</sup>, summed on the 29 data points), and natural logarithm (log<sub>e</sub>) of the likelihood *L*<sub>r1</sub> relative to the best fit. *φ* relates to the epoch of inferior conjunction BJD<sub>c</sub> through BJD<sub>c</sub>=2,457,352.6485+*φ*<sub>orb</sub>, the reference date being chosen so as to minimise the variation of *φ* between the four cases.</caption><thead><tr><th>K<br>(km s<sup>-1</sup>)</th><th>P<sub>orb</sub><br>(P<sub>rot</sub>)</th><th>P<sub>orb</sub><br>(d)</th><th>φ</th><th>BJD<sub>c</sub><br>(2,457,340+)</th><th>RV<sub>0</sub><br>(km s<sup>-1</sup>)</th><th>χ<sup>2</sup><sub>r</sub></th><th>Δχ<sup>2</sup></th><th>log L<sub>r1</sub></th><th>style<br>in Fig. 9</th></tr></thead><tbody><tr><td>0.131±0.020</td><td>18.80±0.23</td><td>13.41±0.16</td><td>0.709±0.026</td><td>8.75±0.35</td><td>0.009±0.014</td><td>0.445</td><td>0</td><td>0.00</td><td>thick green</td></tr><tr><td>0.133±0.021</td><td>15.27±0.14</td><td>10.90±0.10</td><td>0.715±0.024</td><td>9.54±0.26</td><td>0.012±0.014</td><td>0.542</td><td>2.80</td><td>-0.53</td><td>thin magenta</td></tr><tr><td>0.124±0.020</td><td>24.56±0.41</td><td>17.52±0.30</td><td>0.684±0.028</td><td>7.11±0.50</td><td>0.009±0.016</td><td>0.673</td><td>6.61</td><td>-1.85</td><td>dashed blue</td></tr><tr><td>0.107±0.021</td><td>12.76±0.14</td><td>9.11±0.10</td><td>0.724±0.031</td><td>10.14±0.28</td><td>0.018±0.015</td><td>1.079</td><td>18.38</td><td>-6.87</td><td>dotted black</td></tr><tr><td>0</td><td></td><td></td><td></td><td></td><td>0.013±0.014</td><td>2.025</td><td>45.82</td><td>-19.73</td><td></td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Top: Periodograms of the raw (top), filtered (middle) and residual (bottom) RV curves over the whole dataset (black line). The red line represents the 2015 Nov dataset, the green line the 2016 Jan dataset, the blue line the odd data points and the magenta line the even data points. FAP levels of 0.33, 0.10, 0.03 and 0.01 are displayed as horizontal dashed cyan lines. The rotation frequency (1.402 cycles per day) is marked by a cyan dashed line, as well as its first harmonic (2.803 cycles per day) and the frequency that has the smallest FAP (0.06% at 0.075 cycles per day, corresponding to $P_{orb}=13.41$ d). Aliases of the highest peaks, related to the observation window, appear as lower peaks separated by one cycle per day. Bottom: Zoom in the periodogram of filtered RVs.</figcaption></figure>

<table><caption>**Table 5.** Optimal orbital parameters derived with the method described in Sec. 5.2, respectively: semi-amplitude *K*, orbital period *P*<sub>orb</sub> in units of *P*<sub>rot</sub>, orbital period *P*<sub>orb</sub> in days, phase of inferior conjunction *φ* relative to cycle 11.0, BJD of inferior conjunction, χ<sup>2</sup><sub>r</sub>, Δχ<sup>2</sup> summed on 2581 data points, and natural logarithm of the likelihood *L*<sub>r2</sub> relative to the best fit. The case where no planet is taken into account in the model is given for comparison.</caption><thead><tr><th>K<br>(km s<sup>-1</sup>)</th><th>P<sub>orb</sub><br>(P<sub>rot</sub>)</th><th>P<sub>orb</sub><br>(d)</th><th>φ</th><th>BJD<sub>c</sub><br>(2,457,340+)</th><th>χ<sup>2</sup><sub>r</sub></th><th>Δχ<sup>2</sup></th><th>log L<sub>r2</sub></th></tr></thead><tbody><tr><td>0.154±0.022</td><td>15.29±0.15</td><td>10.91±0.11</td><td>0.671±0.035</td><td>9.06±0.38</td><td>0.96824</td><td>0.00</td><td>0.00</td></tr><tr><td>0.144±0.023</td><td>18.78±0.25</td><td>13.40±0.18</td><td>0.685±0.041</td><td>8.43±0.55</td><td>0.96979</td><td>4.00</td><td>-1.34</td></tr><tr><td>0.148±0.025</td><td>12.83±0.12</td><td>9.16±0.09</td><td>0.677±0.038</td><td>9.69±0.35</td><td>0.97180</td><td>9.17</td><td>-3.61</td></tr><tr><td>0</td><td></td><td></td><td></td><td></td><td>0.98631</td><td>46.62</td><td>-21.60</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>Figure 11. $\Delta\chi^2$ map as a function of K and $P_{orb}/P_{rot}$, derived with ZDI from corrected Stokes I LSD profiles at constant information content. Here the phase is fixed at 0.67, i.e., the value of $\phi$ at which the 3D paraboloid is minimum. The outer colour delimits the 99.99% confidence level area (corresponding to a $\chi^2$ increase of 21.10 for 2581 data points in our Stokes I LSD profiles). The minimum value of $\chi_t^2$ is 0.96824.</figcaption></figure>

By fitting the filtered RVs with a Keplerian orbit rather than a circular orbit, we obtain an eccentricity of 0.16± 0.15, indicating that there is no evidence for an eccentric orbit (following the precepts of Lucy & Sweeney 1971). We can thus conclude that the orbit of TAP 26 b is likely close to circular, or no more than moderately eccentric.

## 5.2 Deriving the planetary parameters from the LSD profiles

A second technique, following the method of Petit et al. (2015), consists of taking into account the presence of a planet into the ZDI model. Rather than fitting the measured Stokes I LSD profiles with a synthetic activity jitter directly, we first apply a translation in velocity to each of them, to remove the reflex motion caused by a planet of given parameters, and then apply ZDI to the corrected data set. Practically speaking, we repeat the experiment for a range of values for the orbital parameters ($K, P_{orb}, \phi$) at the vicinity of the minima previously identified in Sec. 5.1 and look for the set of values that yields the best result. The same way as for differential rotation, we derive the error bars on all parameters from the curvature of the 3D $\chi_t^2$ paraboloid around the minimum.

In the present case, since we have two datasets separated by a 45 d gap and we know that intrinsic variability occurred (see Sec. 4 and 5.1), a modification to the method described above was implemented: after correcting the global dataset from the reflex motion, ZDI is applied separately on each dataset, reconstructing two different brightness maps (one for late 2015 and one for early 2016) in order to obtain a more precise reconstruction. The quantity used to measure the likelihood of each set of parameters is therefore a global $\chi_t^2$, computed as a weighted average of both individual $\chi_t^2$,

with respective weights proportional to the number of data points in each set (1424 for 2015 Nov and 1157 for 2016 Jan).

As in the previous section, several minima are found, which are listed in Table 5. We also computed the relative likelihood of each case compared to the best one from the corresponding difference in $\chi_t^2$. We note that the case with no planet yields $\chi_t^2=0.98631$, which leads to a relative probability lower than $10^{-9}$ compared to the case with a 10.91 d period planet.

Figure 11 displays a $\Delta\chi^2$ map around the local minimum $P_{orb}/P_{rot}=15.29$, at $\phi=0.67$, showing the 99.99% confidence area.

## 5.3 Gaussian-Process Regression (GPR)

The third method we used works directly with the raw RVs and aims at modelling the activity jitter and its temporal evolution with GPR, assuming it obeys a priori covariance function (Haywood et al. 2014; Rajpaul et al. 2015). Similarly to the previous method, we fit both the orbit model and the jitter model simultaneously. For a planet with given parameters, we first remove the planet reflex motion from the RVs, then we fit the corrected RVs with a Gaussian process (GP) of pseudo-periodic covariance function:

$$c(t, t') = \theta_1^2 \exp \left[ -\frac{(t-t')^2}{\theta_3^2} - \frac{\sin^2 \left( \frac{\pi(t-t')}{\theta_2} \right)}{\theta_4^2} \right] \quad (3)$$

where $t$ and $t'$ are two dates, $\theta_1$ is the amplitude (in km s$^{-1}$) of the GP, $\theta_2$ the recurrence timescale (in units of $P_{rot}$), $\theta_3$ the decay timescale (i.e., the typical spot lifetime in the present case, in units of $P_{rot}$) and $\theta_4$ a smoothing parameter (within [0,1]) setting the amount of high frequency structures that we allow the fit to include. From a given set of orbital parameters ($K, P_{orb}, \phi$) and of covariance function parameters ($\theta_1$ to $\theta_4$, called hyperparameters), we can derive the GP that best fits the corrected RVs (noted y below) as well as the log likelihood $\log \mathcal{L}$ of the corresponding set of parameters from:

$$2 \log \mathcal{L} = -n \log(2\pi) - \log|C + \Sigma| - y^T(C + \Sigma)^{-1}y \quad (4)$$

where $n$ is the number of data points (29 in our case), $C$ is the covariance matrix of all the observing epochs and $\Sigma$ is the diagonal variance matrix of the raw RVs.

Coupled with a Markov Chain Monte Carlo (MCMC) simulation to explore the parameter domain, this method generates samples from the posterior probability distributions for the hyperparameters of the noise model and the orbital parameters. From these we can determine the maximum-likelihood values of these parameters and their uncertainty ranges. After an initial run where all the parameters are free to vary, we fix $\theta_4$ and $\theta_3$ to their respective best values (0.50±0.09 and 180±60 $P_{rot}=128±43$ d) before carrying out the main MCMC run to find the best estimates of the 5 remaining parameters. We note that the best value found for the decay time is exactly equal to the differential rotation lap time within error bars, and to twice the total span of our data. This decay time corresponds to both the differential rotation lap time and the starspot coherence time, since these are the most influent phenomena on the periodicity of the activity jitter. Such a starspot coherence
<table><caption>**Table 6.** Sets of orbital parameters that allow to fit the corrected RV curve best, using a GP with a covariance function given in Eq. 4, derived from the MCMC run. Respectively: reflex motion RV semi-amplitude *K*, orbital period *P*<sub>orb</sub> in units of *P*<sub>rot*, orbital period *P*<sub>orb</sub> in days, phase of inferior conjunction *φ* relative to rotation cycle 11.00 (ephemeris defined in Eq. 1), BJD of inferior conjunction, natural logarithm of the marginal likelihood *L* and natural logarithm of the relative marginal likelihood *L*<sub>r3</sub> as compared to the best case. The case where no planet is taken into account in the model is given for comparison.</caption><thead><tr><th>K<br>(km s<sup>-1</sup>)</th><th>P<sub>orb</sub><br>(P<sub>rot</sub>)</th><th>P<sub>orb</sub><br>(d)</th><th>φ</th><th>BJD<sub>c</sub><br>(2,457,340+)</th><th>log L</th><th>log L<sub>r3</sub></th></tr></thead><tbody><tr><td>0.163<br>±0.028</td><td>12.61<br>±0.13</td><td>8.99<br>±0.09</td><td>0.766<br>±0.030</td><td>10.54<br>±0.27</td><td>-3.48</td><td>0.00</td></tr><tr><td>0.149<br>±0.026</td><td>15.12<br>±0.20</td><td>10.79<br>±0.14</td><td>0.728<br>±0.033</td><td>9.71<br>±0.36</td><td>-3.73</td><td>-0.25</td></tr><tr><td>0.139<br>±0.026</td><td>18.74<br>±0.34</td><td>13.37<br>±0.24</td><td>0.694<br>±0.042</td><td>8.56<br>±0.57</td><td>-5.60</td><td>-2.12</td></tr><tr><td>0</td><td></td><td></td><td></td><td></td><td>-15.80</td><td>-12.52</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Phase plots of our 5-parameter MCMC run with yellow, red and blue points marking respectively the 1σ, 2σ and 3σ confidence regions. The optimal values found for each parameters are: θ<sub>1</sub>=1.19±0.21 km s<sup>-1</sup>, θ<sub>2</sub>=1.0005±0.0002 P<sub>rot</sub>, K=0.152±0.029 km s<sup>-1</sup>. Several optima are detected for P<sub>orb</sub>: 12.61±0.13 P<sub>rot</sub>, 15.12±0.20 P<sub>rot</sub> and 18.74±0.34 P<sub>rot</sub>, ordered by decreasing likelihood. The corresponding phases φ are: 0.766±0.030, 0.728±0.033 and 0.694±0.042 respectively.</figcaption></figure>

time is consistent with previous studies (Grankin et al. 2008; Bradshaw & Hartigan 2014; Lanza 2006).

As shown in Fig. 12, this method successfully recovers the different minima previously found with the first two techniques, with little correlation between the various parameters thus minimum bias in the derived values. Applying the method of Chib & Jeliazkov (2001) to the MCMC posterior samples, we obtain that the marginal likelihood of the case P<sub>orb</sub>=12.61 P<sub>rot</sub> is larger than that of the case P<sub>orb</sub>=15.12 P<sub>rot</sub> by a Bayes factor of only 1.28, which implies that there is as

yet no clear evidence in favor of either of them. The third most likely case, P<sub>orb</sub>=18.74 P<sub>rot</sub>, has a marginal likelihood which is inferior to the first one by a Bayes' factor of > 8, and the case with no planet has a marginal likelihood which is smaller than that of the first case by a Bayes factor of 2 10<sup>5</sup>. The three most likely sets of parameters are summarised in Table 6.

Trying to fit a non-circular Keplerian orbit to our data, i.e. adding the periapsis argument and the eccentricity *e* to the parameters in our MCMC run, we obtain *e*=0.05±0.18,
<figure><img src="image_14.png" /><figcaption>**Figure 13.** RV curves for a GPR fit of the activity jitter, with parameters $K=0.163$ km s$^{-1}$, $P_{\text{orb}}=12.61$ $P_{\text{rot}}$, $\phi=0.766$, $\theta_1=1.19$ km s$^{-1}$, $\theta_2=1.0005$ $P_{\text{rot}}$, $\theta_3=180$ $P_{\text{rot}}$, $\theta_4=0.50$ $P_{\text{rot}}$. Top panel: raw RVs and their error bars are shown in red, the solid cyan curve is the sum of the activity jitter predicted by GPR and the planet signal, and the dashed cyan lines show the 68.3% confidence intervals about the prediction around this model. Middle panel: filtered RVs and their error bars, resulting from the subtraction of the GP-fitted activity jitter from the raw RVs (in red), and the sine curve corresponding to the assumed planet signal (in cyan). Bottom panel: residual RVs resulting from the subtraction of the planet signal from the filtered RVs, and their error bars. The residual RVs feature a rms value of 29 m s$^{-1}$, i.e. the GP fits the RVs down to $\chi^2=0.151$.</figcaption></figure>

with a marginal likelihood slightly smaller than that of the case of a circular orbit. This further supports that the planet eccentricity is low if non-zero.

The best fit with our third method is shown in Fig. 13, where we see the raw RVs and the modelled RV curve predicted with this method, i.e., the sum of the GPR-fitted activity jitter and of the planet signal. Zooming in shows that this curve presents similarities with the RV jitter curve derived by ZDI (Fig. 8), indicating that, although working only with the RV data points, GPR successfully retrieves a convincing model for the activity. We also note the ability of the GP to model the activity jitter not only during our observing runs, but also during the 45 d gap between them, emphasizing the variability of the RV signal with time. The residual RVs in the case presented here have a rms value of 29 m s$^{-1}$ (close to the instrument RV precision 20-30 m s$^{-1}$) whereas the residual RVs derived with the 1st method yield a rms value of 51 m s$^{-1}$. Though the rms value is 2.5 times smaller than the error bar, GPR only fits 2 parameters, which illustrates its flexibility without decreasing its reliability, since the results are consistent with those found using independent methods (Sec. 5.1, 5.2). This demonstrates that GPR

does a better job at modelling the activity jitter and its temporal evolution than the 2 previous methods, in agreement with the conclusions of Donati et al. (2017) in the case of the wTTS V830 Tau. As a result, we consider the optimal planet parameters derived with GPR as the most reliable ones, and therefore conclude that the orbital periods of 10.8 and 9.0 d are more or less equally likely.

Table 7 summarises the likelihood of the different periods found with each method.

# 6 SUMMARY & DISCUSSION

This paper reports the results of an extended spectropolarimetric run on the wTTS TAP 26, carried out within the framework of the international MaTYSSE Large Programme, using the echelle spectropolarimeter ESPaDOnS at CFHT, spanning 72 d from 2015 Nov 18 to 2015 Dec 03 then from 2016 Jan 17 2016 Jan 29, and complemented by contemporaneous photometric observations from the 1.25 m telescope at CrAO.

Applying Zeeman-Doppler Imaging (ZDI) to our two
<table><caption>**Table 7.** Results yielded by the methods ZDI #1 (Sec. 5.1), ZDI #2 (Sec. 5.2) and GPR (Sec. 5.3), for the two periods ≈15 P<sub>rot</sub> and ≈13 P<sub>rot</sub>. From top to bottom: reflex motion semi-amplitude *K*, phase of inferior conjunction *φ* relative to cycle 11.0, orbital period *P<sub>orb</sub>* in units of *P<sub>rot</sub>*; orbital period *P<sub>orb</sub>* in days, semi-major axis *a*, M sin i, BJD of inferior conjunction BJD<sub>c</sub>, natural logarithm of relative likelihood as compared to the best case *L<sub>r</sub>*; GP amplitude θ<sub>1</sub> and GP recurrence timescale θ<sub>2</sub>. Results are displayed in bold font when the period is found with the highest likelihood using the corresponding method.</caption><thead><tr><th></th><th>ZDI #1</th><th>ZDI #2</th><th>GPR</th></tr></thead><tbody><tr><td>K (km s<sup>-1</sup>)</td><td>0.133±0.021</td><td><strong>0.154±0.022</strong></td><td>0.149±0.026</td></tr><tr><td>φ</td><td>0.715±0.024</td><td><strong>0.671±0.035</strong></td><td>0.728±0.033</td></tr><tr><td>P<sub>orb</sub> (P<sub>rot</sub>)</td><td>15.27±0.14</td><td><strong>15.29±0.15</strong></td><td>15.12±0.20</td></tr><tr><td>P<sub>orb</sub> (d)</td><td>10.90±0.10</td><td><strong>10.91±0.11</strong></td><td>10.79±0.14</td></tr><tr><td>a (au)</td><td>0.0974±0.0032</td><td><strong>0.0975±0.0032</strong></td><td>0.0968±0.0032</td></tr><tr><td>M sin i (M<sub>Jup</sub>)</td><td>1.49±0.25</td><td><strong>1.73±0.27</strong></td><td>1.66±0.31</td></tr><tr><td>BJD<sub>c</sub> (2,457,340+)</td><td>9.54±0.26</td><td><strong>9.06±0.38</strong></td><td>9.71±0.36</td></tr><tr><td>log L<sub>r</sub></td><td>-0.53</td><td><strong>0.00</strong></td><td>-0.25</td></tr><tr><td>θ<sub>1</sub> (km s<sup>-1</sup>)</td><td></td><td></td><td>1.19±0.21</td></tr><tr><td>θ<sub>2</sub> (P<sub>rot</sub>)</td><td></td><td></td><td>1.0004±0.0002</td></tr><tr><td>K (km s<sup>-1</sup>)</td><td>0.107±0.021</td><td>0.148±0.025</td><td><strong>0.163±0.028</strong></td></tr><tr><td>φ</td><td>0.724±0.031</td><td>0.677±0.038</td><td><strong>0.766±0.030</strong></td></tr><tr><td>P<sub>orb</sub> (P<sub>rot</sub>)</td><td>12.76±0.14</td><td>12.83±0.12</td><td><strong>12.61±0.13</strong></td></tr><tr><td>P<sub>orb</sub> (d)</td><td>9.11±0.10</td><td>9.16±0.09</td><td><strong>8.99±0.09</strong></td></tr><tr><td>a (au)</td><td>0.0864±0.0028</td><td>0.0868±0.0028</td><td><strong>0.0858±0.0028</strong></td></tr><tr><td>M sin i (M<sub>Jup</sub>)</td><td>1.13±0.23</td><td>1.56±0.28</td><td><strong>1.71±0.31</strong></td></tr><tr><td>BJD<sub>c</sub> (2,457,340+)</td><td>10.14±0.28</td><td>9.69±0.35</td><td><strong>10.54±0.27</strong></td></tr><tr><td>log L<sub>r</sub></td><td>-6.87</td><td>-3.61</td><td><strong>0.00</strong></td></tr><tr><td>θ<sub>1</sub> (km s<sup>-1</sup>)</td><td></td><td></td><td>1.19±0.21</td></tr><tr><td>θ<sub>2</sub> (P<sub>rot</sub>)</td><td></td><td></td><td>1.0005±0.0002</td></tr></tbody></table>

data sets, we derived the surface brightness and magnetic maps of TAP 26, revealing the presence of cool spots and warm plages totalling up to 12% of the stellar surface (we however caution that this is a lower limit given the insensitivity of ZDI to small spots evenly spread over the stellar surface). The large-scale field of TAP 26 is found to be mainly poloidal and axisymmetric, with a 120 G dipole component tilted at 40° from the rotation axis. The 2015 Nov and 2016 Jan maps are mostly similar, but nonetheless feature some differences that indicate temporal evolution of the surface brightness and the magnetic field, demonstrated by the inability of ZDI to model the whole dataset at noise level, on a timescale comparable to that spanning our sample (72 d). ZDI also enabled us to detect the differential rotation pattern at the surface of TAP 26, with $dΩ=0.0492±0.0010$ rad d$^{-1}$, a value close to that of the Sun, implying a time for the equator to lap the pole by one rotation equal to 128±3 d.

We then applied three different methods to search for a planetary signature in the observed spectra. The first method studies the radial velocities filtered out from the activity jitter predicted by ZDI. Our second method looks for the planet parameters that enable the best fit to the corrected LSD profiles, in a way similar to that used to estimate surface differential rotation. The third method uses Gaussian-Process Regression (GPR) to fit the activity jitter in the raw RVs, and like the second method, searches for the orbital parameters which enable GPR to fit the raw RVs corrected from the reflex motion best. We find that GPR succeeds best at modelling the intrinsic variability occurring at the surface of TAP 26, and is able to fit raw RVs at a rms precision of 29 m s$^{-1}$, i.e., close to the instrumental precision of ESPaDOnS (20-30 m s$^{-1}$, Moutou et al. 2007; Donati et al.

2008) and 30% better than with our first method (yielding a rms precision of 51 m s$^{-1}$). A similarly low rms was reached by GPR in the study of wTTS V830 Tau (35-37 m s$^{-1}$, Donati et al. 2017).

All three methods demonstrate the clear presence of a planet signature in the data, although the gap between both data sets generates aliasing problems, causing multiple nearby peaks to stand out in the periodogram. Of the dominant periods, the 10.8 d one emerges strongly for all three methods. It is the most likely with the second method, and equally likely as other periods when using the first and third methods (13.4 d and 9.0 d respectively). Although the 9.0 d orbital period ranks low (and in particular lower than the 13.4 d period) with our first and second methods, we nonetheless consider it as the second most likely given its first rank with GPR; the most probable explanation for this apparent discrepancy lies in the higher ability of GPR at modelling intrinsic variability of the activity jitter plaguing the RV curve. Allowing ZDI to model temporal evolution of spot distributions and magnetic topologies should bring all methods on an equal footing; this upgrade is planned for a forthcoming study.

Assuming the 10.79±0.14 d period is the true orbital period, and using the values yielded by GPR for *K* and *φ*, we find a circular orbit of semi-major axis *a* = 0.0968 ± 0.0032 au = 17.8 ± 2.7 R<sub>⋆</sub>, epoch of inferior conjunction BJD<sub>c</sub>=2,457,349.71±0.36 and M sin i=1.66±0.31 M<sub>Jup</sub>. If the orbital plane is aligned with the equatorial plane of TAP 26, with an assumed inclination of 55°, we obtain a mass M=2.03±0.46 M<sub>Jup</sub> for TAP 26 b. The 8.99±0.09 d period leads to *a*=0.086±0.003 au, BJD<sub>c</sub>=2,457,350.54±0.27 and M sin i=1.71±0.31 M<sub>Jup</sub>.

With an age of ≈17 Myr, TAP 26 is already an aging
T Tauri star and on the verge of becoming a post T Tauri star, as demonstrated by its complex geometry and weaker dipole field component (consistent with TAP 26 having a mostly radiative interior). Akin to V830 Tau b (Donati et al. 2016), the hJ in a nearly circular orbit that we have discovered in the young system TAP 26 is better explained by type II disc migration than by planet-planet scattering coupled to tidal circularisation. When compared to V830 Tau, a 2 Myr wTTS of similar mass (Donati et al. 2015, 2016, 2017), appears as an evolved version, rotating 4x faster than its younger sister, likely as a direct consequence of its 4x smaller moment of inertia (according to the evolutionary models of Siess et al. 2000).

Regarding the hJs we detected around TAP 26 and V830 Tau and despite their differences (in mass in particular), it would be tempting to claim that, like its host star, TAP 26 b is an evolved version of V830 Tau b. This would actually imply that TAP 26 b migrated outwards under tidal forces from a distance of ≈0.057 au (where V830 Tau b is located) to its current orbital distance of 0.094 au, as a result of the spin period of TAP 26 being ≈15x shorter than the orbital period of TAP 26 b. This option seems however unlikely given the latest predictions of tidal interactions between a young T Tauri star and its close-in hot Jupiter (Bolmont & Mathis 2016), indicating that tidal forces can only have a significant impact on a hJ within 0.06 au of a solar-mass host star (for a typical TTS with a radius of ≈2 R☉). The most likely explanation we see is thus that TAP 26 b:

* ended up its type-II migration in the accretion disc at the current orbital distance, when TAP 26 was still young, fully convective and hosting a large-scale dipole field of a few kG similar to that of AA Tau (Donati et al. 2010), i.e., strong enough to disrupt the disc up to a distance of 0.09 au.

* was left over once the disc has dissipated at an age significantly smaller than 2 Myr, i.e., before the large-scale field had time to evolve into a weaker and more complex topology, and the inner accretion disc to creep in as a result of the decreasing large-scale field and the subsequent chaotic accretion (e.g., Blinova et al. 2016).

Admittedly, this scenario requires favorable conditions to operate; in particular, it needs the accretion disc to vanish in less than 2 Myr, which happens to occur in no more than 10% of single T Tauri stars in Taurus (Kraus et al. 2012). In fact, since both TAP 26 and V830 Tau have the same angular momentum content, it is quite likely that TAP 26 indeed dissipated its disc very early (see Sec. 3). Quantitatively speaking, assuming (i) that the hJ we detected tracks the location of the inner disc when the disc dissipated, (ii) that the spin period at this time was locked on the Kepplerian period of the inner disc (equal to the orbital period of the detected hJ) and (iii) that stellar angular momentum was conserved since then, we derive that the disc must have dissipated when TAP 26 was about three times larger in radius, at an age of less than 1 Myr (according to Siess et al. 2000). Generating a magnetospheric cavity of the adequate size (0.085 to 0.097 au depending on the orbital period) would have required TAP 26 to host at this time a large scale dipole field of 0.3-1.0 kG for mass accretion rates in the range $10^{-9}$ to $10^{-8}$ M☉/yr, compatible with the large-scale fields found in cTTSs of similar masses (e.g., GQ Lup, Donati et al. 2012).

Along with other recent reports of close-in giant planets (or planet candidates) detected (or claimed) around young stars (van Eyken et al. 2012; Donati et al. 2016, 2017; Mann et al. 2016; Johns-Krull et al. 2016; David et al. 2016), our result may suggest a surprisingly high frequency of hJs around young solar-type stars, with respect to that around more evolved stars (≈1%, Wright et al. 2012). However, this may actually reflect no more than a selection bias in the observation samples (as for their mature equivalents in the early times of velocimetric planet detections). Planets are obviously much easier to detect around non-accreting TTSs as a result of their lower level of intrinsic variability; observation samples (like that of MaTYSSE) are thus naturally driven towards young TTSs whose accretion discs vanished early, i.e., at a time when their large-scale fields were still strong and their magnetospheric gaps large, and thus for which hJs had more chances to survive type-II migration. A more definite conclusion must wait for a complete analysis of the full MaTYSSE sample.

More observations of TAP 26, featuring in particular a more regular temporal sampling, are currently being planned to better determine the characteristics of the new-born hJ we detected. Furthermore, analysing thoroughly the full MaTYSSE data set to pin down the frequency of new-born hJs within the sample observed so far will bring a clearer view on how the formation and migration of young giant planets is occurring. Ultimately, only a full-scale planet survey of young TTSs such as that to be carried out with SPIRou, the new generation spectropolarimeter currently being built for CFHT and scheduled for first light in 2018, will be able to bring a consistent picture of how young close-in planets form and migrate, how their population relates to that of mature hJs, and more generally how young hJs impact the formation and early architecture of planetary systems like our Solar System.

## ACKNOWLEDGEMENTS

This paper is based on observations obtained at the Canada-France-Hawaii Telescope (CFHT), operated by the National Research Council of Canada, the Institut National des Sciences de l'Université of the Centre National de la Recherche Scientifique (INSU/CNRS) of France and the University of Hawaii. We thank the CFHT QSO team for the great work and effort at collecting the high-quality MaTYSSE data presented in this paper. MaTYSSE is an international collaborative research programme involving experts from more than 10 different countries (France, Canada, Brazil, Taiwan, UK, Russia, Chile, USA, Ireland, Switzerland, Portugal, China and Italy). We also warmly thank the IDEX initiative at Université Fédérale Toulouse Midi-Pyrénées (UFT-MiP) for funding the STEPS collaboration program between IRAP/OMP and ESO. We acknowledge funding from the LabEx OSUG@2020 that allowed purchasing the ProLine PL230 CCD imaging system installed on the 1.25-m telescope at CrAO. SGG acknowledges support from the Science & Technology Facilities Council (STFC) via an Ernest Rutherford Fellowship [ST/J003255/1]. SHPA acknowledges financial support from CNPq, CAPES and Fapemig.
REFERENCES

Baruteau C., et al., 2014, Protostars and Planets VI, pp 667–689

Batygin K., Bodenheimer P. H., Laughlin G. P., 2016, ApJ, 829, 114

Blinova A. A., Romanova M. M., Lovelace R. V. E., 2016, MNRAS, 459, 2354

Bolmont E., Mathis S., 2016, Celestial Mechanics and Dynamical Astronomy, 126, 275

Bradshaw S. J., Hartigan P., 2014, ApJ, 795, 79

Brown S. F., Donati J.-F., Rees D. E., Semel M., 1991, A&A, 250, 463

Chib S., Jeliazkov I., 2001, Journal of the American Statistical Association, 96, 270

David T. J., et al., 2016, Nature, 534, 658

Davies C. L., Gregory S. G., Greaves J. S., 2014, MNRAS, 444, 1157

Donati J.-F., 2003, in Trujillo-Bueno J., Sanchez Almeida J., eds, Astronomical Society of the Pacific Conference Series Vol. 307, Astronomical Society of the Pacific Conference Series. pp 41–+

Donati J.-F., Brown S. F., 1997, A&A, 326, 1135

Donati J.-F., Collier Cameron A., 1997, MNRAS, 291, 1

Donati J.-F., Semel M., Carter B. D., Rees D. E., Collier Cameron A., 1997, MNRAS, 291, 658

Donati J.-F., Collier Cameron A., Petit P., 2003, MNRAS, 345, 1187

Donati J.-F., et al., 2006, MNRAS, 370, 629

Donati J.-F., et al., 2008, MNRAS, 385, 1179

Donati J., et al., 2010, MNRAS, 409, 1347

Donati J., et al., 2011, MNRAS, 412, 2454

Donati J.-F., et al., 2012, MNRAS, 425, 2948

Donati J.-F., et al., 2014, MNRAS, 444, 3220

Donati J.-F., et al., 2015, MNRAS, 453, 3706

Donati J. F., et al., 2016, Nature, 534, 662

Donati J.-F., et al., 2017, MNRAS, 465, 3343

Feigelson E. D., Jackson J. M., Mathieu R. D., Myers P. C., Walter F. M., 1987, AJ, 94, 1251

Gallet F., Bouvier J., 2015, A&A, 577, A98

Grankin K. N., 2013, Astronomy Letters, 39, 251

Grankin K. N., Bouvier J., Herbst W., Melnikov S. Y., 2008, A&A, 479, 827

Gregory S. G., Donati J.-F., Morin J., Hussain G. A. J., Mayne N. J., Hillenbrand L. A., Jardine M., 2012, ApJ, 755, 97

Haywood R. D., et al., 2014, MNRAS, 443, 2517

Herczeg G. J., Hillenbrand L. A., 2014, ApJ, 786, 97

Huélamo N., et al., 2008, A&A, 489, L9

Johns-Krull C. M., et al., 2016, ApJ, 826, 206

Kenyon S. J., Hartmann L., 1995, ApJS, 101, 117

Kraus A. L., Ireland M. J., Hillenbrand L. A., Martinache F., 2012, ApJ, 745, 19

Kurucz R., 1993, CDROM # 13 (ATLAS9 atmospheric models) and # 18 (ATLAS9 and SYNTHE routines, spectral line database). Smithsonian Astrophysical Observatory, Washington D.C.

Lagrange A.-M., et al., 2010, Science, 329, 57

Landi degl’Innocenti E., Landolfi M., 2004, Polarisation in spectral lines. Dordrecht/Boston/London: Kluwer Academic Publishers

Lanza A. F., 2006, MNRAS, 369, 1773

Loinard L., Torres R. M., Mioduszewski A. J., Rodríguez L. F., González-Lópezlira R. A., Lachaume R., Vázquez V., González E., 2007, ApJ, 671, 546

Lucy L. B., Sweeney M. A., 1971, AJ, 76, 544

Mann A. W., et al., 2016, AJ, 152, 61

Moutou C., et al., 2007, A&A, 473, 651

Pecaut M. J., Mamajek E. E., 2013, ApJS, 208, 9

Pecaut M. J., Mamajek E. E., 2016, MNRAS, 461, 794

Petit P., et al., 2015, A&A, 584, A84

Rajpaul V., Aigrain S., Osborne M. A., Reece S., Roberts S., 2015, MNRAS, 452, 2269

Sallum S., et al., 2015, Nature, 527, 342

Semel M., 1989, A&A, 225, 456

Setiawan J., Henning T., Launhardt R., Müller A., Weise P., Kürster M., 2008, Nature, 451, 38

Siess L., Dufour E., Forestini M., 2000, A&A, 358, 593

Skelly M. B., Donati J.-F., Bouvier J., Grankin K. N., Unruh Y. C., Artemenko S. A., Petrov P., 2010, MNRAS, 403, 159

Stout-Batalha N. M., Vogt S. S., 1999, ApJS, 123, 251

Torres R. M., Loinard L., Mioduszewski A. J., Rodríguez L. F., 2009, ApJ, 698, 242

Valenti J. A., Fischer D. A., 2005, ApJS, 159, 141

Vogt S. S., Penrod G. D., Hatzes A. P., 1987, ApJ, 321, 496

Wright J. T., Marcy G. W., Howard A. W., Johnson J. A., Morton T. D., Fischer D. A., 2012, ApJ, 753, 160

van Eyken J. C., et al., 2012, ApJ, 755, 42

APPENDIX A: ADDITIONAL FIGURES

Images of brightness and magnetic field on the surface of TAP 26, as derived with ZDI using our 29 spectra, are shown in Fig. A1.

APPENDIX B: ACTIVITY PROXIES

In order to investigate whether the detected periodic RV signal may relate to activity, we plotted periodograms of the longitudinal magnetic field $B_l$ and of the $H\alpha$ emission equivalent width (Fig. B1 and B2 respectively). Peak frequencies for these proxies are located at periods of 0.7145±0.0002 d and 0.7132±0.0001 d respectively, as well as their aliases. Given the surface differential rotation parameters measured for TAP 26 (see Sec. 4.2), the values of their respective periods indicate that the longitudinal field traces an average latitude of 46° whereas the bulk of $H\alpha$ emission comes from a lower average latitude of 27° (see Eq. 2). As opposed to the raw RVs, the rotation period $P_{\text{rot}}$ has a higher power than its first harmonic $P_{\text{rot}}/2$ (Fig. 10). No signal is detected at the planet periods found in Sec. 5.

Plotting phase-folded curves of the longitudinal magnetic field and the $H\alpha$ emission equivalent width (where the x-axis indicates the rotation phase as defined in Eq. 1), in Fig. B3 and B4, we observe a decrease in the longitudinal magnetic field around phase 0.77 in 2015 Nov and phase 0.97 in 2016 Jan, which correspond approximately to the phases where the dipole pole points towards the Earth (0.73±0.03 and 0.85±0.03 respectively), causing $B_l$ to have strong negative values and showing the importance of the dipole in the value of $B_l$. Similarly, the increase in emission equivalent width of the $H\alpha$ line between phases 0.6 and 0.9 illustrates the correlation between the lower harmonics of the magnetic field of TAP 26 and this activity proxy.

<footnote n="1">¹ Université de Toulouse, UPS-OMP, IRAP, 14 avenue E. Belin, Toulouse, F-31400 France</footnote>

<footnote n="2">² CNRS, IRAP / UMR 5277, Toulouse, 14 avenue E. Belin, F-31400 France</footnote>

<footnote n="3">³ Department of Physics and Astronomy, York University,</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure A1.** Brightness and magnetic components surface maps when fitting the 2015 Nov and 2016 Jan datasets altogether, at rotation cycle 51.</figcaption></figure>

Toronto, Ontario L3T 3R1, Canada

⁴ CFHT Corporation, 65-1238 Mamalahoa Hwy, Kamuela, Hawaii 96743, USA

⁵ Département de physique, Université de Montréal, C. P. 6128, Succursale Centre-Ville, Montréal, QC, Canada H3C 3J7

⁶ Crimean Astrophysical Observatory, Nauchny, Crimea 298409

⁷ ESO, Karl-Schwarzschild-Str 2, D-85748 Garching, Germany

⁸ SUPA, School of Physics & Astronomy, Univ. of St Andrews, St Andrews, Scotland KY16 9SS, UK

⁹ School of Physics, Trinity College Dublin, the University of Dublin, Ireland

¹⁰ Departamento de Física – ICEx – UFMG, Av. Antônio Carlos, 6627, 30270-901 Belo Horizonte, MG, Brazil

¹¹ Université Grenoble Alpes, IPAG, BP 53, F-38041 Grenoble Cédex 09, France

¹² CNRS, IPAG / UMR 5274, BP 53, F-38041 Grenoble Cédex 09, France

¹³ Institute of Astronomy and Astrophysics, Academia Sinica, PO Box 23-141, 106, Taipei, Taiwan

¹⁴ Kavli Institute for Astronomy and Astrophysics, Peking University, Yi He Yuan Lu 5, Haidian Qu, Beijing 100871, China

¹⁵ LUPM, Université de Montpellier, CNRS, place E. Bataillon, F-34095, Montpellier, France
<figure><img src="image_16.png" /><figcaption>**Figure B1.** Periodogram of the longitudinal magnetic field. The rotation period at 0.7135 d is represented by a dashed vertical cyan line, as well as the orbital period at 10.92 d.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure B2.** Periodogram of the Hα line equivalent width. The rotation period at 0.7135 d is represented by a dashed vertical cyan line, as well as the orbital period at 10.92 d.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure B3.** Folded curve of the longitudinal magnetic field against the rotation phase. 2015 Nov (red upward-pointing triangles) data are fitted with the sum of a sine curve and 1 harmonic (red dashed line) and 2016 Jan (blue downward-pointing triangles) data are fitted with the sum of a sine curve and 2 harmonics (blue dotted line).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure B4.** Folded curve of the equivalent width of Hα against the rotation phase. 2015 Nov (red upward-pointing triangles) and 2016 Jan (in blue) data are fitted with the sum of a sine curve and 2 harmonics (red dashed line and blue dotted line respectively).</figcaption></figure>
APPENDIX C: SIMULATIONS

In order to check that the detected planetary signal does not come from the filtering process, we conducted simulations to test our three methods on two different data sets: one where the presence of a planet was input in the simulation (scenario #1), and one without any planet (scenario #2). Stokes I and V LSD profiles were generated from the brightness and magnetic maps found in Sec. 4, at the same dates of observation as the real data, with a comparable noise level. The added planet signature had the properties of the best fit found with the second method (Sec. 5.2): $K=0.154 \text{ km s}^{-1}$, $P_{\text{orb}}/P_{\text{rot}}=15.29$, $\phi=0.671$. Applying ZDI to these data sets, we reconstructed brightness and magnetic maps as in Sec. 4. For both simulations, the maps we found look similar to the ones reconstructed from the real data, with an information loss amounting to 4% for the spottedness and ≈80 G for the rms magnetic flux, but the main features, such as the polar spot, are recovered. Fig. C1 shows the brightness maps for simulation #1, at both epochs.

As in Sec. 5.1, synthetic RV curves are shown in Fig. C2, C3 for simulations #1 and #2 respectively. While a signal is detected in the filtered RVs of simulation #1 (rms 107 m s⁻¹), no significant signal is detected in the filtered RVs in simulation #2 (rms 58 m s⁻¹). Table C1 summarises the characteristics of the best fit to the filtered RVs for both scenarios, in comparison with a ($K=0 \text{ km s}^{-1}$, $RV_0=0 \text{ km s}^{-1}$) curve. The periodograms of the filtered RV curves, displayed in Fig. C4 and Fig. C5, further confirm this, with simulation #2 yielding no significant signal at the frequencies found with our three methods in Sec. 5.

We note that changing the noise pattern can make the FAP of the filtered RVs highest peak vary between 4% and

<table><caption>**Table C1.** Results found with the 1st method on both simulation datasets. The first column shows the results on the scenario with a planet and the second column shows the results on the scenario without planet. For each, a comparison is made between the best sine fit to the filtered RVs and a fit by a constant value, with the reflex motion semi-amplitude *K*, the orbital period *P*<sub>orb</sub> in units of *P*<sub>rot</sub>, the phase of inferior conjunction *φ* relative to cycle 11.0, the mean RV *RV*<sub>0</sub> and *χ*<sup>2</sup><sub>r</sub>. Differences in *χ*<sup>2</sup> (summed on 29 data points) and in logarithmic (log<sub>e</sub>) likelihood are given in the last row.</caption><thead><tr><th>Scenario #1</th><th>Scenario #2</th></tr></thead><tbody><tr><td><i>K</i>=0.122±0.020 km s<sup>-1</sup><br><i>P</i><sub>orb</sub>/<i>P</i><sub>rot</sub>=15.35±0.16<br><i>φ</i>=0.647±0.026<br><i>RV</i><sub>0</sub>=0.018±0.014 km s<sup>-1</sup><br><i>χ</i><sup>2</sup><sub>r</sub>=0.540<br><i>K</i>=0 km s<sup>-1</sup><br><i>χ</i><sup>2</sup><sub>r</sub>=1.893<br><i>Δχ</i><sup>2</sup>=39.2<br><i>Δ(log L<sub>r</sub>)</i>=-16.57</td><td><i>K</i>=0.036±0.021 km s<sup>-1</sup><br><i>χ</i><sup>2</sup><sub>r</sub>=0.436<br><i>K</i>=0 km s<sup>-1</sup><br><i>χ</i><sup>2</sup><sub>r</sub>=0.529<br><i>Δχ</i><sup>2</sup>=2.7<br><i>Δ(log L<sub>r</sub>)</i>=-0.49</td></tr></tbody></table>

<table><caption>**Table C2.** Characteristics of the best fit found with the 2nd method on simulation #1 (top row), compared to a fit with a no-planet model (i.e. $K=0 \text{ km s}^{-1}$, middle row). Differences in $\chi^2$ (summed on 2581 data points) and in logarithmic (log<sub>e</sub>) likelihood are given in the last row.</caption><thead><tr><th>Scenario #1</th></tr></thead><tbody><tr><td><i>K</i>=0.155±0.022 km s<sup>-1</sup><br><i>P</i><sub>orb</sub>/<i>P</i><sub>rot</sub>=15.32±0.14<br><i>φ</i>=0.662±0.036<br><i>χ</i><sup>2</sup><sub>r</sub>=0.95226</td></tr><tr><td><i>K</i>=0 km s<sup>-1</sup><br><i>χ</i><sup>2</sup><sub>r</sub>=0.97529</td></tr><tr><td><i>Δχ</i><sup>2</sup>=30.0<br><i>Δ(log L<sub>r</sub>)</i>=-10</td></tr></tbody></table>

26% for simulation #2 (no planet). It can also change the relative power of the different orbital periods in the filtered RVs for simulation #1: most of the time the ≈10.8 d period is recovered as the highest peak with a FAP<0.5%, but the ≈9.0 d period reaches a smaller $\chi^2$ in one case out of five. This sheds light on why the different methods do not always favour the same orbital periods in our analysis.

The second method also recovers the different orbital periods from simulation #1, the ≈10.8 d one being the most likely, with a $\Delta\chi^2$ of 7.37 compared to the ≈13.4 d period and 7.67 compared to the ≈9.0 d period. Fig. C6 shows the $\chi^2$ map around the minimum $P_{\text{orb}}$≈10.8 d for $\phi=0.67$ (value at the 3-D local minimum), with the white colour bounding the 99.99% confidence region. We chose not to apply this method on simulation #2 (without planet) because it runs computations for orbital parameters close to the local minima found with the first method, and no such significant minima were found for simulation #2.

GPR also successfully recovers the input planet period in simulation #1, as shown in Fig. C7 on which the raw RVs, the planet signal and the residual RVs at $P_{\text{orb}}=15.31P_{\text{rot}}$
<table><caption>**Table C3.** Characteristics of the best fit (first row) found with the 3rd method on simulation #1 (left) and #2 (right), compared to a fit with a no-planet model (i.e. K=0 km s⁻¹, middle row). Differences in logarithmic (log_e) likelihood are given in the last row.</caption><thead><tr><th>Scenario #1</th><th>Scenario #2</th></tr></thead><tbody><tr><td>K=0.138±0.027 km s<sup>-1</sup><br>P<sub>orb</sub>/P<sub>rot</sub>=15.31±0.21<br>φ=0.646±0.038<br>θ<sub>1</sub>=1.14±0.21 km s<sup>-1</sup><br>θ<sub>2</sub>=1.0002±0.0002 P<sub>rot</sub><br>log L<sub>r</sub>=-6.20</td><td>K=0.060±0.052 km s<sup>-1</sup><br>log L<sub>r</sub>=-5.99</td></tr><tr><td>K=0 km s<sup>-1</sup><br>log L<sub>r</sub>=-21.42</td><td>K=0 km s<sup>-1</sup><br>log L<sub>r</sub>=-5.48</td></tr><tr><td>Δ(log L<sub>r</sub>)=-15.22</td><td>Δ(log L<sub>r</sub>)=-0.51</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>**Figure C1.** Brightness maps reconstructed from the simulation #1 data, for the 2015 Nov data subset (top) and the 2016 Jan data subset (bottom). Both maps feature a spot coverage of ≈8%.</figcaption></figure>

are seen. The results of the MCMC runs are shown in the phase plots of Fig. C8 and C9, demonstrating that no orbital period stands out in simulation #2, i.e., in the activity jitter synthesised from the reconstructed brightness maps, whereas several orbital periods are detected in simulation #1, ≈10.8 d and ≈9.0 d being respectively the most likely and the second most likely, with a Bayes factor of only 1.25 between them. The comparison with a MCMC

run where no planet is subtracted is shown in table C3 for both scenarios, demonstrating that, for scenario #1, taking a planet into account in the model results in a significant increase in the likelihood of the best fit, whereas it is not the case for scenario #2.

We conclude that all three methods enable us to recover the planet signal (scenario #1), and that the detected periods in the observational filtered RVs are not artifacts of the numerical process (scenario #2). Furthermore, for our particular observation window, the noise pattern can change the relative likelihood of the different detected peaks, as can the choice of the method to use.

This paper has been typeset from a TeX/LATeX file prepared by the author.
<figure><img src="image_21.png" /><figcaption>**Figure C2.** Simulation #1: raw, filtered and residual RV curves as derived with the method described in 5.1. The residual RVs feature a rms value of 59 m s⁻¹.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure C3.** Simulation #2: raw and filtered RV curves as derived with the method described in 5.1. The filtered RVs feature a rms value of 58 m s⁻¹.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure C4. Simulation #1: periodograms of the raw (top), filtered (middle) and residual (bottom) RVs.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure C5. Simulation #2: periodograms of the raw (top) and filtered (bottom) RVs.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure C6.** Simulation #1: $\Delta\chi^2$ map at $\phi=0.67$ as derived with the method described in 5.2, centered around the minimum $P_{\text{orb}}=15.3 P_{\text{rot}}$. Parameters values are found to be $K=0.155\pm0.022 \text{ km s}^{-1}$, $P_{\text{orb}}=15.32\pm0.14 P_{\text{rot}}$. The minimum value of $\chi^2$ is 0.95226.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure C7.** Example of GPR for simulation #1. Top: raw RVs (red dots) and GP fit+planet (cyan), middle: curves as derived with the method described in 5.3, for the local minimum $P_{\text{orb}}=15.31\pm0.21 P_{\text{rot}}$, $K=0.138\pm0.027$, $\phi=0.646\pm0.038$. The rms of the residual RVs amounts to $47 \text{ m s}^{-1}$.</figcaption></figure>
<figure><img src="image_27.png" /></figure>

**Figure C8.** Simulation #1: phase plots of the MCMC run as described in 5.3. We find $\theta_1=1.15\pm0.19$ km s$^{-1}$, $\theta_2=1.0002\pm0.0002$ $P_{\text{rot}}$, $K=0.140\pm0.026$ km s$^{-1}$, and the dominant orbital periods $P_{\text{orb}}=15.31\pm0.21$ $P_{\text{rot}}$ and $P_{\text{orb}}=12.74\pm0.13$ $P_{\text{rot}}$, with the corresponding phases being $\phi=0.646\pm0.038$ and $\phi=0.699\pm0.036$ respectively.

<figure><img src="image_28.png" /></figure>

**Figure C9.** Simulation #2: phase plots of the MCMC run as described in 5.3. We find $\theta_1=1.09\pm0.19$ km s$^{-1}$, $\theta_2=1.0003\pm0.0002$ $P_{\text{rot}}$ and $K=0.060\pm0.053$ km s$^{-1}$.