The effect of M dwarf starspot activity on low-mass planet detection thresholds

Abstract

In light of the growing interest in searching for low mass, rocky planets, we investigate the impact of starspots on radial velocity searches for earth-mass planets in orbit about M dwarf stars. Since new surveys targeting M dwarfs will likely be carried out at infrared wavelengths, a comparison between V and Y band starspot-induced jitter is made, indicating a reduction of up to an order of magnitude when observing in the Y band. The exact reduction in jitter is dependent on the photosphere to spot contrast ratio, with greater improvements at smaller contrasts.

We extrapolate a model used to describe solar spot distributions to simulate the spot patterns that we expect to find on M dwarfs. Under the assumption that M dwarfs are near or fully convective, we randomly place starspots on the stellar surface, simulating different levels of spot coverage. Line profiles distorted by spots are derived and are used to investigate the starspot-induced jitter. By making assumptions about the degree of spot activity, detection limits for earth-mass planets in habitable zones are simulated for between 10 and 500 observation epochs. We find that ≤50 epochs are required to detect 1–2 planets (with <1 per cent false alarm probability) orbiting slowly rotating 0.1 and 0.2 M_⊙ stars. This sensitivity decreases when typical rotation velocities and activity levels for each stellar mass/spectral type are considered. No detections of below 20 planets are expected for ≤500 observations for the most active stars with v sin i≥ 20 km s⁻¹ and dark spots.

techniques: radial velocities, stars: activity, stars: atmospheres, planetary systems, starspots

1 INTRODUCTION

Nearly 500 extrasolar planets in over 400 planetary systems have been discovered,¹ and yet <6 per cent of those with known or estimated host star masses are M dwarf planetary systems. Since M dwarf masses may be as low as one-tenth of a solar mass, it is not unreasonable to expect that their most massive planets might also be correspondingly smaller (Ida & Lin 2005) and thus more difficult to detect than the ubiquitous gas giants found orbiting F, G and K dwarf stars. On the other hand, planets of a fixed mass orbiting lower mass stars become easier to detect, owing to the reduced mass ratio of the system. Detecting low-mass planets (which are predicted in significant numbers by Ida & Lin 2005) is however no easy task, since the planet-induced stellar velocity amplitudes are small, even for M dwarfs. Additionally, in the case of M dwarfs, there is an observational bias against detecting planets at optical wavelengths (where most searches have so far been carried out) owing to the relatively low fluxes compared with redder wavelengths. Observations carried out at near-infrared wavelengths, where the host M stars are several magnitudes brighter than at shorter wavelengths, open up the possibility of carrying out larger surveys of lower mass stars that are capable of detecting low-mass planets.

Radial velocity (RV) surveys are still by far the most successful means of detecting planets. Despite the difficulties of observing in the infrared, there are a number of fledgling projects aimed at carrying out planet searches. Both Ramsey et al. (2008) and Steinmetz et al. (2008) showed that high precision (<10 m s⁻¹) RV measurements could be made at short infrared wavelengths via observations of the Sun. Seifahrt & Käufl (2008) demonstrated that similar precision was possible on short time-scales using the cryogenic high-resolution infrared échelle spectrograph (CRIRES), operating at L-band wavelengths, at the Very Large Telescope (VLT). More recently, an ammonia gas cell has been employed to achieve ∼5 m s⁻¹ precision (Bean et al. 2010b) over 6 month time-scales with CRIRES operating at K-band wavelengths. This survey (Bean et al. 2010a) has not found evidence for the massive planet claimed to be orbiting the low-mass M dwarf VB 10 (Pravdo & Shaklan 2009; Zapatero Osorio et al. 2009).

With infrared multi-order spectrographs working at high resolution (R > 50 000), precision RVs of order 1 m s⁻¹ are expected to be achieved in the next few years. Instrumental precision of 1 m s⁻¹ is a desirable goal when we consider that a 1 planet orbiting in the habitable zone of 0.1 M_⊙ star induces a ∼1.7 m s⁻¹ RV amplitude. Rivera et al. (2005) have nevertheless already estimated ≲1.5 m s⁻¹ RV jitter for the M4 dwarf, GJ 876 (at optical wavelengths).

A major cause of RV jitter is the presence of magnetic activity induced starspots that distort absorption line profiles (Saar & Donahue 1997). Any time variant asymmetries in line profiles due to starspots will lead to biased measurements of dynamically induced RVs from orbiting bodies such as planets. If the RV signal of a planet is of similar magnitude or smaller than spot-induced jitter, any candidate planet signal may not be recovered with reasonable allocations of observing time. In this paper, we simulate semirealistic starspot patterns for M dwarfs. By making use of starspot size distribution models for the Sun and by considering indirectly observed starspot distribution patterns for M dwarf stars, we model the subsequent line profiles from which RV jitter is determined via cross-correlation. This approach extends similar work carried out by Desort et al. (2007) and Reiners et al. (2010) to more realistic spot distributions that are solar spot distributions extrapolated to active stars. The method is more analogous to the recent work by Lagrange, Desort & Meunier (2010) that investigates the effect that solar sunspot activity would have on detection thresholds of an earth-like planet orbiting in the habitable zone. Here, we use our models to determine the detection thresholds for earth mass habitable zone planets orbiting M dwarf stars that exhibit different equatorial rotation velocities and starspot activity levels. In Section 2 we discuss evidence for M dwarf starspot distributions and introduce our extrapolated solar model. In Section 3, we briefly describe the methods used to generate line profiles. A comparison of the relative spot amplitudes for different photosphere-to-spot temperature contrast ratios (T_p/T_s) is made in Section 3.2. The RVs induced by stars that exhibit the spot distributions introduced in Section 2.2 are then investigated in Section 3.3. Here, we investigate the effects of spots for three stellar masses at extremes of the expected T_p/T_s values. We then carry out detection threshold simulations for earth-mass planets in Section 4 before finally considering possibilities that generalize our specific model cases further in Section 5.

2 M DWARF SPOT MODELLING

2.1 Observed starspot distributions

Compared with earlier spectral types, our knowledge of starspot distribution patterns on M dwarfs is less complete. For earlier spectral types, publications that predict and report observations of these distributions are numerous. The solar analogue is our chief reference point for which we observe spots appearing chiefly at 0°–40° latitudes. There is also considerable evidence that other stars exhibit similar starspot activity to the Sun. Observations made by the Mount Wilson Survey (Baliunas et al. 1995), which has observed changes in the rotation periods of a number of stars over decades (Donahue, Saar & Baliunas 1996), suggest that similar magnetic dynamo processes are at work. The appearance of spots within defined (low) latitude bands has been attributed to the interface dynamo process by which the magnetic fields responsible for the spots are generated. A magnetic dynamo located at the boundary between the radiative core and convection zone together with radial transport of flux through convection is able to explain the appearance of photospheric flux at low to mid latitudes only (Moreno-Insertis, Schüssler & Ferriz-Mas 1992). For more rapidly rotating stars which can be indirectly imaged via the Doppler imaging process, this scenario (Schüssler et al. 1996) is less successful at predicting the starspot distributions (Barnes et al. 1998). However spots still often appear within specific latitude regions in Doppler images.

Unlike G and K stars, owing largely to their inherent faintness, M dwarfs are less well studied. By mid-M, it is predicted that stars become fully convective and thus the standard αΩ interface dynamo process can no longer operate. There is nevertheless a great deal of evidence indicating that fully convective M dwarf stars are magnetically active. A study by Gizis et al. (2000) found that the percentage of objects with Hα appearing in emission peaks at late spectral type while Mohanty & Basri (2003) found that the equatorial rotation velocity at which Hα activity saturates appears at higher values for mid-late M spectral types than for early to mid M spectral types. This was however later ruled out by Reiners & Basri (2010). The topic is too large to review at length here where we are concerned primarily with starspot distributions.

Observations of rapidly rotating early M dwarfs via indirect Doppler imaging techniques (Barnes & Collier Cameron 2001; Barnes, James & Cameron 2004) reveal that they are more uniformly covered with spots than earlier spectral types. Donati et al. (2008) and Morin et al. (2008) have found, via magnetic Doppler imaging using Stokes V measurements, that the magnetic topology of M dwarfs change at approximately spectral type M4V. Although a switch to a fully convective dynamo might be expected to lead to a more uniformly spotted star, the results of Donati et al. (2008) and Morin et al. (2008) appear to counter-intuitively indicate that fields become more dipolar. This phenomenon has also been investigated by Reiners & Basri (2009) who find that in fact more than 85 per cent of the magnetic flux is stored in magnetic fields that are invisible to Stokes V. However, using all Stokes components is necessary for a full description of a star's magnetic topology. Additionally, the field from the darkest magnetic regions (e.g. spots) is not visible in Stokes V owing to the large contrast seen at optical wavelengths. The high average fields of a few kG found on M dwarfs (Johns-Krull & Valenti 1996; Reiners & Basri 2007, 2010) are unlikely to be concentrated in small spots because very high local fields would be necessary. They are perhaps evenly distributed across the surfaces of the stars. The traditional picture of spots as the only concentrations of magnetic flux may therefore be misleading in these stars. It is thus clear that neither magnetic imaging nor surface brightness imaging at optical wavelengths are able to give a complete picture of the starspot/magnetic field topology of stars.

While starspot patterns (e.g. the latitudes at which spots appear) on more rapidly rotating G and K stars may vary as a function of rotation, it is not clear whether such changes take place among fully convective M dwarf stars that are expected to generate magnetic fields via a turbulent dynamo process. Hence those surface brightness images derived for M dwarfs (Barnes & Collier Cameron 2001; Barnes et al. 2004) may also be representative of slower rotators. Moreover, at later spectral types, M7–M9.5 for example, Reiners & Basri (2010) find no correlation between rotation and magnetic flux generation. Further evidence that starspots may be more uniformly distributed comes from observations of reduced starspot-induced light-curve variability. Messina et al. (2003) have shown that already by spectral type K6–M4, the maximum starspot-induced photometric variability is around a factor of 2 lower. Other individual studies of mid-M and late-M (M5 and M9) dwarfs (Rockenfeller, Bailer-Jones & Mundt 2006a; Rockenfeller et al. 2006b) also show peak-to-peak amplitudes of order 0.05 or less in G, R and I photometric bands. A more uniform distribution of spots (rather than clustering of spots in one or two active regions) would be one explanation for the observed reduction in photometric-induced light variations in M dwarf stars.

A series of investigations into spot coverage factors (Neff, O'Neal & Saar 1995; O'Neal, Saar & Neff 1996; O'Neal, Neff & Saar 1998; O'Neal et al. 2004) have been carried out using TiO as a tracer of cooler temperatures on G and K stars of different spectral class, including dwarfs. These studies indicate typical spot coverage of 20–50 per cent for active stars. This may seem surprising, but indicates that light-curve analyses and Doppler imaging studies (that typically find of order 10 per cent spot coverage) are not sensitive to an underlying spot distribution. It is again unclear whether these spots arise from a boundary dynamo or whether a turbulent dynamo is responsible. Based on these findings, investigations of high levels of spot coverage seem warranted. A threshold level of spot coverage may be expected (for uniform spot coverage) at which any RV jitter effects no longer increase in magnitude.

With a largely unknown starspot pattern in moderately rotating M dwarfs and the general prediction that a distributed dynamo should produce randomly distributed spots, we carry out a number of simulations to assess the detectability of planets around M stars. Without strong evidence to the contrary, we assume that more rapidly rotating stars are more spotted. Synthetic starspot models are used to generate line profiles and to investigate the RV amplitudes resulting from non-uniform line profiles. These RVs are then used to determine our ability to detect low-mass planets which are in habitable zone orbits around M dwarf stars.

2.2 Modelling randomly distributed spots

We synthesize spot maps which follow a lognormal size distribution on the surface of an immaculate star. The Doppler imaging code ‘Doppler Tomography of Stars’ (dots; Collier Cameron 1997) was then used to place the spots on the surface of a model star. The input parameters to dots for modelling spots were first introduced in Jeffers (2005). They are, where x is a random number in the range (0 ≤x≤ 1),

• longitude: randomly distributed between 0° and 360°;

• latitude: −(π/2) < θ < (π/2), following with 0 ≤x≤ 1, to eliminate an artificial concentration of spots at the pole;

• spot radius: computed using the previously described lognormal distribution as tabulated in Table 1;

• spot brightness and spot sharpness: modelled to obtain an umbral to penumbral ratio of 1:3 (Solanki 1999).

Table 1

Tabulation of the input parameters to the lognormal size distribution of starspots given by equation (1). The parameters derived by Bogdan et al. (1988) for the Sun are data sets 1 and 2, and those calculated by Solanki (1999) for active stars are data sets 3–6 (A _1/2⊙≡ 2πR²_⊙). The spot filling is determined directly from dots and differs from the umbral spot coverage since we also include penumbral regions. The spot coverage and filling are the same for models 1 and 2 since the spot sizes are close to the pixel resolution in these cases. R′_HK values corresponding to the solar min and max cases are given (Lagrange et al. 2010) and the upper limit for starspot coverage derived from TiO studies (e.g. O'Neal et al. 1998) is indicated.

Model	1	2	3	4	5	6
σA(× 10−6 A 1/2⊙)	3.8	5.0	6.8	9.2	12.2	15.8
(d N/d A)max	5	25	65	125	205	305
Umbral spot coverage	0.03 per cent	0.3 per cent	1.6 per cent	6.1 per cent	18 per cent	48 per cent
Total spot filling	0.03 per cent	0.3 per cent	1.9 per cent	9.0 per cent	29.5 per cent	62.4 per cent
Solar:	Min	Max
R′HK:	−5.0	−4.85			–> TiO obsns.

The key difference when compared with the solar case is that spots are allowed at all latitudes. The solar spot size distribution has been determined by Bogdan et al. (1988) from direct observations taken from the Mount Wilson white-light plate collection covering the period 1917–1982. Bogdan et al. (1988) show that the number of sunspots, N, with umbral area, A, is given by

1

where the constants 〈A〉 and σ_A are the mean and geometric standard deviation of the lognormal distribution of the sunspot areas, and is the maximum value reached by the distribution. For the case of the Sun, these values are tabulated in Table 1 where model 1 is for an inactive Sun and model 2 is for an active Sun. Cool, young potentially planet hosting stars are expected to be more active than the Sun. For these stars we use the extrapolations of Solanki (1999) also shown in Table 1. All starspots are modelled, following Solanki (1999) with circular umbral areas. A starspot is defined such that the umbral region is at the temperature of the spot (tabulated in Table 2). Following the assumptions of Solanki (1999), we also include penumbral regions with umbral to penumbral areas of 1:3 (i.e. radii ratios of 1:2). The penumbral regions are defined such that their intensity is equal to half the difference of the photosphere and spot intensities.

Table 2

Simulated stellar mass, photospheric temperature, T_p, spot temperatures (to the nearest 50 K) and K. Also tabulated are the periods of orbiting habitable zone planets for each stellar mass and the stellar RV amplitudes induced by orbiting planets of mass 1, 2, 5, 10 and 20 M⁠.

Stellar mass	Tp			PPlanet	K* for planet masses of
					1 M	2 M	5 M	10 M	20 M
(M⊙)	(K)	(K)	(K)	(d)	(m s−1)
0.1	2750	1800	2550	4.96	1.74	3.48	8.69	17.41	34.76
0.2	3250	2100	3050	13.02	0.79	1.59	3.97	7.94	15.87
0.5	3750	2450	3550	35.6	0.31	0.62	1.54	3.08	6.16

We have included the solar values for R′_HK at solar min and solar max, as derived by Lagrange et al. (2010) who used the relationships of Noyes et al. (1984) and Lockwood et al. (2007). These values also give an indication of the kinds of activity levels that are routinely found in optical planet survey targets (Tinney et al. 2002). We also indicate the models (i.e. numbers 5–6) that correspond to the upper limit of ∼50 per cent spot filling fraction determined in active G and K stars (O'Neal et al. 1998, 2004). Consideration of such heavily spotted stars is of interest, especially if the contrast ratio, T_p/T_s, is also much lower in mid-late M stars (see Sections 2.1 and 3.3). This may enable more active stars to be included in RV surveys. The six models listed in Table 1 are shown in Fig. 1.

Figure 1

The distribution of spots for models 1 (top left) to 6 (bottom right) for an M dwarf. Models 1 and 2 are analogous to solar min and solar max activity levels while model 3 represents a high solar activity case. Models 4 to 6 are included for completeness (see Table 1).

3 GENERATING SYNTHETIC RADIAL VELOCITY DATA

The imaging code, dots, enables us to generate line profiles using a 3D stellar model with arbitrarily placed spots. RV measurements are then made directly from the line profiles. dots works with two intensity profiles for each image pixel, representing the photospheric and spot temperatures. The degree of spot filling for each pixel is represented by a value in the range 0.0–1.0, where 1.0 represents complete filling (i.e. an image pixel at the spot temperature).

Relative fluxes for a star with a given temperature are determined from the absolute magnitude and radii determined by Baraffe et al. (1998) and Chabrier et al. (2000) for low-mass stars and substellar objects. Hence we are able to determine the fluxes for a star of any mass with a specified model photospheric temperature and an attributed spot temperature. We use two specific photosphere/spot temperature contrast extremes in later sections of the paper. However, in the following section we begin by examining different photosphere/spot contrast ratios. Throughout this paper, we use non-linear limb-darkening values from Claret (2000) to represent the radial variation of intensity for both photospheric and spot temperatures.

3.1 Choice of near-infrared passband

Our choice of infrared passband for the simulations was determined by the observation that some near-infrared wavelength regions contain a greater degree of spectral line information, in the sense that more sharp features are found in these regions. This was noted for example by Reiners et al. (2010) who found that the highest RV precision was achieved in the Y band when compared with the J and H bands. The Y band outperformed the J and H bands in terms of achievable spectral precision for all M dwarf spectral types (see their fig. 5). Only for the early to mid M spectral types, owing to the much higher spectral information content, does the V band still outperform the infrared bands despite much lower flux. As we shall see in the next section, the main advantage of observing at infrared wavelengths is the reduction in starspot-induced jitter. For the following simulations, we therefore focus on the Y band which appears to offer the best chances of detecting low-amplitude signatures that may arise from the reflex motion of orbiting planets.

Figure 5

Detection FAPs versus number of observation epochs for habitable zone planets orbiting M = 0.1 M_⊙ (squares connected by solid/red lines) 0.2 M_⊙ (circles connected by dashed/green lines) and 0.5 M_⊙ (triangles connected by dotted/blue lines) stars. The five curves for each stellar mass represent planetary masses of 20, 10, 5, 2 and 1 decreasing in a left-to-right sense (labelled for the 0.5 M_⊙ curves). The horizontal line indicates the 1 per cent FAP (FAP = 0.01) with the grey region representing undetected planets. All points with FAP < 0.01 are considered as detections of the planet. The plots are for activity model 2 with v sin i = 2 km s⁻¹, instrumental precision = 1.5 m s⁻¹ and with starspot contrast ratios of T_s = 0.65T_p (high contrast) on the left and T_s=T_p− 200 K (low contrast) on the right.

3.2 Relative visible and near-infrared spot induced amplitudes

We first investigate the effect of contrast variations on the amplitude of the RVs. Observations made at two different wavelengths, the V- and Y-band regions (centred at 5450 and 10 350 Å, respectively), are simulated. We begin by placing a single spot with radius 10° on the equator of a star inclined at 90°. The stellar line profile is calculated at a number of phases for which the starspot is visible. The only differences between phases are the intensity (due to foreshortening angle) and RV of the spot contribution to the integrated line profile. Since we are interested in the relative amplitude of starspot jitter only in this section, the simulated v sin i is not important (we used v sin i = 10 km s⁻¹), providing it is greater than the instrumental resolution.

Using the above parameters, we show in Fig. 2 (top), the amplitude induced as a function of wavelength for the case where T_p = 3250 K. The curves represent high-contrast (T_s = 2000 K) and low-contrast (T_s = 3000 K) scenarios. The decrease in amplitude as a function of increasing wavelength is more pronounced for the low-contrast case. Here, a secondary effect, namely the relative equivalent width of the lines at the photospheric and spot temperatures, also becomes important and is responsible for the increase in amplitude at λ = 10 350 and 13 000 Å when compared with shorter and longer wavelengths. This effect is discussed further below. The results are in broad agreement with those presented in Reiners et al. (2010).

Figure 2

Top: RV amplitude as a function of wavelength for an equatorial spot of radius 10°. The curves are shown for high- and low-contrast scenarios with T_p = 3250 K. Bottom: the RV amplitude induced by spots in the visible (V band) compared with the infrared (Y band). For the labelled photospheric temperatures in the range 2000 ≤T_p≤ 4000 K, spot amplitudes ratios are calculated for 1500 ≤T_s≤T_p− 250 K.

In Fig. 2 (bottom) we show the relative V-band/Y-band spot induced amplitudes, K_V/K_Y, for various spot contrast ratios. There are a number of factors that characterize the magnitude of the starspot RV-induced jitter in a given passband, or specific wavelength range, the most important being (1) the contrast ratio between starspot and photosphere and (2) the normalized equivalent width of the local intensity profile. Additional, secondary effects come from changes in limb-darkening coefficient which we also model using non-linear limb-darkening values (see fig. 3 of Claret 2000). The local intensity profile may be considered as the local unbroadened spectrum of a non-rotating star. The combination of the Doppler shifted local intensity profiles from all points on the star, or in our case, all pixels on the model star, are summed up to give the observed profile [it is this summation of local intensity profiles that is also employed to derive indirectly resolved Doppler images (Vogt & Penrod 1983; Vogt, Penrod & Hatzes 1987; Collier Cameron 2001) of spotted stars]. The effects of (1) and (2) are investigated and discussed in detail in Reiners et al. (2010). Here, we summarize those effects in order to illustrate the trends seen in Fig. 2. The Doppler shifted local intensity corresponding to the location of a cool spot on a star is much less compared with the equivalent photospheric contribution. The corresponding photosphere/spot intensity ratio is typically greater than an order of magnitude for standard solar-like spots that exhibit T_p−T_s∼ 1000 K. In this instance, because the spot contribution to the Doppler broadened profile is so much smaller than the photospheric contribution, the spot profile EW and shape are relatively insignificant. In other words, spots contribute little flux, both in the line and the continuum, when the contrast between photosphere and spots is large. Hence the missing light at a given Doppler shifted position is relatively independent of wavelength. This is illustrated in Fig. 2 (bottom) by the observation that for any T_p, K_V/K_Y is close to unity for the smallest values of T_s. The value of K_V/K_Y increases as the spot temperature (and hence intensity) increases (i.e. as T_p/T_s decreases). This trend arises simply because the contrast ratio for fixed T_p/T_s is much lower in the infrared than in the optical. The wavelength dependence of the RV signature is illustrated for blackbody fluxes in Reiners et al. (2010). Here we use fluxes derived from the models of Baraffe et al. (1998) and Chabrier et al. (2000) for 5-Gyr stars. The top-right point on each curve in Fig. 2 (bottom) indicates T_p−T_s = 250 K. This small temperature difference leads to starspot-induced jitter that is several times smaller in the Y band as compared with the V band.

Figure 3

The Y-band rms RV variation for each of the six simulated spot models. Pairs of values (connected by a line) are shown for each stellar mass (0.1, 0.2 and 0.5 M_⊙ cases are labelled for model 2) and stellar v sin i values of 1, 2, 5, 10, 20 and 50 km s⁻¹. Each pair of values represents the upper and lower variation for the two extreme F_p/F_s contrasts listed in Table 2.

For the T_p = 3000 and 4000 K cases, reductions in starspot-induced jitter of an order of magnitude are achieved. These augmented reductions are due largely to the behaviour of the relative line intensities in the spot and photosphere at V- and Y-band wavelengths. While the continuum contrast is low, the relative equivalent width and depth of the lines is greater in the V band, compared with the Y band, leading to the large observed ratios in Fig. 2. Again, this effect is discussed by Reiners et al. (2010) in section 4.1.2 of their paper. In essence, for low contrast ratios, the spot-induced jitter is much smaller since there is little difference between the intensity of the photospheric and spot contributions. However, the ratio of line depths at T_p and T_s then becomes important. If the continuum contrast ratio (1) is small but the line depth of the spectral lines in the spot is greater than in the photosphere, the two effects may cancel. In other words, the emission bump created by the cooler spot spectrum is cancelled by additional absorption due to deeper lines (see Reiners et al. 2010, fig. 11). We derive a mean line depth for each passband region from the models of Brott & Hauschildt (2005). For T_p = 3250 K, we see in Fig. 2 (bottom) that K_V/K_Y is much smaller at low T_p/T_s. This arises because of a significant growth in equivalent width of lines in the Y band below this temperature. For the cases where the continuum contrast ratio is small, the centre of the lines are then much deeper, leading to higher contrast in the lines. As we shall see in Section 3.3, this leads to larger starspot-induced velocity amplitudes for our simulated 0.2 M_⊙ star in low-contrast scenarios.

Our results in Fig. 2 are in close agreement with those presented in Reiners et al. (2010). A marked difference between the change in RV-induced jitter as a function of wavelength is seen when moving from their ‘toy’ model results of fig. 11 to the model atmosphere results (analogous to our results) of fig. 12. When comparing K_V/K_Y, Reiners et al. (2010) find an approximate decrease in the starspot-induced jitter of ∼two and ∼five times for T_p = 3700 and 2800 K, respectively, with T_p−T_s = 200 K. We similarly find starspot-induced jitter ratios of ∼5 and ∼6 for T_p = 3750 and 2750 K, respectively, with T_p−T_s = 250 K. The discrepancies in the results likely arise from differences in the models, the line depths used to represent the spot and photosphere local intensities and to a smaller degree the differences in limb-darkening models used. Nevertheless, the results are in broadly close agreement, indicating a reduction in RV-induced jitter in the Y band over the V band of order five for T_p−T_s = 250 K.

Fig. 2 shows a clear advantage of observing at near-infrared wavelengths, when the contrast ratio between spots and photosphere is not large. At higher contrasts, the advantage is less obvious. Obtaining an estimate for reasonable values of T_s versus T_p is not straightforward, but clearly has important implications for further simulations. This issue is considered further in the following section.

3.3 Radial velocities induced by random spot distributions

The absolute RV amplitudes depend not only on the starspot distributions but also on the photosphere/spot contrast ratio. Following Reiners et al. (2010), we have investigated two extreme scenarios for M dwarf stars with three different masses. In Table 2, we list the chief parameters for the simulated stars, including the photospheric and spot temperatures.

Stars with axial inclination i = 90° are simulated according to the parameters tabulated in Table 2 and for each of the six stellar models described above. A total of 36 line profiles (i.e. every 10°) with a range of v sin i values of 1, 2, 5, 10, 20 and 50 km s⁻¹ are generated for a complete stellar rotation period. The local intensity profile used to represent the rotation profile of a non-rotating star is a synthetic Voigt profile with an instrumental resolution that is appropriate (R = 70 000) for high-resolution IR spectroscopy (e.g. Jones et al. 2009).

Fig. 3 plots the starspot-induced rms RV variation for each spot model (introduced in Section 2.2) and v sin i value. The pairs of values for each stellar mass of 0.1, 0.2 and 0.5 (points connected by lines in Fig. 3) represent the (greater rms RV) and K (lesser rms RV). As expected, there is little variation due to the stellar mass for the point since the relative contrast ratio is unchanged. The lower point in each case shows more variation for a given model since the ratio K is not fixed. Moreover, for the T_p/T_s combination of 3250 K/3050 K, we see that the RV variation is between 1.5 and 6 times larger than for the corresponding 3750 K/3550 K case. This difference is due to the effects discussed in the Section 3.2, where the equivalent width of the lines grows for temperatures close to, but below 3250 K. The large range in factors of between 1.5 and 6 is likely due to the spot sizes and number that define each model and their relationship with the Doppler/spatial resolution limit which improves (i.e. smaller spots resolved) with increasing v sin i. This will ensure that for low-contrast spots on our 0.2 M_⊙ model, the starspot jitter will be greater than for the 0.1 and 0.5 M_⊙ models.

For model 1, with v sin i = 2 km s (i.e. a slow rotator analogous to the Sun) the approximate range of velocities (considering all three stellar masses) for the two T_p/T_s extremes are ∼2–6 cm s⁻¹, while for model 2 the range of observed values is ∼4 cm s⁻¹–27 m s⁻¹. This is an order of magnitude smaller than the values observed for the Sun by Lagrange et al. (2010) which show equivalent RV jitter of solar minimum and maximum periods of up to 0.6 and 2 m s⁻¹, respectively. The reduction in contrast ratio in the optical versus the infrared, and the fact that spots may appear at all latitudes, act to reduce the simulated RV jitter for models 1 and 2. This latter point is important since placing only a few small spots randomly at all latitudes, as opposed to restricting them to solar-like low-latitude bands may lead to significant reductions in the measured jitter (i.e. spots at low latitude induce the greatest jitter). By the time a star reaches v sin i = 10 km s⁻¹, we find rms velocities of ∼4–38 cm s⁻¹ and ∼17 cm s⁻¹–1.9 m s⁻¹ for models 1 and 2, respectively. Clearly, the solar minimum and solar maximum analogue cases (models 1 and 2) will enable the greatest possibility of detecting lower mass planets since they yield the lowest rms RV variations. As expected, the rms jitter peaks at model 6 which corresponds to a star with a 62.4 per cent starspot filling fraction.

For those early M dwarfs that are fast rotators, the spot sizes derived from Doppler images (Barnes & Collier Cameron 2001; Barnes et al. 2004) are much larger than is seen for the Sun. It is not clear however whether the spots seen in Doppler images are actual individual spots or indeed unresolved spot groups. Indeed if we take model 6 and reconstruct the surface spot distribution for a moderately rotating case, we derive images resembling the Doppler images of Barnes & Collier Cameron (2001) and Barnes et al. (2004). We have also investigated the effect of larger spots on the RVs by scaling the sizes used for models 1–6 by factors of 2, 5 and 10. In other words, a model identical to each of those shown in Fig. 1 was created, but with all spots scaled up in radius. The starspot-induced rms jitter maximum in Fig. 3 is found to increase by a factor of ∼2 for a spot size scalefactor of 10. The main effect however is a change in the model number at which the peak jitter occurs. For scalefactor 2, the rms jitter is less peaked for model 6 when compared with model 5. For scalefactor 5, a peak occurs between models 3 and 4, while the peak is at model 3 by scalefactor 10. However, even for the more modest scalefactor 2 scenario, the planet detection thresholds (see Section 4) for model 2 will more closely resemble those of model 3 where the spot scalefactor was unaltered. In the following sections, we do not scale the spot sizes, but use the spot size distributions as previously defined in Section 2.2.

3.4 Stellar mass versus v sin i and spot filling

The effect of v sin i is an important factor that must be considered when determining the RV jitter. While Fig. 3 shows that model 1 does not yield jitter of >1 m s⁻¹, even at v sin i = 20 km s⁻¹, a very active star such as that represented by model 6 only yields jitter of order 1 m s⁻¹ or less for the low-contrast case when v sin i = 1 km s⁻¹. There is approximately a 1:1 correspondence between jitter and v sin i for models 3–6. In other words, doubling the v sin i doubles the rms jitter. However, for models 1 and 2, the increment is small owing to the finite spectral resolution of 70 000 which corresponds to 4.3 km s⁻¹.

Since the rotation–activity relation seen at earlier spectral types is also seen to persist in M dwarfs, through observation of Hα (Mohanty & Basri 2003), it seems reasonable to assume that the slower rotators will be the least active stars while the more rapid rotators exhibit higher degrees of spot coverage. Unfortunately, since there have been no successful studies of spottedness of M dwarfs (O'Neal et al. 2005), relating v sin i to the spot filling fraction is not possible. We must therefore make assumptions about the degree of spottedness as a function of rotation velocity. Based on solar observations and on the spot filling fraction measured on other stars (O'Neal et al. 1996, 1998, 2004), we assume that models 1 and 2 are appropriate for v sin i = 2 km s⁻¹, and models 3–6 are appropriate for v sin i = 5–50 km s⁻¹.

Jenkins et al. (2009) have recently determined the v sin i values for a number of M dwarfs, thereby extending the sample of measured rotation values. This study reinforced the observation (Delfosse et al. 1998; Mohanty & Basri 2003; Reiners & Basri 2008) that amongst field M dwarfs, the earlier spectral types exhibit slower rotation compared with mid and late M spectral types. In other words, while earlier spectral types are more likely to be slower rotators, enabling more precise RVs to be determined, the higher mass of the star sets a lower limit on the reflex motion due to an orbiting planet of given mass. Conversely, the detectable mass limit for an orbiting planet is adversely affected by the more moderate rotation of mid to late M dwarf stars, while the lower mass is more favourable for detecting low-mass planets.

For spectral types M2V and M4V (0.5 and 0.2 M_⊙, respectively), Jenkins et al. (2009) (see their fig. 9) find the median v sin i is 3–5 km s⁻¹ whereas for spectral types M6V (0.1 M_⊙), v sin i≃ 8 km s⁻¹. However, velocities from 1 to several 10s of km s⁻¹ are found for all stars in the spectral range considered, albeit with fewer slow rotators found at later spectral types [i.e. for masses ≤0.1 M_⊙ (M6V)]. We therefore simulate a range of rotation velocities for each spectral type, but urge the reader to bear in mind the trend of increasing mean v sin i with decreasing stellar mass.

4 DETECTION THRESHOLDS FOR M DWARF STARS WITH STARSPOTS

In this section, we simulate detection thresholds for low-mass planets orbiting at the centre of the habitable zone of M dwarf stars. We use the starspot models from the preceding section to derive activity induced jitter. Since the parameter space for characterizing detection thresholds is so large, we have chosen to fix a number of parameters as in previous sections. We have simulated the detection thresholds for 1, 2, 5, 10 and 20 planets that orbit in the habitable zones of 0.1, 0.2 and 0.5 M_⊙ stars. Table 2 lists the periods of the planets for each stellar mass and the corresponding stellar RV amplitudes, K_*, for each planet. In this section we have simulated all stellar and planetary orbit inclinations with i = 90°, and used a range of v sin i values as discussed in Section 3.4.

For each star/planet combination, we generate RV points for a range of observation epochs. We simulate 10, 20, 50, 100, 200 and 500 epochs, with one observation made every night for simplicity. Jitter from the two sources (i.e. starspots and instrumental/measurement precision) is then added to the planetary RVs. The starspot models 1–6 are used to add the stellar activity jitter to each planetary RV point. This is achieved by sampling, at a random observation phase, the line profile that arises from the particular starspot model we are interested in. We use model estimates to approximate the instrumental/measurement precision. This is an important further consideration, since v sin i affects the best precision that can be achieved. For an assumed Y band signal-to-noise (S/N) ratio of 90–100 and R∼ 70 000, fig. 7 of Reiners et al. (2010) indicates the appropriate accuracy that may be achieved in the Y band for a 3000-K atmosphere. This varies from ∼1.5 m s⁻¹ at v sin i = 2 km s⁻¹ to ∼11 m s⁻¹ at v sin i = 50 km s⁻¹. An example of a simulated planet induced stellar RV signal is shown in Fig. 4 before (solid/red) and after (dashed/green) jitter is added. This planet is detected with FAP < 0.01 in Fig. 7 (upper right-hand panel, green square at 100 epochs).

Figure 4

Example of a simulated stellar RV curve for a 2 planet orbiting a 0.2 M_⊙ star (red/solid line). The RV curve is shown with added jitter for an instrumental resolution (IP) of 2 km s⁻¹, stellar v sin i = 5 km s⁻¹ and starspot model 2.

Figure 7

As for Fig. 5 for activity models 3 (top) and 6 (bottom) with v sin i = 10 km s⁻¹ (instrumental precision = 3 km s⁻¹).

We then carry out a Lomb Scargle periodogram analysis (Press et al. 1992) on each RV curve in an attempt to recover the planetary signature. Although this is adequate in our simulations, we note that we have not included the effects of orbital eccentricity (i.e. e = 0 in all simulations), which will modify our detection thresholds (Ford 2005) in extreme cases. The results of our periodogram analyses are plotted in Figs 5–8 as detection thresholds. We have selected to illustrate a range of scenarios that represent varying v sin i and activity levels, following the preceding discussion.

Figure 6

As for Fig. 5 for activity models 2 (top) and 4 (bottom) with v sin i = 5 km s⁻¹ (instrumental precision = 2 km s⁻¹).

Figure 8

As for Fig. 5 for activity models 3 (top) and 6 (bottom) with v sin i = 20 km s⁻¹ (instrumental precision = 6 km s⁻¹).

The false alarm probability (FAP) is plotted against the total number of observations (see caption of Fig. 5 for full details). The left-hand plots show the results for the cases where T_s=T_p− 200 K (low contrast) while the right-hand plots show the cases for which T_s/T_p = 0.65 (high contrast). As expected, our ability to detect a planetary signature decreases with increasing v sin i and stellar activity level.

Fig. 5 represents our low v sin i and activity (model 2 – active solar analogue) level. Since both models 1 and 2 exhibit significantly <1 m s⁻¹ starspot jitter (see Fig. 3) with v sin i = 2 km s⁻¹, we illustrate the results for model 2 only since Reiners et al. (2010) predict that the highest precision achievable in the Y band with v sin i = 2 km s⁻¹ is ∼1.5 m s⁻¹. The precision therefore dominates the noise in this scenario rather than the starspot jitter. As few as 20–30 epochs of observations are required to detect ≤2 planet orbiting a 0.1 M_⊙ star, while 50 epochs are required to detect a 1 planet. However for a 0.5 M_⊙ star, 500 epochs are required to detect a 1 planet. Obtaining such a large number of observations is unlikely to be feasible, at least with limited telescope allocations and a sufficiently large target sample. There is little difference between the low- and high-contrast scenarios, again because all starspot-induced jitter is well below the achievable precision level.

Fig. 6 illustrates results for v sin i = 5 km s⁻¹ and two different activity models (models 2 and 4 – with up to 9 per cent spot filling). In this instance, for model 2, the changes are slight compared with the preceding case where v sin i = 2 km s⁻¹. This observation illustrates a further important limitation to precision, namely the instrumental resolution. A greater degree of spottedness is required (i.e. model 4 cases in Fig. 6) before a noticeable change in detection thresholds begins to take place. Here, while 80 epochs in the low-contrast case will enable 1 planets to be detected, ∼500 epochs are required to detect 1 planets in the high-contrast scenario. Jenkins et al. (2009) report a median v sin i for 0.5 and 0.2 M_⊙ stars that most closely match the v sin i = 5 km s⁻¹ used in Fig. 6. This simulation may therefore be deemed to represent the detectability of earth-mass to few-earth-mass planets around average 0.5 and 0.2 M_⊙ stars.

In Fig. 7, we show the expected detectability of planets around stars with v sin i = 10 km s⁻¹. This rotation velocity is the median for later M spectral types, such as the 0.1 M_⊙ star. By 10 km s⁻¹, the less active model 3 scenario still enables 2 planets orbiting 0.1 M_⊙ stars to be detected in the low-contrast case, whereas 5 planets may be detected in the same number of epochs for the high-contrast case. Only planets with mass ≥10 may be detected with <100 epochs of observations for the high-contrast case. With a highly active star (model 6), only one planet (20 ⁠) remains detectable in the high-contrast regime, albeit requiring 500 epochs of observations. For model 6, the starspot jitter contribution is so large (compared with the instrumental/measurement precision) that it dominates the noise contribution. We also noted in Section 3.3 (see also Fig. 3) that the jitter is greatest for the 0.2 M_⊙ model with low spot contrast spots. As a result, in this scenario, our model predicts that detection of a planet orbiting a 0.5 M_⊙ star becomes easier than for a planet orbiting a 0.2 M_⊙ star.

The global trend seen in Figs 5–7 continues in Fig. 8, where only the most massive planets may be detected for very active stars. If the spot filling factors of up to 50 per cent, reported at earlier spectral types for stars rotating with ∼20 km s⁻¹ (O'Neal et al. 1998) are to be found in M dwarfs, then a lower planet detection limit of order one Neptune mass may be expected. By 50 km s⁻¹ (not shown), only 20 planets may be detected in orbit around a 0.1 M_⊙ star with low-contrast spots (model 6) in less than 100 epochs of observations (for model 4, only ≥10 planets may similarly be detected in ≤ 100 epochs).

4.1 The effects of activity and observing strategy on planet detectability

We have demonstrated how the relative contrast ratio F_p/F_s is an important factor that determines the detectability of planets, particularly when spot activity dominates the jitter. However, it must also be realized that the observing strategy simulated in this paper will determine the detectability of planets. In particular, the orbital period of the planet may be either somewhat shorter, or longer than the span of the observations. Table 2 indicates that the orbital periods of the simulated habitable zone planets are ∼5, 13 and 36 d for 0.1, 0.2 and 0.5 M_⊙ stars, respectively. Hence, detection of a periodic signal on time-scales shorter than the period is less likely. In other words, it should be easier to detect habitable zone planets orbiting lower mass stars owing to the shorter periods. Since we only simulate a minimum of 10 epochs, this mostly applies only to the 36 d period experienced by a planet orbiting a 0.5 M_⊙ star. It can be seen from Figs 5–8 that no planets are detected in orbit about a 0.5 M_⊙ star until over half the period has been sampled. Therefore if targeting such stars for habitable zone planets, a different observing strategy would be needed to minimize the number of observation epochs. In reality, a RV survey would wish to search for planets with a range of orbital radii and periods, especially as shorter period planets are more likely to be detected.

The effect of activity on the number of epochs required to detect a planet can also be assessed by completely removing the starspot jitter from all simulations. In our low activity simulation (v sin i = 2 km s⁻¹ instrumental precision of 1.5 m s⁻¹ and model 2 starspot coverage), removing the starspot activity does not have any significant effect on the detection curves plotted in Fig. 5. This has already been discussed in the preceding section and is a consequence of the instrumental precision at this rotation velocity (and implied instrumental precision) dominating the jitter. The effect of removing spot activity from the simulations presented in Figs 6–8 however leads to significant changes in the number of epochs necessary for detection, as might be expected. In our most extreme scenario simulated in Fig. 8 (high contrast), although no planets are detected, complete removal of the starspot jitter leaves a 6 m s⁻¹ precision floor due to the rotation velocity of v sin i = 20 km s⁻¹. Here, planets of 5–20 may then be detected with 300–40 epochs, respectively.

5 SUMMARY AND CONCLUSIONS

We have used more realistic starspot models than previous studies to determine the effect of activity-induced jitter in precision RV studies of M dwarf stars. We have demonstrated that with only several tens of epochs, habitable zone earth-mass planets can be detected around low-activity stars. Since the contrast ratio between photosphere and spots is uncertain, we opted to simulate two extreme cases. In reality, the true contrast ratio likely falls between the two extreme cases simulated. Rockenfeller et al. (2006a) required temperature differences of only a few hundred K in order to fit their light curves of M5V and M9V stars. However, their modelling assumed only a single spot was present. A more uniformly spotted surface, especially if highly spotted, could produce a similar light curve but would be expected to require a higher contrast between photosphere and spot to achieve the same photometric amplitude.

The factors that determine the starspot-induced jitter have important consequences for the estimated detection threshold limits. In agreement with other studies (Reiners et al. 2010) we showed (Fig. 1) that once the contrast ratio between spot and photosphere is sufficiently high, the improvement in jitter, when moving from V band to Y band, is less impressive. In this case, the spots are simply dark enough that the radial jitter does not show such a strong decrease with increasing wavelength. Conversely, the decrease in jitter as a function of increasing wavelength is much more pronounced at lower contrast ratios. At low contrast, the exact normalized line strengths and relative line strengths of the photosphere and spot play an important role in determining the exact relative jitter at different wavelengths. This additional factor leads to the irregularity in relative RV amplitudes as a function of increasing T_p, as shown in Fig. 2 (see Section 3.2). Although we do not show the full jitter amplitude ratio for other photometric bands, Fig. 2 (top) for T_p = 3250 K indicates that the relationship between wavelength and jitter is not as straightforward as a simple blackbody model might suggest. Further estimates of the wavelength dependence of jitter may be obtained by observation of the results of Reiners et al. (2010) who plot RV amplitude as a function of wavelength for the 5000–18 000 Å range (their fig. 12). At longer wavelengths than Y band, equivalent width effects at low photosphere/spot contrast apart, there is little further gain in precision. A further important consideration is the number of lines available for cross-correlation, which may be fewer at longer wavelengths.

We have assumed that activity scales with rotation (Browning et al. 2010) so that less spotted stars are slow rotators and more spotted stars are fast rotators. However, as discussed in Section 2.1, Reiners & Basri (2010) find that the relation between rotation and activity is weaker in M dwarfs. Knowledge of the true spottedness of these stars is clearly vital for an accurate estimation of the detectability limits for low-mass planets orbiting the lowest mass M dwarfs. In addition, spot coverage derived from other spectroscopic methods using temperature sensitive (TiO) lines (O'Neal et al. 1998) indicate that there may be large discrepancies with results from Doppler images. For active stars, up to ∼50 per cent spot coverage has been derived using the TiO method, while Doppler images of similar stars typically derive <10 per cent coverage. The difference between these methods is that the TiO procedure derives an unresolved mean spot coverage while the Doppler imaging results are limited by the amount of resolvable information in the broadened rotation profile of the star. It seems likely that active stars are therefore more spotted than Doppler images show, and that they possibly exhibit smaller spots than are derived via this method. In summary, high spectral resolution Doppler imaging surveys (in the red optical or infrared), possibly combined with other global spot coverage surveys are needed for a more accurate picture of spot patterns on mid-late M dwarfs.

It is clear that there are many physical factors that will determine the true detection threshold for planets orbiting M dwarfs; we have only simulated planets in circular orbits for instance. Kennedy & Kenyon (2008) simulated planet scattering and found that for low-mass stars only, planets with long circularization times on eccentric orbits could form. Although we have not simulated the effects of eccentricity in full, we have already noted that in extreme cases, detection thresholds will be raised and more observations will be required (Ford 2005). We find that for longitude of periastron, ω = 0°, v sin i = 10 km s⁻¹ (instrumental precision = 3 m s⁻¹), the number of nights required to make a detection increases by approximately 2–2.5 times for eccentricity e = 0.5 and by ∼5 times for e = 0.9.

In this paper, we have shown that earth-mass or near-earth-mass planets in the habitable zones of late to mid M dwarfs can be detected when the number of starspots matches those seen on the Sun at extremes of activity. Moderate rotation quickly increases the number of observations that are required to make a detection such that by the time v sin i = 10 km s⁻¹, ≥100 observations are required to detect planets with masses ≤5 ⁠. Variations in v sin i and contrast ratio are therefore important factors in particular in determining the detection thresholds. Given the evidence for a general increase in v sin i with decreasing mass (Delfosse et al. 1998; Mohanty & Basri 2003; Jenkins et al. 2009) among M dwarfs, the balance of these factors may be important. Obtaining estimates of starspot coverage from further modelling and photometric observations will provide vital information that will enable a link between starspots and more traditional chromospheric activity indicators (Delfosse et al. 1998; Mohanty & Basri 2003; Reiners, Basri & Browning 2009; Browning et al. 2010) to be made.

1

http://www.exoplanet.eu

We would like to thank David Parker (Monk's Walk School, Herts, UK) and Robert Sanders (Richard Hale School, Herts, UK) for their contributions to this work. We would also like to thank the referee, Ansgar Reiners, for constructive comments that have resulted in an improved final version of the manuscript. SVJ currently acknowledges support from De Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO).