# Recovering planet radial velocity signals in the presence of starspot activity in fully convective stars

J. R. Barnes, S. V. Jeffers, G. Anglada-Escudé, C. A. Haswell, H. R. A. Jones, M. Tuomi, F. Feng, J. S. Jenkins, P. Petit

► To cite this version:

J. R. Barnes, S. V. Jeffers, G. Anglada-Escudé, C. A. Haswell, H. R. A. Jones, et al. Recovering planet radial velocity signals in the presence of starspot activity in fully convective stars. Monthly Notices of the Royal Astronomical Society, 2017, 466, pp.1733-1740. 10.1093/mnras/stw3170. insu-04836441

HAL Id: insu-04836441

https://insu.hal.science/insu-04836441v1

Submitted on 13 Dec 2024

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
<figure><img src="image_2.png" /></figure>

# Recovering planet radial velocity signals in the presence of starspot activity in fully convective stars

J. R. Barnes,<fnref n="1" />* S. V. Jeffers,<fnref n="2" />* G. Anglada-Escudé,<fnref n="3" />* C. A. Haswell,<fnref n="1" />

H. R. A. Jones,<fnref n="4" /> M. Tuomi,<fnref n="4" /> F. Feng,<fnref n="4" /> J. S. Jenkins<fnref n="5" /> and P. Petit<fnref n="6" />,<fnref n="7" />

Accepted 2016 December 2. Received 2016 December 1; in original form 2016 November 4

## ABSTRACT

Accounting for stellar activity is a crucial component of the search for ever-smaller planets orbiting stars of all spectral types. We use Doppler imaging methods to demonstrate that starspot-induced radial velocity variability can be effectively reduced for moderately rotating, fully convective stars. Using starspot distributions extrapolated from sunspot observations, we adopt typical M dwarf starspot distributions with low contrast spots to synthesize line profile distortions. The distortions are recovered using maximum entropy regularized fitting and the corresponding stellar radial velocities are measured. The procedure is demonstrated that for a late-M star harbouring an orbiting planet in the habitable zone. The technique is effective for stars with $v \sin i = 1–10\,\text{km s}^{-1}$, reducing the stellar noise contribution by factors of nearly an order of magnitude. With a carefully chosen observing strategy, the technique can be used to determine the stellar rotation period and is robust to uncertainties such as unknown stellar inclination. While demonstrated for late-type M stars, the procedure is applicable to all spectral types.

**Key words:** techniques: radial velocities – techniques: spectroscopic – planets and satellites: detection – stars: activity – stars: low-mass.

## 1 INTRODUCTION

Efforts aimed at determining planet occurrence rates are now focusing on low-mass planets orbiting the lowest mass stars. Instrumental precision, which is limited by stability, and intrinsic stellar activity are the main factors that determine the number of observations needed to reliably recover planetary signals. The most stable spectrometers amongst the current generation of dedicated radial velocity (RV) instruments regularly achieve ~1 $\text{ms}^{-1}$ (Mayor et al. 2003; Cosentino et al. 2012). However, obtaining this level of precision outside the standard optical range (0.4–0.7 µm), at red-optical (0.6–1.0 µm) and near-infrared (>1 µm) wavelengths is also required if the lowest mass, fully convective stars, with $M_* < 0.35 M_\odot$ (spectral type M4V or later), are to be surveyed efficiently. Barnes et al. (2014) demonstrated precision down to ~2.5 $\text{ms}^{-1}$ is possible with current technology operating at red-optical wavelengths,

and Gao et al. (2016) have demonstrated similar precision down to ~2 $\text{ms}^{-1}$ in the near-infrared K band. There are also dedicated RV surveys, including the Habitable Zone Planet Finder (Mahadevan et al. 2014) and SpectroPolarimètre Infra-Rouge (SPIRou; Thibault et al. 2012) that will target the lowest mass stars at the bottom of the main sequence with 1 $\text{ms}^{-1}$ instrumental RV precision. For M dwarfs, although the red-optical contains more Doppler information in the thousands of available absorption lines, enabling greater RV precision to be achieved (Reiners et al. 2010), the spectral energy distribution of low-mass stars peaks at near-infrared wavelengths. The Calar Alto high-Resolution search for M dwarfs with Exoearths with Near-infrared and optical Échelle Spectrographs survey (CARMENES; Quirrenbach et al. 2010) is addressing this issue by covering both spectral regions simultaneously.

Despite the expectation that fully convective stars cannot possess a solar-like dynamo in the absence of a convective-radiation boundary where strong shearing is believed to take place, M stars are still observed to possess strong magnetic fields. Moreover, by mid-M, field stars in the solar neighbourhood have not spun down and typically possess significant rotation (Jenkins et al. 2009); for

<footnote n="1">¹Department of Physical Sciences, The Open University, Walton Hall, Milton Keynes MK7 6AA, UK</footnote>

<footnote n="2">²Institut für Astrophysik, Georg-August-Universität, Friedrich-Hund-Platz 1, Friedrich-Hund-Platz 1, D-37077 Göttingen, Germany</footnote>

<footnote n="3">³School of Physics and Astronomy, Queen Mary, University of London, 327 Mile End Rd., London E1 4NS, UK</footnote>

<footnote n="4">⁴Centre for Astrophysics Research, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK</footnote>

<footnote n="5">⁵Departamento de Astronomía, Universidad de Chile, Camino del Observatorio, 1515 Las Condes, Santiago, Chile</footnote>

<footnote n="6">⁶Université de Toulouse, UPS-OMP, Institut de Recherche en Astrophysique et Planétologie, F-31000 Toulouse, France</footnote>

<footnote n="7">⁷CNRS, Institut de Recherche en Astrophysique et Planétologie, 14 Avenue Édouard Belin, F-31400 Toulouse, France</footnote>

<footnote>* E-mail: john.barnes@open.ac.uk (JRB); jeffers@astro.physik.uni-goettingen.de (SVJ); guillem.anglada@gmail.com (GA-E)</footnote>
an M5V star, the mean equatorial rotation, $v \sin i \sim 6$ km s$^{-1}$. Newton et al. (2016) have however suggested that $v \sin i$ upper limits measured from low resolution ($R < 40000$) spectra significantly overestimate rotation rates following a similar previous claim by Reiners, Joshi & Goldman (2012). By measuring rotation periods in a sample of nearly 400 M stars, Newton et al. find that the mid-low-mass M dwarfs tend to exhibit short rotation periods of $P < 10$ d or periods of $P > 70$ d, with a dearth of intermediate rotators. They suggest that this dichotomy exists because stars maintain fast rotation for a few Gyr before rapid spin down to slow rotation (100 d) by an age of 5 Gyr.

Fully convective stars show activity all the way down to spectral type M9V, despite an overall decline in activity (Mohanty & Basri 2003; Reiners & Basri 2010). There is a clear correlation between rotation period and Ca II chromospheric emission index, $\log R'_{HK}$, down to mid-M (Suárez Mascareño et al. 2015), while for 15 M5-M9 dwarfs, Barnes et al. (2014) found a correlation between r.m.s. RV and $H\alpha$ emission strength. Doppler imaging studies by Morin et al. (2008a), Phan-Bao et al. (2009) and Barnes et al. (2015) have shown that line profile variability in fully convective stars can be interpreted as starspots. This interpretation is valid, even amongst the latest spectral types (Barnes et al. 2015), where $\sim 100$ ms$^{-1}$ RV variation is observed for the rapidly rotating M9 dwarf, LP 944-20 ($v \sin i = 31$ km s$^{-1}$). Fully convective field M dwarfs are thus moderate rotators that display activity that causes significant RV variations. Crucially, monitoring only slow rotators amongst a population of moderate rotators will favour low-inclination (pole-on) stars for which planet detection is less favoured. If we wish to conduct unbiased RV surveys to determine planet occurrence rates for fully convective stellar hosts, we must find a means to account for the stellar activity component.

Modelling and removal of the stellar activity component in RV measurements, which can swamp planetary signals, can be achieved in a number of ways. Traditional techniques involving measuring the line bisector are most effective for simple cases where a single spot, or spot group, on the star induces the line asymmetries. Other methods have recently been developed that enable a more sophisticated approach that identifies the stellar activity signatures directly from the line profiles rather than the RV measurement. By fitting Gaussian profiles to residual time series spectra, Moulds et al. (2013) were able to reduce starspot-induced RVs by more than 80 per cent. The technique is applicable to stars with significant rotation: i.e. $v \sin i = 10–50$ km s$^{-1}$. Dumusque, Boisse & Santos (2014) have shown that RV jitter can be removed by modelling activity regions on the stellar surface, but this has only been applied to stars with spectral type similar to the Sun. Petit et al. (2015) have applied maximum-entropy techniques to simultaneously recover planetary RVs with spot signatures on stars. The method was shown to work well for $v \sin i$ values of $\ge 20$ km s$^{-1}$. Below this $v \sin i$, it is unclear whether Doppler imaging techniques can simultaneously recover the planet and spot distribution.

A large number of variables must be considered for the detection of a planetary RV signature in the presence of stellar activity. Planetary orbital elements include the orbital period, $P_{orb}$, and the stellar reflex motion velocity amplitude of the star, $K_*$. The stellar activity includes a potentially unknown distribution of spots with some characteristic photospheric-spot contrast that distort the line profiles in a quasi-periodic manner, associated with the stellar rotation period, $P_{rot}$. The observation strategy is important for recovery of the planet and modelling of the spots. Specifically good phase sampling on time-scales that minimizes activity evolution and our ability to obtain spectra with sufficient S/N ratios will determine the

effectiveness of our ability to model the line profiles. In this paper, we apply Doppler imaging techniques to show that the effects of cool starspots on line profiles can be modelled for stars with $v \sin i = 1–10$ km s$^{-1}$. We then show that this enables effective reduction of starspot-induced RVs. Reconstruction of the line profiles using maximum entropy regularization has the advantage that no prior assumption of the spot locations and sizes is required. In Section 2, we describe the spot models and line profile simulation and recovery. The technique is illustrated in Section 3 using a representative case. We also investigate the ability of the technique to recover parameters, including stellar rotation period and the importance of stellar axial inclination. Sensitivities for a broader set of cases are explored in Section 4 before a brief summary and discussion in Section 5.

## 2 TECHNIQUE

### 2.1 Spot models

In Barnes et al. (2011), we used extrapolated solar models to generate spot size distributions for greater than solar activity levels, following the work of Jeffers (2005). A key difference of our models compared with solar spot observations is that rather than confining spots to low latitudes, we allowed spots to be located at all latitudes and longitudes. Evidence for spots distributed more uniformly than FGK dwarf stars comes from previous findings from Doppler imaging studies of early M stars (Barnes & Collier Cameron 2001; Barnes, James & Cameron 2004) and more recent studies of stars at or below the fully convective boundary (Morin et al. 2008a; Phan-Bao et al. 2009; Barnes et al. 2015, 2016). The recovered spot distributions, along with the low contrasts required to model the spots, confirmed the low photometric amplitude variability observed during earlier monitoring of the latest M dwarf (Rockenfeller, Bailer-Jones & Mundt 2006) and in the much larger MEarth Project sample (Newton et al. 2016).

Fig. 1 shows the three starspot models that we use in this paper. These models are identical to those presented in Barnes et al. (2011).

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Spot Models 2–4 based on Jeffers (2005) and adapted to M dwarf simulations in Barnes, Jeffers & Jones (2011). Model 2 is analogous to solar max activity level while Model 3 represents a high solar activity case. The respective mean spot filling factors are 0.3, 1.9 and 9.0 per cent.</figcaption></figure>
and are based on the models of Jeffers (2005), but also include umbral and penumbral regions. We retain the same model numbering for consistency: Model 2 represents solar maximum activity, while Models 3 and 4 are for extrapolated activity levels with higher degrees of spot coverage.

## 2.2 Line profile modelling and RV recovery

We used our spot models to investigate the effect of activity on precision RV measurements of M dwarfs. We calculated synthetic line profiles using spotted 3D stellar models and calculated the resulting RVs by cross-correlating the profiles. We tested our ability to recover planet-induced RVs by combining our simulated spot-induced RVs with various planetary RV signals. In addition to investigating the effect of spot-induced jitter on M dwarfs, we have also investigated the effect of stellar activity on rapidly rotating G and K dwarfs (Jeffers et al. 2014). These simulations incorporated the effects of convective blueshift and modelled facular regions in addition to starspot activity.

Here, we restrict simulations to starspots alone in order to demonstrate the feasibility of removing spot-induced distortions in slowly rotating stars. Doppler Tomography of Stars (DOTS) is a Doppler imaging code that utilizes a two-temperature model to recover surface brightness distributions of active stars (Collier Cameron 2001). DOTS uses maximum entropy regularization to obtain image solutions that minimize artefacts in the presence of noise; it has been used extensively to recover spot distributions on single stars and binaries. A brief overview has also been given recently in Barnes et al. (2015), with application in particular to late M dwarfs. For an image with *i* pixels, spot filling factors *f*<sub>*i*</sub> are obtained with DOTS, assuming a two-temperature model representing the spots and photosphere. For stars that are rotating slowly and for which insufficient Doppler resolution elements can be obtained across the stellar line profile, the derived images are not particularly informative. The regularized fitted line profiles nevertheless offer the potential for distinguishing between distortion due to spots and line centroid shifts due to an orbiting extrasolar planet.

We cross-correlate the fitted line profiles recovered with DOTS, but crucially, do not incorporate the small RVs induced by the planet in the recovery procedure. Similarly, we cross-correlate the input synthetic profiles, which contain both the spot-induced RV variations and the planet-induced RV variations. Any correlation between the simulated and recovered RV can be attributed to stellar activity modulated at the stellar rotation period and can be subtracted from the observed RVs. This approach will be most effective when the $P_{\text{rot}}$ and $P_{\text{orb}}$ are distinct, but further analysis of line profile moments (Anglada-Escudé et al. 2016; Berdiñas et al. 2016) may be useful for distinguishing more difficult cases. Similarly, the method is applicable when starspots are the main contributor to line profile shape variability above the noise level, since the model is designed to fit absorption line distortions (due to spots) and not line shifts (due to an orbiting planet).

## 3 SIMULATIONS

We investigate whether for typical M stars with significant rotation (i.e. $v \sin i \le 10 \text{ km s}^{-1}$) and un-evolving starspot patterns, we are able to distinguish between both starspot signals and stellar reflex motion due to an orbiting planet with modest observational effort. The issue of starspot evolution is discussed further in Section 5. We first begin by simulating a specific case in order to demonstrate the feasibility of the technique.

<table><caption>**Table 1.** System parameters for simulation based on the Proxima Centauri b system, but with a higher $v \sin i$ and thus $P_{\text{rot}}$.</caption><thead><tr><th></th><th>Star</th><th></th></tr></thead><tbody><tr><td>$v \sin i$</td><td>[km s<sup>-1</sup>]</td><td>5</td></tr><tr><td>Stellar radius [$R_{*}$]</td><td>[R<sub>☉</sub>]</td><td>0.14</td></tr><tr><td>Rotation period [$P_{\text{rot}}$]</td><td>[d]</td><td>1.23</td></tr><tr><td>Axial inclination [$i$]</td><td>[°]</td><td>60</td></tr><tr><td>$T_{\text{phot}}$</td><td>[K]</td><td>3000</td></tr><tr><td>$T_{\text{spot}}$</td><td>[K]</td><td>2700</td></tr><tr><td></td><td>Planet</td><td></td></tr><tr><td>Orbital period [$P_{\text{orb}}$]</td><td>[d]</td><td>11.2</td></tr><tr><td>Stellar reflex ampl. [$K_{*}$]</td><td>[ms<sup>-1</sup>]</td><td>2</td></tr></tbody></table>

### 3.1 Proxima centauri analogue model

We illustrate the procedure by adopting a model with the stellar parameters of Proxima Centauri and the orbital period of Proxima Centauri b reported in Anglada-Escudé et al. (2016). While Proxima Centauri possesses a very low $v \sin i < 0.1 \text{ km s}^{-1}$, we investigate projected equatorial rotation velocities in the $v \sin i = 1–10 \text{ km s}^{-1}$ range following the finding that fully convective M dwarfs are on average moderate rotators (Jenkins et al. 2009). We have chosen to simulate a planetary RV signal with a $P_{\text{orb}} = 11.2$ d orbit that induces a stellar semi-amplitude of $K_{*} = 2 \text{ ms}^{-1}$ (slightly higher than the $1.4 \text{ ms}^{-1}$ reported for Proxima Centauri b). The model parameters are summarized in Table 1. In this section, we use an axial inclination of $i = 60^{\circ}$ and projected equatorial rotation velocity of $v \sin i = 5 \text{ km s}^{-1}$. A rotation period of $P_{\text{rot}} = 1.23$ d was calculated from the estimated $R_{*} = 0.14 R_{\odot}$ of Proxima Centauri using $P_{\text{rot}} = 2\pi R_{*} \sin i / v \sin i$.

### 3.2 Simulated line profiles and planet recovery

#### 3.2.1 S/N and number of observations

The maximum entropy procedure requires high S/N ratios to enable recovery of spot features from the line distortions. This can be achieved either by considering many lines simultaneously or by applying least squares deconvolution (Donati et al. 1997) to typically several thousand lines in each spectrum to obtain a single high S/N ratio line profile. Deconvolution has a computational benefit compared with simultaneous fitting of multiple lines. Unlike a straightforward cross-correlation profile, the least squares deconvolved profile accounts for line blending effects. We applied our implementation of least squares deconvolution (Barnes et al. 1998, 2012) to the 0.64–1.03 µm region of M dwarf spectra (observed with UVES) with mean S/N ratios of $25 \lesssim S/N_{\text{obs}} \lesssim 140$. This yielded deconvolved profiles with $2000 \lesssim S/N_{\text{decon}} \lesssim 12000$ and a mean $S/N_{\text{decon}} \sim 5300$ (see table 1 of Barnes et al. 2014). Although $S/N_{\text{obs}}$ in the above range can be achieved for the brightest mid-M dwarfs in the reddest orders with HARPS, red-optical and near-infrared surveys are required to enable sufficient wavelength coverage to achieve $S/N_{\text{decon}} \ge 2000$ on a large sample of mid-late M dwarf targets.

The simulations presented here consider typical line profile $S/N_{\text{decon}}$ ratios of 2000, 5000 and 10000 that would be expected from spectra with respective mean $S/N_{\text{obs}} \sim 25, 60$ and 120. We used the least squares deconvolved profile of the slow rotator, GJ 1061, observed with spectral resolution, $R \sim 100000$ (Barnes et al. 2014), to represent the unbroadened local intensity profile during modelling. Line profiles with starspot asymmetries and an additional
sinusoidal velocity variation specified by the planet observables, $P_{\text{orb}}$ and $K_*$ (i.e. a circular orbit is assumed), were generated with DOTS assuming one observation per night for $N_{\text{obs}} = 15, 30$ and 60 nights. A 3 h random time shift on each observation was applied to minimize day aliases.

### 3.2.2 Radial velocities

The simulated line profiles containing a planet RV signal and line profile distortions due to spots were cross-correlated using the mean profile (i.e. sum of all profiles) as a reference. For the Proxima Centauri analogue model, with $S/N_{\text{decon}} = 5000$, the resulting simulated RVs, $RV_{\text{sim}}$, are shown by the red triangles in Fig. 2 (left-hand panels) using Models 2, 3 and 4 with 15 (top), 30 (middle) and 60 nights (bottom) of observations. For each data set, we used DOTS to fit the simulated line profiles, but without allowing for the sinusoidal variation of the line centroids due to an orbiting planet (i.e. in contrast to Petit et al. 2015 who simultaneously recovered the planet RVs in the presence of spots). The simulated line profiles were fit using the input parameters, $v \sin i, i$ and $P_{\text{rot}}$. The fitted or recovered line profiles thus only contain the line profile distortions due to starspots. These were then cross-correlated against the mean recovered line profile to obtain the recovered velocities, $RV_{\text{rec}}$, shown by the green dashed lines in Fig. 2 (left-hand panels). For Models 3 and 4, the line asymmetries induced by the spots are the dominant contributor to the RVs. Since the spots in Model 2 are too small to be resolved, the fitting procedure is not able to model the line profiles adequately. As a result, for Model 2, the $RV_{\text{rec}}$ points do not well match the $RV_{\text{sim}}$ points. The corrected RVs are obtained by subtracting the recovered RVs from the simulated RVs, where $RV_{\text{corr}} = RV_{\text{sim}} - RV_{\text{rec}}$ (blue circles). The middle panels of Fig. 2 show the recovered RVs ($RV_{\text{rec}}$) versus the simulated RVs ($RV_{\text{sim}}$) showing a linear trend in all instances. The shallow slope for Model 2 is a reflection of the inability of the modelling process to recover the small spots.

### 3.2.3 Periodograms and false alarm probabilities

The right-hand panels of Fig. 2 show the periodograms of the simulated RVs and the corrected RVs, $RV_{\text{sim}}$ and $RV_{\text{corr}}$, for each simulation. Tick marks indicate the simulated periods, $P_{\text{rot}} = 1.23$ d and $P_{\text{orb}} = 11.2$ d, and 0.1, 1 and 10 per cent false alarm probability (FAP) levels are shown. A number of additional peaks arising from the window function sampling and aliasing are present (e.g. $P_{\text{rot}}/2 = 0.615$ d peak and the 1 d sampling beat period with $P_{\text{rot}}$ at 5.3 d). However, it is noteworthy that the solar maximum activity levels from Model 2 are not sufficient to mask the $P_{\text{orb}} = 11.2$ d period, even with only 15 observations. This is expected since in fig. 3 of Barnes et al. (2011), we showed that for the low-contrast case, $<1 \text{ ms}^{-1}$ jitter is expected from Model 2 for $v \sin i = 5 \text{ km s}^{-1}$ in the absence of limiting photon noise. For Models 3 and 4, the $P_{\text{orb}}$ peak only appears at low significance in the periodograms of the simulated RVs (red dashed curves). For the Model 3 and 4 spot corrected $RV_{\text{corr}}$ data, the $P_{\text{rot}}$ peak is no longer present in the periodograms (blue solid curves). At the same time, the $P_{\text{orb}}$ peak significance is boosted. Other peaks that coincide with beating between the 1 d sampling and $P_{\text{orb}}$ however appear at close to ~1 d, with ~2 d aliases. In the case of the $N_{\text{obs}} = 60$ epoch simulation, these peaks are of comparable or greater significance than $P_{\text{orb}}$. Thus, while we can effectively remove the spot jitter using only relatively few observations, the $P_{\text{orb}} = 11.2$ d peak is still of

low significance, reaching only ≳10 per cent FAP for Models 2 and 3 with 30 observations. Fig. 3 (left-hand panel) shows that for our simulated models with $v \sin i = 5 \text{ km s}^{-1}$, $P_{\text{orb}}$ is recovered with <0.1 per cent FAP for Models 2 and 3 with $N_{\text{obs}} \sim 54$, while $N_{\text{obs}} = 58$ is sufficient to recover $P_{\text{orb}}$ with <1 per cent FAP for Model 4.

## 3.3 Identification of the stellar rotation period

The recovered spot RVs in Section 3.2.2 were obtained by fitting the simulated line profiles with the simulated (i.e. known) rotation period, $P_{\text{rot}}$. However, identification of the rotation period is not necessarily straightforward; a $P_{\text{rot}} = 1.23$ d peak is not always significant, or the most significant peak in the periodograms in Fig. 2 (red/dashed curves). The $P_{\text{rot}} = 1.23$ d peak strength depends on $N_{\text{obs}}$ and the number of spots. The window function, aliasing and possible beat periods result in peaks at other periods (see Section 3.2). Rather than relying on Lomb–Scargle periodogram analysis, we used DOTS to search for the simulated stellar rotation period. Since maximum entropy regularized line fitting is an iterative process, the $\chi^2$ achieved after a fixed number of iterations can be used as a criterion for recovering the best-fitting parameters. This procedure works well for the stellar rotation period since spots are periodically visible at the same location on the star as viewed by an observer. Fig. 4 shows that even with a 15 epoch simulation, $P_{\text{rot}} = 1.23$ d is recovered with the minimum $\chi^2$ (i.e. proportional to the maximum likelihood). With $N_{\text{obs}} = 60$, confidence in the period is increased relative to other local $\chi^2$ minima. $P_{\text{rot}}$ can thus be determined with a modest number of observations. Because the periodicities are determined directly from comparing information in the line profiles with a physical model, it is possible to recover periodicities more confidently than with a straightforward Lomb–Scargle periodogram search for sinusoidal signals. Nevertheless, using an armoury of diagnostic tools, including periodogram analyses (i.e. Fig. 2 generally shows a peak at $P_{\text{rot}}$) and photometry, will offer the most effective means of recovering $P_{\text{rot}}$.

## 3.4 Rotation and stellar inclination

In addition to recovering $P_{\text{rot}}$, the projected equatorial rotation velocity, $v \sin i$, must be recovered. It is common practice in Doppler imaging to optimize $v \sin i$, which correlates with line strength when using a fixed line profile to represent the local stellar spectrum (e.g. see Barnes et al. 2000 and Collier Cameron 2001). This procedure is often needed for each data set to optimize the fit and minimize image artefacts. These parameters are generally easily determined by $\chi^2$ minimization and are orthogonal to periodic signatures. In addition, the stellar axial inclination, $i$, must be specified during the starspot modelling and line profile recovery. The stellar radius cannot generally be measured directly from observation, which means that $i$ is not generally known. We performed a $\chi^2$ search for a range of stellar axial inclinations, $i_{\text{rec}}$, to determine the effect of an incorrectly determined axial inclination. A minimum $\chi^2$ is always attained at the simulated combination of $v \sin i$ and $P_{\text{rot}}$, irrespective of adopted axial inclination for $30^{\circ} \le i_{\text{rec}} \le 90^{\circ}$. These parameters are therefore effectively independent of the choice of $i_{\text{rec}}$. However, for the trial $i_{\text{rec}}$ values, only 2 per cent r.m.s. variation in $\chi^2_{\text{min}}$ was seen, with no discernible trend in $\chi^2$ versus $i_{\text{rec}}$ in the $30^{\circ}-90^{\circ}$ range above the scatter. Thus, although the simulated inclination of $i = 60^{\circ}$ cannot be recovered for $v \sin i = 5 \text{ km s}^{-1}$, this does not preclude identification of $v \sin i$ and $P_{\text{rot}}$, which are correctly recovered at all trial $i_{\text{rec}}$ values.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Simulated RV curves comprising a $P_{\text{orb}} = 11.2$ d planet that induces a stellar semi-amplitude of $K_{*} = 2$ ms⁻¹ in the presence of stellar modulations from starspot Models 2–4 (see Fig. 1). Stellar parameters are $v \sin i = 5$ km s⁻¹, $P_{\text{rot}} = 1.23$ d and $i = 60^{\circ}$. Simulations for $N_{\text{obs}} = 15$ (top), 30 (middle) and 60 nights (bottom) are plotted. Left: the simulated spot + planet RVs and maximum entropy recovered RVs are respectively denoted $RV_{\text{sim}}$ and $RV_{\text{rec}}$. The corrected RVs are $RV_{\text{corr}} = RV_{\text{sim}} - RV_{\text{rec}}$. The model planetary RVs are shown by the (grey) sine curve. Centre: plots showing the correlation between $RV_{\text{rec}}$ and $RV_{\text{sim}}$. Right: Lomb–Scargle periodograms for $RV_{\text{sim}}$ (simulated) and $RV_{\text{corr}}$ (corrected) RVs. Dashed vertical lines indicate $P_{\text{rot}}$ and $P_{\text{orb}}$. False alarm probability levels, FAP = 0.001, 0.01 and 1 (i.e. 0.1, 1 and 10 per cent) are indicated by the horizontal lines (top to bottom, respectively).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** FAP of $P_{\text{orb}}$ for Models 2–4 after removing the modelled starspot contribution. Left: FAP versus number of observations, $N_{\text{obs}}$. Right: for the $N_{\text{obs}} = 60$ case, $P_{\text{orb}}$ FAP versus the stellar axial inclination used to recover the starspot contribution, $i_{\text{rec}}$. The input simulation of $i = 60^{\circ}$ is indicated with the vertical dashed line.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Stellar rotation period recovery with DoTS for Model 3 with 15 (red/lower curve) and 60 (blue/upper curve) observations. $\chi^2$ normalized to the minimum recovered value is plotted versus stellar rotation period. $P_{\text{rot}} = 1.23$ d is shown by the solid vertical line. Solid tick marks denote $P_{\text{rot}/2}$, $2P_{\text{rot}}$, $3P_{\text{rot}}$ and $4P_{\text{rot}}$. The dashed vertical line and tick marks indicate $P_{\text{orb}} = 11.2$d, $P_{\text{orb}/2}$, $P_{\text{orb}/3}$ and $P_{\text{orb}/4}$.</figcaption></figure>

A robustly determined rotation period, $P_{\text{rot}}$ and $v \sin i$ still enable $i$ to be estimated for an assumed stellar radius. The typical 10 per cent uncertainty for M dwarf estimates of $R_*$ (Stassun et al. 2012; Williams et al. 2015) will dominate uncertainty in the axial inclination determination. For instance, uncertainties of 10 per cent in $\sin i$ yield $i = 30^{+3.4}_{-3.3}^{\circ}$, $45^{+6.1}_{-5.5}^{\circ}$, $60^{+12.3}_{-8.8}^{\circ}$, $75^{+15}_{-9.6}^{\circ}$ and $90^{+0}_{-25.8}^{\circ}$ (upper limits at $75^{\circ}$ and $90^{\circ}$ arise from requiring $\sin i \le 1$). For RVsim data using $i = 60^{\circ}$ and $N_{\text{obs}} = 60$, Fig. 3 (right) shows the FAP for RVcorr data points using $i = 30^{\circ}-90^{\circ}$ (in $15^{\circ}$ steps). In fact, there is a tendency towards improved FAPs when using higher inclinations (i.e. $i_{\text{rec}} = 75^{\circ}$ and $90^{\circ}$) for recovery. This is likely related to the lack of resolution elements, when instrumental resolution and stellar $v \sin i$ are of comparable magnitude.

## 4 SENSITIVITIES

As every object and data set combination is unique, a full exploration of parameter space is beyond the scope of this paper. Nevertheless, a more complete picture of the range of sensitivities can be gained by exploring the S/N and $v \sin i$ values that are expected empirically. Fig. 5 shows the reduction in starspot noise achieved, $\Gamma = \rho_{\text{sim}}/\rho_{\text{corr}}$ (where $\rho_{\text{sim}}$ and $\rho_{\text{corr}}$ are the r.m.s. for the RVsim and RVcorr data points), as a function of $v \sin i$ for the simulated S/N ratios. The lowest values of $\Gamma$ are found for Model 2, where $\Gamma \sim 1$. Again,’ this is due to an inability to recover the small spots, and in some combinations,

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Starspot noise reduction factor, $\Gamma$, as a function of $v \sin i$ for S/N = 2000, 5000 and 10 000 (bottom to top).</figcaption></figure>

the noise may be marginally increased by ~5 per cent. The procedure is most efficient for $v \sin i = 5$ km s$^{-1}$, where a maximum $\Gamma = 8.6$ is achieved. At $v \sin i = 1$ km s$^{-1}$, the efficiency is generally lower, as fitting distortions below the instrumental resolution is less effective. For $v \sin i = 10$ km s$^{-1}$, although starspot distortions are modelled well, the broader lines due to higher $v \sin i$ limit the efficiency of the procedure by increasing the RV uncertainties.

For projected equatorial rotation velocities, $v \sin i = 1, 2, 5$ and $10$ km s$^{-1}$, we expect $P_{\text{rot}} = 6.13, 3.10, 1.23$ and $0.613$ d for a $0.14$ R⊙ star with $i = 60^{\circ}$. With $P_{\text{orb}}$ fixed at 11.2d, we have estimated the minimum $K_*$ recovered with $\le 0.1$ FAP using $N_{\text{obs}} = 60$ before and after subtracting the fitted starspot contributions. A noise floor of $1$ ms$^{-1}$ is assumed and we simulated recovery FAPs for a number of discrete values of $K_*$ up to $50$ ms$^{-1}$. The increase in sensitivity with increasing S/N ratio and decreasing $v \sin i$ that one intuitively expects is seen in Fig. 6. The connected points indicate the sensitivity for each model after removing the spot contribution. The vertical arrows indicate the change in sensitivity from the RVsim points to the RVcorr points. The limiting precision is also shown by the grey dot-dashed line, and is derived from line profiles generated with an unspotted model. The absorption lines are more effectively cleaned for Model 4 since these spots induce the largest amplitudes and are most easily fit. The procedure is ineffective for Model 2 because the spots are not resolved. For Model 2, when the S/N ratio is low and $v \sin i$ is $10$ km s$^{-1}$, the modelling adds noise and decreases the sensitivity slightly.

For the fixed $P_{\text{orb}} = 11.2$d, in the case of the slowest rotators with $v \sin i = 1-2$ km s$^{-1}$, $N_{\text{obs}} = 60$ should enable detection of $1-2$ M$_{\oplus}$ planets, and is limited to $2$ M$_{\oplus}$ only for S/N = 2000.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Minimum detected amplitude and planet mass with FAP $\le 0.001$ (0.1 per cent) for $S/N = 2000, 5000$ and $10000$ (bottom to top). For Models 2, 3 and 4, the minimum detected $K_*$ and equivalent planet mass, $m_{\text{plan}}$, is plotted for the simulated $v \sin i$ values after applying spot correction (with an assumed $i = 60^\circ$, $K_*/m_{\text{plan}} \sim 1$ for the adopted $M_*$ and $P_{\text{orb}}$). The vertical arrows indicate the change in sensitivity achieved by subtracting the fitted starspot RVs. Before correction, the sensitivity limits for $v \sin i = 10\text{km s}^{-1}$ are $K_* = 50 \text{ms}^{-1}$ (i.e. $M_{\text{plan}} \sim 50 M_\oplus$). The grey dot-dashed line is the precision limit for an unspotted star.</figcaption></figure>

For $v \sin i = 5\text{km s}^{-1}$ and the highest $S/N_{\text{decon}} = 10000$, $1 M_\oplus$ planets can potentially still be detected if Model 2 spot levels are realistic. More typical limits of $2-4 M_\oplus$ are achieved for Models 3 and 4 after removal of spot RVs. Once $v \sin i$ reaches $10 \text{km s}^{-1}$, only $\ge 4 M_\oplus$ planets are likely to be detected for the simulated observations, assuming Models 3 and 4 are representative of spot coverage.

## 5 SUMMARY AND DISCUSSION

Our ability to recover and remove stellar activity induced RVs is demonstrated for late-M dwarfs with distributed starspot patterns. We have not attempted a complete exploration of parameter space. None the less, we are confident that a typical habitable zone planet, bounded by the runaway and maximum greenhouse conditions, with $a = 0.042-0.082 \text{au}$ and $P_{\text{orb}} \sim 9.1-24.5 \text{d}$ (Delfosse et al. 2000; Kopparapu et al. 2013a,b), could be detected orbiting an M6V analogue of Proxima Centauri (Anglada-Escudé et al. 2016). The size of the line distortions induced by the spots, the resolving power and the S/N of the observed spectra are the most important factors that determine the efficacy of the method. Our adopted spot models are extrapolated solar spot size distributions, with spot contrasts and surface distributions based on observations of rapidly rotating M dwarfs, although we do not know the exact form of the spot distributions for moderately rotating, fully convective M

dwarfs. Nevertheless, with reasonable levels of spot coverage, we are confident that spot activity can be reduced to a level that enables $1 M_\oplus$ signals to be detected in spectra observed with typical $S/N_{\text{obs}} \sim 60$ (corresponding to the simulated deconvolved profiles with $S/N_{\text{decon}} \sim 5000$ in this paper).

Observations that aim to detect planets and determine planet occurrence rates for fully convective M dwarfs will only be unbiased if moderate rotators are surveyed. This means removing spots is crucial. We find that for moderately spotted stars, $2 M_\oplus$ planets can be detected for $v \sin i = 5 \text{km s}^{-1}$, while $4 M_\oplus$ planets can be detected at $v \sin i = 10 \text{km s}^{-1}$. Without corrections, the mass limits are up to a factor of $\Gamma = 8.6$ or higher in the I band. Applying the technique to M dwarfs with $v \sin i > 5 \text{km s}^{-1}$ at $S/N_{\text{obs}} \le 25$ (yielding $S/N_{\text{decon}} = 2000$) will not necessarily be advantageous if spot coverage is modest (i.e. Models 2 and 3 with filling factors of <2 per cent). Clustering of spots in groups potentially results in larger amplitude starspot distortions (i.e. where spots within a spot group are not resolved in the line profile) that are more easily fit. Although Doppler images of rapidly rotating M dwarfs (i.e. with $v \sin i \ge 20 \text{km s}^{-1}$) reveal distributed starspot patterns, it is less clear whether spot activity in stars with $v \sin i \le 10 \text{km s}^{-1}$ will show similar distributed spot patterns, or spots located in active regions, as seen on the Sun and other earlier G and K dwarfs.

Obtaining rotation periods for the latest M stars is troublesome owing to the low contrast and distributed nature of the spots. As demonstrated in Section 3.3, identification of the planetary period from other short-period alias peaks in Fig. 2 is likely to require careful consideration of a combination of photometry, activity parameters and line diagnosis tools. We have shown that using maximum entropy fitting offers one such method for reliable determination of the rotation period, where periodogram analysis alone does not yield an unambiguous peak. Further, for real observations with additional systematics, traditional Lomb–Scargle periodogram analyses do not offer the optimal path to recovery of planet-induced periodicities. Frequentist and Bayesian approaches that incorporate the stellar signal into the period analysis (Baluev 2013; Tuomi et al. 2014; Anglada-Escudé et al. 2016), using additional instrumental and atmospheric noise priors, will provide a better means of assessing the periodicities.

For slowly rotating stars, where $v \sin i \ll R$ (the instrumental resolution), the maximum entropy method presented here will not be effective. For $v \sin i < 1 \text{km s}^{-1}$ even if spot distributions similar to Model 3 or 4 were seen, it is unlikely that they will be resolved since micro- and macro-turbulence will dominate the local intensity profile width. Additional tools for assessing line shape behaviour can also be applied to assess correlations with instrumental systematics and activity (Anglada-Escudé et al. 2016; Berdiñas et al. 2016). For instance, although late-M field stars typically exhibit significant rotation, Proxima Centauri, with a probable $v \sin i < 0.09 \text{km s}^{-1}$, has demonstrated line profile (Anglada-Escudé et al. 2016) variability on an ~80 d time-scale, matching the photometric modulation. The origin of the variability in M dwarf stars is also not clear. The variability in the second central moment (essentially the width of the spectral lines), was shown by Anglada-Escudé et al. (2016) to be anti-correlated with the photometry. In other words, the line width is greatest when the stellar flux is lowest. This implies potential effects from plage regions, which are likely associated with spot activity, and which will also modify the line equivalent width. Further modelling to enable the recovery of line variations due to plage may therefore be appropriate as an extension to the assumed two-temperature (photosphere + spot) simulations that we present in this work.
Recovery of surface temperature inhomogeneities will be further complicated by the fact that starspots are not static but constantly evolving with time. Spot groups on dwarf stars are nevertheless stable on time-scales of months (Barnes et al. 1998; Bradshaw & Hartigan 2014), while coherence of individual spots on day-week time-scales has been used to measure differential rotation (e.g. Barnes et al. 2005; Collier Cameron & Donati 2002). This is also true for fully convective active M dwarfs that show spot stability over several days (Barnes et al. 2016). Differential rotation with small amplitude or consistent with solid body rotation has been found by Morin et al. (2008a,b) and Barnes et al. (2016) in these targets, and may be responsible for the stability of spots on even longer time-scales. RV surveys that target stars intensively, with observations obtained in relatively short campaigns are thus likely to be the most successful. The campaign that identified Proxima Centauri b (Anglada-Escudé et al. 2016) obtained ~60 observations during a 2-month timespan. On short time-scales, peak aliasing is minimized, while stellar activity variability can more easily be monitored and accounted for. This contrasts with typical multi-epoch surveys that often obtain a few observations per season over several years. We have demonstrated that the Doppler imaging method works well, even with only 15 observations that sample the range of rotation phases, at relatively high cadence over a short observing baseline. Splitting a 60 epoch campaign into subgroups of 15 observations would enable the effects of spots to be mitigated, while further minimizing the effects of spot evolution.

The simulations presented indicate that for stars with distributed starspots at low contrast, removal of starspot-induced RVs is challenging. Nevertheless, the technique presented in this paper demonstrates that we can model and correct for starspot-induced RVs in stars with $v \sin i = 1–10$ km s⁻¹. With higher spot contrasts, the technique can be applied to less active and more slowly rotating stars. Since no prior assumptions about the spot distributions are required, RV surveys aimed at detecting earth-mass planets at all spectral types will benefit from this technique. We have clearly demonstrated that high cadence observations on the time-scale of the stellar rotation period are essential for reliable RV detection of planets orbiting active stars.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for taking the time to review this manuscript. JRB and CAH were supported by the Science and Technology Facilities Council (STFC) under the grant ST/L000776/1. SVJ acknowledges research funding by the Deutsche Forschungsgemeinschaft (DFG) under grant SFB 963/1, project A16. HRAJ MT and FF are supported by the Leverhulme Trust grant, RPG-2014-281, and the STFC grant, ST/M001008/1. JSJ acknowledges funding by Fondo Nacional de Desarrollo Científico y Tecnológico (FONDECYT) through grants 1161218 and 3110004, and partial support from CATA-Basal (PB06, Comisión Nacional de Investigación Científica y Tecnológica (CONICYT)), the GEMINI-CONICYT FUND and from the Comité Mixto ESO-GOBIERNO DE CHILE.

## REFERENCES

* Anglada-Escudé G. et al., 2016, Nature, 536, 437

* Baluev R. V., 2013, MNRAS, 429, 2052

* Barnes J. R., Collier Cameron A., 2001, MNRAS, 326, 950

* Barnes J. R., Collier Cameron A., Unruh Y. C., Donati J. F., Hussain G. A. J., 1998, MNRAS, 299, 904

* Barnes J. R., Collier Cameron A., James D. J., Donati J.-F., 2000, MNRAS, 314, 162

* Barnes J. R., James D. J., Cameron A. C., 2004, MNRAS, 352, 589

* Barnes J. R., Cameron A. C., Donati J.-F., James D. J., Marsden S. C., Petit P., 2005, MNRAS, 357, L1

* Barnes J. R., Jeffers S. V., Jones H. R. A., 2011, MNRAS, 412, 1599

* Barnes J. R. et al., 2012, MNRAS, 424, 591

* Barnes J. R. et al., 2014, MNRAS, 439, 3094

* Barnes J. R., Jeffers S. V., Jones H. R. A., Pavlenko Y. V., Jenkins J. S., Haswell C. A., Lohr M. E., 2015, ApJ, 812, 42

* Barnes J. R., Jeffers S. V., Jones H. R. A., Pavlenko Y. V., Jenkins J., Haswell C., Lohr M., 2016, in Feiden G., ed., The 19th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun (CS19). Uppsala, Sweden

* Berdiñas Z., Amado P., Anglada-Escudé G., Rodríguez-López C., Barnes J., 2016, MNRAS, 459, 3551

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

* Collier Cameron A., 2001, in Boffin H. M. J., Steeghs D., Cuypers J., eds., Lecture Notes in Physics, Vol. 573, Astrotomography – Indirect Imaging Methods in Observational Astronomy. Springer-Verlag, Berlin, p. 183

* Collier Cameron A., Donati J.-F., 2002, MNRAS, 329, L23

* Cosentino R. et al., 2012, Proc. SPIE, 8846, 1

* Delfosse X., Forveille T., Ségransan D., Beuzit J.-L., Udry S., Perrier C., Mayor M., 2000, A&A, 364, 217

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

* Dumusque X., Boisse I., Santos N. C., 2014, ApJ, 796, 132

* Gao P. et al., 2016, PASP, 128, 104501

* Jeffers S. V., 2005, MNRAS, 359, 729

* Jeffers S. V., Barnes J. R., Jones H. R. A., Reiners A., Pinfield D. J., Marsden S. C., 2014, MNRAS, 438, 2717

* Jenkins J. S., Ramsey L. W., Jones H. R. A., Pavlenko Y., Gallardo J., Barnes J. R., Pinfield D. J., 2009, ApJ, 704, 975

* Kopparapu R. K. et al., 2013a, ApJ, 765, 131

* Kopparapu R. K. et al., 2013b, ApJ, 770, 82

* Mahadevan S. et al., 2014, Proc. SPIE, 9147, 1

* Mayor M. et al., 2003, The Messenger, 114, 20

* Mohanty S., Basri G., 2003, ApJ, 583, 451

* Morin J. et al., 2008a, MNRAS, 384, 77

* Morin J. et al., 2008b, MNRAS, 390, 567

* Moulds V. E., Watson C. A., Bonfils X., Littlefair S. P., Simpson E. K., 2013, MNRAS, 430, 1709

* Newton E. R., Irwin J., Charbonneau D., Berta-Thompson Z. K., Dittmann J. A., West A. A., 2016, ApJ, 821, 93

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

* Phan-Bao N., Lim J., Donati J.-F., Johns-Krull C. M., Martín E. L., 2009, ApJ, 704, 1721

* Quirrenbach A. et al., 2010, Proc. SPIE, 7735, 13

* Reiners A., Basri G., 2010, ApJ, 710, 924

* Reiners A., Bean J. L., Huber K. F., Dreizler S., Seifahrt A., Czesla S., 2010, ApJ, 710, 432

* Reiners A., Joshi N., Goldman B., 2012, AJ, 143, 93

* Rockenfeller B., Bailer-Jones C. A. L., Mundt R., 2006, A&A, 448, 1111

* Stassun K. G., Kratter K. M., Scholz A., Dupuy T. J., 2012, ApJ, 756, 47

* Suárez Mascareño A., Rebolo R., González Hernández J. I., Esposito M., 2015, MNRAS, 452, 2745

* Thibault S. et al., 2012, Proc. SPIE, 8446, 30

* Tuomi M., Jones H. R. A., Barnes J. R., Anglada-Escudé G., Jenkins J. S., 2014, MNRAS, 441, 1545

* Williams P. K. G., Berger E., Irwin J., Berta-Thompson Z. K., Charbonneau D., 2015, ApJ, 799, 192

This paper has been typeset from a TeX/LATeX file prepared by the author.