© ESO, 2015

1. Introduction

Thanks to the radial velocity (RV) technique, more than 500 exoplanets have been discovered in the two last decades (http://exoplanets.eu). First limited to Jupiter-like planets on short period orbits, today the RV technique allows the detection of Neptune-like planets (10−40 M_Earth) and Super Earths (1.2−10 M_Earth) on longer period orbits. This progress was made possible thanks to the important improvements in the sensitivity and stability of instruments, as well as to observational strategies to average out known sources of stellar signals (Dumusque et al. 2011a). For example, the High-Accuracy Radial velocity Planet Searcher (HARPS) spectrograph located on the ESO 3.6 m telescope in La Silla (Pepe et al. 2002) gives a precision of 1 m s^-1 for a solar-type star in average conditions, compared to 5−10 m s^-1 with previous instruments. Furthermore, future instruments, such as ESPRESSO on the Very Large Telescope (VLT) or G-CLEF on the Giant Magellan Telescope (GMT) are expected to reach a precision down to 0.1 m s^-1 (Mégevand et al. 2010; Fżrész et al. 2014). This level of precision will theoretically give access to lower planetary masses far from their host stars such as Earth-like planets in the habitable zone (HZ).

However, low-amplitude RV planet signals (such as RV signatures of Earth-mass planets in the HZ) are much more sensitive to stellar perturbations than giant planet RV signatures. In the case of a low-mass planet RV signature with an amplitude in the 0.1−1 m s^-1 range, the stellar noise or “jitter” is high enough to either mask or mimic the planet-induced RV variations, even in the case of a chromospherically nearly inactive star (Isaacson & Fischer 2010; Lagrange et al. 2010; Meunier et al. 2010).

The so-called stellar jitter mainly comes from three different sources: stellar oscillations or pulsations, granulation, and stellar magnetic activity. On short timescales, stellar pulsations or oscillations dominate the RV jitter. For FGK dwarfs, they are mainly driven by accoustic or pressure waves (p-modes). Pressure waves are commonly believed to originate from turbulent convective motions occurring in the stellar outer layers and to propagate through the star. They induce RV shifts with an amplitude of a few cm s^-1 to one m s^-1. The oscillation periods range from a few minutes (e.g. five minutes in the case of the Sun) to a few tens of minutes. The amplitude and the period of the oscillations increase with stellar mass. Long exposure times and observational strategies have shown to be efficient in averaging the RV noise induced by pulsations in the case of late-type stars (Santos et al. 2004; Dumusque et al. 2011b).

The photospheric granulation phenomenon accounts for the convective plasma motions occuring in the outer envelope of solar-mass stars. The quick rise and fall of the plasma results in a pattern of bright granules and darker lanes at the stellar surface, and in rapidly evolving RV shifts, with amplitudes up to a few km s^-1 locally. When integrated over the entire stellar disc, upflows and downflows average out, leaving a residual RV jitter at the level of the m s^-1. So far, there have been few attempts to estimate the impact of granulation on the RV. Dumusque et al. (2011b) made a first estimation for different spectral type stars based on the study of asteroseismology measurements. Cegla et al. (2013) developed a four-component model of granulation, building absorption line profiles from 3D magnetohydrodynamic solar simulations. This study was continued by Cegla et al. (2015), who derived the corresponding RV time series and other time series of observables such as the bisector amplitude. They concluded first that granulation has a strong impact on RV at the cm s^-1 to m s^-1 levels and would be a potentially significant limitation to low-mass planets detection. Second, they found correlations between the RV time series and other observables (e.g. bisector curvature or bisector inverse slope) that already allow a partial correction of the granulation signal and reduce the corresponding jitter by one third. Meunier et al. (2015) simulated a collection of both granules and supergranules on the solar visible hemisphere over a whole solar cycle and derived the induced RV and photometric time series. These authors concluded that the granulation noise is hard to average out even over an entire night and will therefore have a significant impact on detection limits.

On timescales from a few days to several years, the stellar jitter is dominated by the so-called stellar magnetic activity. Active regions such as dark spots and bright faculae (colder and hotter than the quiet photosphere, respectively) grow and evolve on the stellar surface, inducing changes in the stellar irradiance and RV variations (because the flux loss or excess in the active regions distorts the CCF of the stellar spectra). The stellar rotation causes these active regions to be seen as moving across the stellar disc, thus inducing an apparent Doppler shift. The first studies of the impact of starspots on RV and line bisectors were made by Saar & Donahue (1997) and Hatzes (2002). Saar & Donahue (1997) built a semi-empirical law that directly bound the spot-induced RV semi-amplitude, the spot filling factor and the stellar rotational velocity for cool F-G stars. Saar (2003, 2009) made similar studies for bright faculae. The first spot model was developed by Desort et al. (2007). The authors computed synthetic stellar spectra and applied a black-body law to take into account the contribution of a colder spot. They accurately quantified the RV amplitude, RV bisector and photometric variations induced by a starspot for a range of different stellar and spot properties.

Such simple spot models have proved to be useful and efficient when simulating the RV signature of one or a few active regions during a timescale of the order of the stellar rotational period. They are nicely adapted to stars that host a single main spot or a few spots, such as active young dwarf stars (e.g. HD 189733, see Dumusque 2014, and references therein). They are, however, not able to reproduce the total activity-induced RV signal for Sun-like stars owing to the presence of multiple evolving active regions on the stellar disc, and owing to the strong contribution of faculae because of slow rotation (Meunier et al. 2010; Dumusque et al. 2014). A complete model of the stellar activity pattern for longer timescales (approximately the activity cycles) is needed. For most stars though, we have no detailed information on the activity properties as active structures are not directly observable. Doppler and Zeeman-Doppler imaging technics allow the largest active structures (or active structure clusters) to be recovered only for fast-rotating stars (i.e. young, active ones) but not for old, slow-rotating solar-type stars (Rice 2002). The structure and convection properties are mostly unknown and we rely only on indirect and global estimators such as the calcium (Ca) index (see e.g. Baliunas et al. 1983; Noyes et al. 1984; Hall & Lockwood 2004).

The Sun is an exceptional star for which we have a lot more information. Indeed, both solar dark spots and bright features have been amply observed and the various solar activity properties are thus profusely described. If seen as a moderately active star, then the Sun is an ideal prototype for studying how activity affects the detectability of low-mass planets. In a series of papers, we modelled with a great accuracy the solar activity pattern using detailed observations (i.e., spots catalogues and magnetograms) extending over the full solar cycle 23. We then rebuilt the induced photometric variations and the corresponding spectra, and deduced from the latter the RV time series. We first considered only cold spots in Lagrange et al. (2010; hereafter Paper I), studying their impact on the detectability of an Earth-mass planet in the HZ of a solar-like star. In Meunier et al. (2010; hereafter Paper II), we extended this study by considering in addition the contribution of bright faculae and, for the RV, the attenuation of the convective blueshift in magnetically active regions. Combining the effect of the three activity components, we found that the attenuation of convective blueshift in faculae (hereafter the convective component) dominated the induced RV variations (this happens because the Sun is a slow rotator: for larger vsini, the photospheric effect of spots becomes dominant). We showed that without correction tools (still to be identified and tested) it would be impossible to detect an Earth-like planet in the HZ. Using the same simulation as in Paper II, Lagrange et al. (2011; hereafter Paper III) estimated the astrometric effect of stellar magnetic activity. Finally, in Meunier & Lagrange (2013; Paper IV), we used the strong correlation between the convective component and the Ca index (as both are directly related to the chromospheric plage or photospheric facula filling factor) as a tool for correcting RV time series from the convection-induced jitter. We showed that an Earth-like planet in the HZ would become detectable given that: i) there is an excellent signal-to-noise ratio (S/N) on the Ca index data; and ii) the temporal coverage of the RV observations over one activity cycle is very dense (typically one night out of four during the cycle).

The objective of the present paper is to extend the previous study to a Sun-like star seen under any inclination. To this end, we fully parameterize the activity pattern of a Sun-like star. A few activity parametrizations have already been developed. Barnes et al. (2011) undertook the first attempt to parameterize dark spot distributions to low- and high-activity levels in the case of M-type stars. A similar simple model was used by Jeffers et al. (2014), this time for young and active G and K dwarfs. A parametrization of an activity pattern was made by Dumusque et al. (2011a), who simulated dark spot distributions over the stellar surface for different levels of activity of a solar-like star. The authors derived detection limits for different activity levels and observational strategies. This study was limited to dark spots, however, and did not include the flux effect of bright features (faculae and network), nor the inhibition of the convective blueshift in active structures.

Using the Sun as a template presents several important advantages: i) we have much more information on solar activity than for any other star; ii) the validity of such a model can easily be asserted by a comparison with the previous papers; and iii) such an activity model can be adapted to other spectral types and other activity levels than the Sun’s, if the active region configuration is similar. In addition, it allows the impact of stellar inclination on the results to be studied, contrary to the solar observations for which only one hemisphere is observed at any given time.

In Sect. 2, we build the spot and facula distributions on the Sun surface over a full solar cycle. We describe in detail the laws and parameters that we used to make these simulations. To test our approach, we then compare the obtained structure distributions with the observed solar distributions used in Paper II. In Sect. 3, we compute the induced RV variations, taking into account both the photometric contribution of spots and faculae and the attenuation of the convective blueshift in active structures. Then we compare our results with those previously obtained using observed solar magnetic features. In Sect. 4, we study the influence of the stellar inclination on the resulting RV and photometry. The astrometric time series will be treated in a separate paper. We also compute the corresponding detection limits. Finally, we discuss our results and the use of our model for other stars in Sect. 5.

2. Building the activity pattern

2.1. Approach

We consider the solar activity pattern during cycle 23. We use the Sun as a template and as a means of comparison with the observed activity pattern used in Papers I and II. We parameterize an extended range of activity levels and timescales. We start from the general distribution of the activity level over the cycle, then we parametrize the spatial and temporal distributions of active structures and their dynamics, and we finally model the individual behaviour of the structures, including the dark spots, the bright faculae around these spots and the network smaller features. The two outputs of our simulations are the lists of the dark spots and of the bright features, respectively. Each file gives, for each time step of the simulation, the structure sizes (in millionths of hemisphere, hereafter μHem), and latitudes and longitudes. For a given set of parameters, these files will then represent the inputs of our simulation tool to produce the spectra and the corresponding RV at each time step, as was done in Papers I and II with observed solar structure lists.

2.2. Building the spot and facula catalogues: input parameters

All input parameters are summarized in Table 1. We detail them here.

• 1.

  Global activity level. The sunspot-induced solar activityfollows a cycle of12 ± 1 yr. To mimic as accurately as possible the global activity level and its distribution over such a cycle, we used the relative sunspot number, or Wolf number (noted R), as a proxy. For this study, we based our simulations on solar cycle 23. This allowed us both to build a realistic activity pattern and to compare our results with those previously obtained in Paper II. The daily Wolf number R for cycle 23 was recovered at the Solar Influences Analysis Data Center¹ (SIDC) and is given in Fig. 1 (upper panel). We first smoothed the daily R data to obtain a long-term reference R_smth (Fig. 1, two top panels) for the cycle shape.

  Fig. 1 Top: daily (dots) and smoothed (bold line) Wolf number over solar cycle 23. Middle: same smoothed Wolf number (black line), and when randomly adding a second-order polynomial dispersion (red dots). Bottom: power spectra of the observed daily (black) and simulated (red) Wolf numbers (respectively R and R′). The solar rotation period (blue solid line) and its first harmonic (blue dashed line) are also shown for comparison.

  The observed daily dispersion of the Wolf number comes from two contributions: the first comes from the rotational modulation and will naturally be present in the model (hence the smoothing); the second occurs because appearing spots do not follow a smooth curve and present some dispersion originating in the dynamo process. To take this contribution into account, we calculated a second-order polynomial dispersion d to the smoothed curve R_smth, expressed as $\mathit{d}\mathrm{=}{\mathit{d}}_{\mathrm{0}}\mathrm{+}{\mathit{d}}_{\mathrm{1}}{\mathit{R}}_{\mathrm{smth}}\mathrm{+}{\mathit{d}}_{\mathrm{2}}{\mathit{R}}_{\mathrm{smth}}^{\mathrm{2}}$(1)where d₀, d₁, and d₂ are derived empirically from observations. We then randomly added this dispersion d to the smoothed data R_smth to obtain a daily input Wolf number R′ (Fig. 1, middle panel). The Wolf number (R, R′) is a combination of the number of individual spots and of the number of spot groups. As the input of our simulation is the total number of spots N (which accounts both for the individual spots and for those in groups), we have derived a calibration for deducing this total number of spots from an activity level expressed in Wolf number units, $\mathit{N}\mathrm{=}{\mathit{n}}_{\mathrm{0}}\mathrm{+}{\mathit{n}}_{\mathrm{1}}\mathit{R}\prime$(2)where n₀ and n₁ are derived empirically from observations. Using Eq. (2), we obtained the theoretical total number of sunspots expected on the solar surface at each time step. The conversion n_0,1 and dispersion d_0,1,2 coefficients used in our simulations are given in Table 1. By subtracting the number of sunspots already existing at this step, one can determine the number of new sunspots to be generated at this step in the simulation. We finally show in Fig. 1 (bottom panel) the power spectra of the observed R and simulated input R′ Wolf number for comparison purposes. The periodograms are quite different, especially at the rotational period and its first harmonic. This is expected as we do not introduce the rotational modulation (which is naturally present in the observed Wolf number) in our input Wolf number R′, but rather at the next step of the simulations (see above). For the high-amplitude (and short-period) variations of the observed Wolf number in the high-activity solar period, we consider that these variations are not due to a lack of spots at some moments in our model compared to the observations, but rather to the presence at these times of a group of large spots or an active cluster in the solar observations (see below for the spot size and spot filling factor distributions). Our model would require a more complex spot distribution function to reproduce such occasional appearances of large spot groups and we decided not to complicate it at this stage. We assumed an activity cycle length of 12.5 yr (length of solar cycle 23) and a time step of one day, leading to a total number of 4566 iterations or days. Furthermore, we added a random scattering of the time steps, with an amplitude of 4 h around the regular daily period, to mimic real observations.

• 2.

  Spatial and temporal distributions of the structures. At the beginning of the activity cycle, sunspots appear at medium latitudes, with some dispersion. As the activity level increases, the mean appearance latitude of sunspot slowly decreases in absolute value. This migration of active regions towards the equator continues during the whole cycle leading to a low mean appearance latitude at the end of the cycle. In our simulations, the mean latitudes of sunspot appearance at the beginning and at the end of the cycle were fixed at ±22° and ± 9°, respectively. These values are derived from the spot catalogues we used in Papers I and II. A linear law was then used to fit the decrease of the sunspot appearance latitude during the cycle, and we added a latitude scattering to reproduce the butterfly diagram. We also had the possibility of creating an asymmetry between the activity patterns in the northern and southern hemispheres, i.e. by adding more spots in one of them. As such a significant asymmetry has not been detected on the Sun during past solar cycles, we decided to put the same proportion of spots in each hemisphere. It is well known that spot groups do not to appear at random longitudes, but around the so-called active longitudes. Such preferred longitudes have been detected on the Sun and on several other stars (Berdyugina & Usoskin 2003; Ivanov 2007; Lanza et al. 2009, 2010). Their origin is still partly unexplained and is probably due to the dynamo process. Some trends seem to emerge from observations : there are generally two (seldom three) persistent active longitudes per hemisphere, shifted by 180°. They form a rigid structure, although they cannot be fixed in a reference frame because of the differential rotation (Berdyugina & Usoskin 2003). On the Sun, magnetic fields and new activity seem to emerge at locations where activity is already present (see e.g. Harvey & Zwaan 1993). Thus, instead of putting active longitudes with fixed longitude values in our model, we decided to build active longitudes from the position of already present sunspots. A fraction of new sunspots then appears in a restricted area in longitude around already present sunspots, while the remnant fraction is uniformly distributed in longitude. The fraction of new sunspots appearing in theses active areas and the longitude extension of the latter are reported in Table 1.

• 3.

  Large-scale dynamics. We describe here the global motion of magnetic features at the solar surface during their lifetime. This motion can be decomposed into two components, namely a longitudinal and a latitudinal motion. The Sun is well known to exhibit a differential rotation in latitude: active structures rotate faster at the equator than at the poles. The longitudinal motion or rotation rate of the active structures ω (in degrees per day) is therefore a function of the latitude θ, according to the following equation (e.g. Ward 1966; Meunier 2005): $\mathit{\omega }\mathrm{=}{\mathit{\omega }}_{\mathrm{0}}\mathrm{+}{\mathit{\omega }}_{\mathrm{1}}\mathit{.}{\sin }^{\mathrm{2}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{+}{\mathit{\omega }}_{\mathrm{2}}\mathit{.}{\sin }^{\mathrm{4}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{.}$(3)The latitudinal motion or meridional flow M of magnetic features is derived from Komm et al. (1993). It is a poleward motion in each hemisphere, which is also a function of the latitude : $\mathit{M}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathit{\alpha }\mathit{.}\sin \mathrm{\left(}\mathrm{2}\mathit{\theta }\mathrm{\right)}\mathrm{+}\mathit{\beta }\mathit{.}\sin \mathrm{\left(}\mathrm{4}\mathit{\theta }\mathrm{\right)}\mathit{.}$(4)

  Table 1

  Input parameters.

• 4.

  Spot properties. Spots appear preferentially with a given size, increase in area rapidly before slowly decreasing. The spot growing phase is very quick (on average 10 to 11 times faster than the decay phase, see Howard 1992). It is then often shorter than our one-day timestep for the small and mid-sized sunspots and of a few days for the largest ones. As we focus here on long-period planets, and as the active structure growing phase is smaller than or equal to our simulation time step, we decided not to include growing spots in our model for the moment and assumed that spots appear with their maximum size and then decrease during their lifetime. We will, however, include a description of the growing phase in future works to better describe the variations at small timescales. Thus, the spot distribution can be described with only the initial size distribution and a decay law. According to the La Laguna classification (see Martinez Pillet et al. 1993, and references therein), one can distinguish between isolated spots (La Laguna type 3) and spots belonging to complex groups (La Laguna type 2). In compliance with this classification, we built two distinct spot distributions. Both obey similar evolution laws, but with different parameters. The initial spot size distribution can be fitted with a log-normal law, according to Baumann & Solanki (2005). The two input parameters are therefore the mean spot size and the standard spot size deviation. We took the parameters used in our simulations from the snapshot model developed in Baumann & Solanki (2005). The spot decay law was parameterized by Martinez Pillet et al. (1993) as a log-normal law, with two input parameters: the mean and median decay values. We also put an upper size limit for the two types of spots, and a lower threshold of 10 μHem for all spots, chosen to be in good agreement with the observed spot distributions used in Paper II.

• 5.

  Facula properties. On the Sun, spots and faculae form active regions, spots being surrounded by large faculae. The ratio between the surface covered by faculae and spots (hereafter the size ratio) has been studied (Chapman et al. 1997, 2001, 2011). The authors found an average size ratio ranging between about 13 and 45 over solar cycles 22 and 23, depending on the method used to estimate the facular contrast. However, this corresponds to an average ratio measured for structures that can be at any state of evolution and any size. For our simulation, we need instead an initial size ratio. Each time we add a spot, we add a facula at the same place. The facula size distribution then depends only on the initial size ratio. As there is little in the literature about this quantity, it is one of our few free parameters. We assume a log-normal law to describe the initial facula-to-spot size ratio. Facula decay is then described by two processes:

  • First the “classical” decay is described as for dark spots, with alog-normal law. This decay corresponds to a certain surfacebeing lost by the facula at each time step of the simulation.

  • Second, a given proportion of this facula lost surface is converted into network (the rest being completely lost, for example by flux cancellation or submergence): facula fragments break away from active regions and diffuse over the surface.

  As the facula distribution is based on the spot distribution, the facula growing phase is not reproduced for the moment in our model. The facula growth and decay are longer than the spot values on average, owing to their longer lifetimes (Howard 1992), but the facula growth is still much shorter than its decay by a factor of ~3 (see e.g. Howard 1991). We decided not to reproduce this parameter to remain coherent with the spot distribution and as we focus mainly on timescales between a fraction of the rotational period and a complete activity cycle, and on long-period planets.

• 6.

  Network properties. We chose to make the network originate in the decay of faculae (see point 5). Network behaviour is then managed by its diffusion coefficient and by its decay rate. The network diffusion coefficient has been studied in many papers, which provide a wide range of values (Cadavid et al. 1999, and references therein).

  Fig. 2 Upper panel: observed (black line) and simulated (red line) spot sizes. Lower panel: the same for bright feature sizes.

  As these authors pointed out, there is a large gap between the diffusion coefficient values found by modelling (which are usually higher than 500 km² s^-1) and the values found with magnetogram analysis (which span from a few tens of km² s^-1 to less than 300 km² s^-1). They also raised the possibility that the diffusion coefficient is time-dependent.

  We finally decided to keep the value of 300 km² s^-1 used by Schrijver (2001) as it stands well between the two extreme groups of values. We assume the diffusion process to be isotropic (Cadavid et al. 1999). As for the decay rate, we simply described it by fixing the network remainder fraction at each time step (it is one of the model’s free parameters). We finally put a minimal size of 3 μHem for all bright structures (faculae or network), again in compliance with Paper II.

2.3. Comparison with the observed activity pattern

Here we compare the structures generated with our parameterized model to the observed structures used in Paper II. To do this, we use the spot and bright feature lists over cycle 23 that we retrieved from sunspot catalogues and MDI/SOHO magnetograms, respectively (Paper II). These data sets have some gaps, leading to a coverage of 3586 days, for a total duration of 4171 days, whereas our simulations have a length of 12.5 yr with a daily sampling, i.e. 4566 days (length of cycle 23). We therefore applied the same calendar as in Paper II to our outputs so as to make a relevant comparison. The full time series will be studied in Sect. 4.

2.3.1. Size distributions

The size distributions (in μHem) of our simulated structures are given in Fig. 2 and are compared with the observed values from Paper II. We note that very few of the input parameters detailed in Table 1 are free parameters. For most of them we retrieved the input laws and parameters from the literature so as to build a fully parameterized model of solar-like activity.

For the dark spot properties, we have no free parameters (i.e. they all originate from the literature). Remarkably, we find a good agreement between the spot size distribution simulated this way and the observed one. The lack of a few large spots (about 40 in all) for the simulated distribution is explained by the fact that the log-normal law usually used in the literature slightly underestimates the number of very large spots (see also Figs. 2 and 3 in Baumann & Solanki 2005). The short-term differences reported above in the Wolf number distributions are then likely to originate in the occasional appearance of these very large spots (and not to a global lack of spots in our model, as the simulated distribution nicely fits the observed one). For smaller spots, the observed distribution is discretized, which partly explained the gap (the remaining difference is due to slight differences between the observed distribution and the log-normal law). This small discrepancy in the distribution for small spots will not influence significantly the RV signal because the RV effect of small spots is negligible compared to the larger ones.

For the bright features (i.e. faculae and network), there are three free parameters (the facula-to-spot initial size ratio distribution, the facula fraction recovered in network, and the network decay rate). These parameters are poorly documented, so we decided to adjust them to better fit the observed distribution. Overall, the simulated size distribution is in good agreement with the observed distribution (Fig. 2). However, there are some remaining differences that are inherent to the model parameters. The size distribution for all observed bright features does not follow a log-normal law. For the simulated bright features, we note that we chose to make the facula initial size distribution directly dependent on the spot distribution, and the network originate in the decrease of the faculae. This explains thus the “flattened-S” shape of the simulated distribution, which corresponds to the addition of the facula and network respective distributions. Reproducing even more accurately the shape of the observed distribution would require making our model much more complex, for example by making the bright feature decay size-dependent or cycle-dependent. Given our goals, such a level of complexity is not necessary. As larger faculae have the largest influence on the RV, we decided to fit the distribution in priority for the highest sizes (i.e. for sizes greater than 5000 μHem) by adjusting our free parameters. There is still a discrepancy between the observed and simulated distributions for the smallest sizes, but we consider that its effect on our observables will be mostly negligible at our level of precision.

2.3.2. Latitude distributions

Another way to compare our simulated pattern to the observed one is to study the active structure latitude distributions. The latter are given in Fig. 3. By doing so we will be able to validate the large-scale behaviour of our model. The latitude distributions of observed and simulated structures are indeed in good agreement. The kind of oscillations that we can distinguish in the observed bright feature latitude distribution at high latitudes and with a one-year period comes from a seasonal effect.

The solar observations are not made exactly in the plane of the solar equator, and the active structure distribution seems to be slightly shifted towards one visible solar hemisphere or the other, depending on the observation time, owing to the noise level in the observed magnetograms². This effect is also visible on the bright feature latitude histogram, where the observed distribution is spread a little more towards higher latitudes.

2.3.3. Filling factors

Moreover, we compare the filling factors of the projected active structures over the full solar cycle for the observed and simulated patterns in Fig. 4. We find the observed and simulated filling factor time series to be in very good agreement. The main difference between the observed and simulated distributions is found during the high-activity period: the activity peak present in the observations between approximately JD 2 452 200 and JD 2 452 400 is not accurately reproduced in the simulations. This difference does not come from a global lack of active structures in the high-activity period of our model because i) the simulated daily spot number comes from Wolf number observations and ii) the facula distribution directly depends on the spot distribution (Sect. 2.2). The difference instead originates in an occasional concentration of a few great structures in the observed pattern that is not reproduced in our model. Thus, this discrepancy is not a bias of our model but can rather be considered as statistical noise. In the case of the bright features only, we also note a slight discrepancy between the observed and the simulated filling factors during the low-activity period (the simulated filling factor being higher than the observed value). We attribute this to the larger number of very small bright features injected in our model (see above).

Finally, we compare the ratio of the bright feature and spot filling factors over the solar cycle in Fig. 5 to estimate the evolution of the facula-to-spot size ratio. We find the ratios for the observed and simulated time series to be in very good agreement. One can also note that the size ratio reaches high values (from 50 to 150) during the low-activity period (i.e. at the beginning and at the end of the solar cycle), whereas it is much lower (around 20) during most of the cycle because the number of dark spots is very low (equal or close to zero) during periods of low activity, whereas there is always some “background noise” due to bright features, and especially the network.

We find the average value of the size ratios over the total cycle to be very close to each other, with values of 25.7 and 26.0 for the observed and simulated structures, respectively. These values are in agreement with the studies of Chapman et al. (1997, 2001, 2011). We also note that the size ratio time series displayed by Chapman et al. (2011) look very much like ours over cycle 23, with peaks above 100 during low-activity periods and minima around 25 during high-activity periods.

Overall, we can conclude that our simulated solar activity pattern closely matches the observed one and is thus a reliable model of a solar-like star magnetic activity. However, to definitely establish its validity, we compare the resulting RV and photometric time series. In particular, we check the impact of the differences between the two patterns on the time series.

3. Comparison with the Sun observations

In this section, we first describe the way we simulate the time series for the different observables and how we take into account the stellar properties (limb-darkening, active structures temperature contrasts, impact on convective blueshift, physical and geometrical stellar properties). We then compare the RV time series simulated in the case of an edge-on solar-like star to the RV time series obtained with the observed activity pattern described in Paper II to establish the validity of our model.

3.1. Description of the simulations

3.1.1. Limb-darkening

In our previous papers, the center-to-limb-darkening ratio was linear with respect to μ (with μ = cos(θ), θ being the angle to the centre of the solar disc, Desort et al. 2007), $\mathit{I}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{+}\mathit{ϵ}\mathit{\mu }$(5)with ϵ = 0.6. In this case, the limb-darkening was not temperature-dependent and therefore applied to both the inactive photosphere and the active structures indiscriminately. Here, we decided to change our limb-darkening law, so as i) to have a more accurate law; and ii) to use a temperature-dependent law, that could be adaptable to different types of stars. We took the non-linear limb-darkening law from Claret & Hauschildt (2003), $\frac{\mathit{I}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}}{\mathit{I}\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{=}\mathrm{1}\mathrm{-}\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{\mathrm{4}}{\mathit{a}}_{\mathit{k}}\left(\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\frac{\mathit{k}}{\mathrm{2}}}\right)$(6)where I(1) stands for the intensity at the centre of the stellar disc. The four Claret limb-darkening coefficients a_k are the bolometric coefficients taken from ATLAS models (see Claret & Hauschildt 2003). These coefficients are temperature-dependent and the effective temperature is one of the input parameters, which is why such a law is better adapted to our simulation tool, and also why and such a law can be extrapolated not only for the Sun, but for any star for which the effective temperature is known. This also means that the limb-darkening coefficients applied to the spotted stellar surface will be slightly different from the ones applied to the inactive photosphere. We therefore use a different limb-darkening value for the photosphere and for the dark spots, depending on their respective temperatures. As faculae show a very strong contrast variation due to more complex processes depending on their position on the disc, we directly define their contrast C_pl(μ) with respect to the limb-darkened photosphere (see below) without using the Claret law according to their temperature.

3.1.2. Spot and bright feature contrast and photometric time series

To estimate the contrast of an active structure (dark spot or bright feature) with respect to the quiet photosphere, we use the procedure described in Paper II. This procedure is an independent, preliminary step to our simulations (i.e. to the generation of the spectra and of the RV and photometric time series), and is aimed at determining the active structure contrasts that we will use as input parameters. Starting from a given range of values for the spot and bright features contrasts, and based on observed structure patterns, we build (for each set of contrast values) two time series representing the respective relative contributions of dark and bright active structures to the stellar irradiance. A χ² minimization between the sum of these contributions (quiet photosphere, dark and bright features) and the observed total solar irradiance (TSI) taken from Fröhlich & Lean (1998) is then performed over solar cycle 23. We adopt the same value for the quiet Sun reference as in Paper II, i.e. 1365.46 W m^-2, which is very close to the average value found by Crouch et al. (2008) with their TSI model over twelve solar cycles. Contrary to Paper II, we include in the procedure the influence of the centre-to-limb-darkening, so as to take into account the Claret limb-darkening law we will now use in our simulations. We end up with a spot temperature deficit ΔT_sp = −605 K and a facula contrast³C_pl = 0.131618−0.218744μ + 0.104757μ². We then use these contrasts as our input parameters in the simulations. As in Paper II, we finally compare the sum of the photometric contributions of spots and bright features obtained with our simulations with the observed TSI of Fröhlich & Lean (1998), so as to check the validity of our activity model. The comparison is done on 2263 points, ranging from 1996 to 2003.

The irradiance obtained with our simulations and the observed TSI of Fröhlich & Lean (1998) are shown in Fig. 6. Our values match the observed irradiance quite closely, except for some peaks in the high-activity period that are not reproduced in the model because our model does not reproduce well the occasional appearance of very large active structures or large active structure clusters (see Sects. 2.3.1 and 2.3.3). Consequently, the distribution of the simulated TSI is slightly narrower than the observed TSI value, and the averaged rms of the simulated TSI is about 15% lower in comparison (see Table 3). However, the temporal evolution of the simulated TSI rms closely matches the observed evolution. As already stated, we consider that these differences do not affect the significance of our simulations as a model.

3.1.3. Attenuation of the convective blueshift

For all the following simulations we adopt the same value for the attenuation of the convective blueshift as in Paper II, i.e. 190 m s^-1. We refer to this paper for justification.

3.1.4. Building the spectra

We build the spectra from the input structure lists as described in Desort et al. (2007) and in Papers I and II. We assume a stellar mass of 1 M_⊙, a temperature T_eff of 5800 K and a rotational velocity of 1.9 km s^-1 at the equator. In this section we keep the stellar inclination i to 90°, corresponding to a star seen edge-on. Briefly, we use a synthetic spectrum from the Kurucz models (Kurucz 1993) corresponding to a G2V star and apply it to each cell of the visible stellar 3D hemisphere divided into a grid. The spectrum is shifted to the cell RV. In the case of the presence of an active structure, it is weighted with a black-body law, taking into account the active structure contrast with respect to the quiet photosphere. To include the effect of the attenuation of the convective blueshift, the cell spectrum is redshifted by 190 m s^-1. This attenuation has other effects on the spectral lines; however, Dumusque et al. (2014) show that only considering the RV shift is sufficient in order to estimate the RV effect of active regions.

All cell spectra are then balanced by the cell’s projected surface and limb-darkening. We finally sum up all cell contributions to obtain the stellar spectrum.

3.1.5. Computation of the RV time series

As in Paper II, we compute the RV using our Software for the Analysis of the Fourier Interspectrum Radial velocities (SAFIR, see Galland et al. 2005) on the built spectra as if they were actual observed spectra. We use only the wavelength range corresponding to the order ♯ 31 of the HARPS spectrograph, as done in Papers I and II.

We obtain three distinct time series: two due to the respective photometric contributions of dark spots and bright features (hereafter the spot and facula time series, respectively); and the third due to the partial inhibition of the convective blueshift in the active structures (hereafter the convection time series). For the last of these, we consider only the inhibition of the convective blueshift (i.e. this effect is not weighted by the active region flux). We sum the three time series to obtain the total RV variations. Here we consider this sum to be a good approximation of the real RV variations given that the convective blueshift is dominant (see below).

3.2. Validation of the activity model

3.2.1. RV time series based on the observed patterns

In this section, we first check if the change of limb-darkening law and the corresponding change of the active structure contrasts have a significant impact on the resulting RV time series. To do so, we compute new RV time series corresponding to the observed solar activity patterns used in Paper II, but this time with the new limb-darkening law and structure contrasts. We then compare these new RV time series to the ones obtained in Paper II, using the same temporal sampling. We find the RV time series to be closely correlated and thus in good agreement. We display in Fig. 7 the RV time series computed with the Claret limb-darkening versus the RV time series computed with the linear limb-darkening. We find the slope of the fits and the correlation of the RV time series to be close to 1 in both of the cases: the spot+facula RV time series and the total RV series. Despite the simpler limb-darkening law and a different active structure contrast, the results shown in Paper II are still valid for estimating the RV effect.

3.3. Comparison between the RV time series based on the observed and simulated patterns

We compare here the RV time series based on the observed and simulated activity patterns (both computed with the Claret limb-darkening) so as to assess the validity of our activity model. We compare the RV time series for the different contributions (spot, facula, convection, and total). All the time series and their dispersion over the cycle and their histograms, are displayed in Fig. 8. We find the amplitudes of the observed and simulated RV time series to be in very good agreement for the spot and facula time series, with very similar histograms. The convection RV time series is closely related to the facula filling factor (with a Pearson correlation coefficient of 0.97) and widely dominates the total RV signal, as in Paper II. The main visible difference comes from the activity peak at around JD 2 452 200−2 452 400 in the RV reconstructed from observations, which is not echoed in the simulated RV. We have already discussed the origin of this difference in Sect. 2.3. In the low-activity period, the averaged amplitude of the simulated convection time series is about 20% higher than for the observed one (see Table 3). We attribute this small discrepancy to the excess of very small bright features in the model that we discuss Sect 2.3.

To have a better understanding of the RV signature of our simulated activity pattern, we also display the RV rms computed over 30-day intervals over the cycle in Fig. 8. This gives an idea of the temporal evolution of the RV time series dispersion. We find the evolution of the RV dispersion during the solar cycle to be in good agreement for the observed and simulated time series. We also provide the RV rms for the different components taken over the complete cycle and for low- and high-activity periods in Table 2, and the RV amplitudes in Table 3. The low- and high-activity periods are the same as in Paper II. We note that:

• for the spot and facula simulated time series, the rms is slightlyhigher than for the observed time series in the low-activity period(about two times higher in the case of brightfeatures). However, this is not the case for theconvective time series;

• in the case of the convective component, the rms over the complete cycle is about 24% lower for the simulated data than for the observed data. When looking at the temporal evolution of the RV rms in Fig. 8, we can see that this discrepancy comes mostly from a few high-activity peaks in the observed time series.

We finally compare the Lomb-Scargle periodograms of the RV time series based on the observed and simulated activity patterns. The periodograms are shown in Fig. 9 in the cases of the “photometric” contributions (spots and faculae) and of the sum of all contributions. The respective periodograms for the observed and simulated patterns remarkably present the same characteristics:

• For the “photometric” component, the power is mostlyconcentrated at the solar rotation period and its two firstharmonics.

• For the sum of all contributions, there is still power at the rotation period, but there is much more power at longer periods (between 300 and 1000 days) due to the long-term activity cycle.

Nonetheless, we note a significant difference: for the spot+facula RV signal, the peaks at half the rotation period are emphasized compared to the peaks at the rotation period in the case of the observed pattern, whereas it is the contrary for our simulated pattern. We already noticed and discussed in Paper II the predominance of the power at half the rotation period in the RV time series derived from observed solar patterns. We explained it by the presence of two symmetrically active longitudes on the solar surface. Since here we do not impose two active longitudes separated by 180° over the whole time series but rather variable active longitudes, this effect may not be very important. An alternative explanation could be that the difference between the periodograms corresponding to observed and simulated patterns may originate in the active structure quick growing phase (which is at present not reproduced in our model). This may also explain the smaller RV and photometric rms when compared to the observations.

As detailed above, small discrepancies in amplitude and dispersion on short timescales are present both in the photometry and in the RV between the time series derived from observed and simulated activity patterns. We consider that these differences likely originate in two sources: first, the spot size distribution, where the few larger spots are hard to model; and second, the active structure growing phase. Overall, we consider that the temporal evolution over the complete cycle of both the amplitude and dispersion of the simulated time series closely match the time series derived from the observed activity pattern. Hence, we conclude that: i) these short-term discrepancies are not significant enough to question our model reliability and ii) the comparison between the RV time series derived from observed and simulated patterns then assess the overall validity of our model.

4. Photometric and RV time series of inclined solar-type stars

For solar-type stars of the same age as the Sun, it is commonly assumed that active structures are mostly concentrated in a belt around the stellar equator (even if bright features are much more dispersed in latitude). A different inclination of the stellar rotation axis should then have a significant effect on the various time series that need to be investigated. In this section, we perform the same simulations as above, but for different inclinations of the stellar rotation axis. We consider inclinations between i = 10° (i.e. for a star seen nearly pole-on) and i = 90° (star seen edge-on), with a sampling of 10°. The time series are studied with their original temporal sampling, i.e. 4566 days and no gaps.

4.1. Photometric time series

The impact of the inclination i of the stellar rotation axis on the long-term solar irradiance variations was studied by Knaack et al. (2001). The authors computed the solar irradiance corresponding to a three-component model (quiet Sun, dark spot, and bright facula) at activity extrema and for a variable inclination. For their active region distributions, they used simple active latitude belts with contrast corresponding to the given active structure (i.e. no individual structure was introduced in their model). They found that the apparent active structure surface coverage decreases with a decreasing i. They nonetheless expected an increase of the TSI with a decreasing i since bright features are limb-brightened while the contrast of dark spots is roughly independent of the limb-darkening. Thus, when going from an edge-on to a pole-on configuration, the impact of dark spots would decrease following their apparent covered surface. On the contrary, the decrease of the apparent surface covered by bright features would be at least compensated by their increased contrast since we would see them mainly on the limb. In Fig. 10, we show the evolution of the spot and facula projected filling factors with i. We find the projected spot filling factor to decrease by about 35% when going from edge-on (i = 90°) to nearly pole-on (i = 10°), and the projected facula filling factor to decrease by about 43%. These results are quite similar to those found by Knaack et al. (2001). We note that the filling factor does not tend towards 0 for i ≃ 0°.

The relative variation of the TSI during the cycle (which is approximately 0.1% of the quiet Sun irradiance) versus i is displayed in Fig. 11 (upper panel). We find that it increases by only 14% when going from i = 90° to i = 10°. This is much smaller than the 40 ± 10% increase predicted by Knaack et al. (2001) with comparable input parameters but a much simpler model. However, it points towards the same verdict proposed by Knaack et al. (2001), i.e. that photometric variations of seemingly inactive Sun-like stars cannot be explained by an inclination effect. We conclude that for a solar-like star with a similar activity configuration, the effect of an inclined rotation axis on the long-term variations of the total irradiance (approximately the solar cycle length) is relatively small. On the contrary, the TSI short-term variations (which were indeed not studied by Knaack et al. 2001) are strongly affected by the inclination. As illustrated in Fig. 11 (lower panel), we find the TSI short-term dispersion to be decreased by a factor of ~6 when going from i = 90° to i = 10°. A possible explanation for such a decrease is the following:

• First, for smaller inclinations, we mainly see the effect of brightfeatures because they are more extended towards the higherlatitudes than the dark spots.

• Then, due to the activity configuration (where active structures are mainly located on two belts on both sides of the solar equator), for small inclinations we see the same structures during all of the rotation period and not during half a period.

Therefore, in the case of a star seen nearly pole-on, the short-term dispersion of the irradiance (approximately a few rotation periods) does not originate when the structure crosses the visible hemisphere for each half rotation period. It originates only when the structure appears and decays. To confirm this, we display in Fig. 12 the Lomb-Scargle periodograms of the simulated TSI for representative inclinations. For a star seen edge-on, the periodogram is dominated by power at the rotation period of the star. When going towards smaller inclinations, the signal at the rotation period gradually decreases until it disappears completely for a star seen nearly pole-on. On the contrary, the signal at a much longer period (which is likely induced by cycle-related long-term periodicities equal to the cycle length) becomes increasingly important for smaller inclinations. A reason for which the signals at long-term periods (in the 600- to 1500-day range) become increasingly dominant with decreasing i in the TSI (and to a lesser degree in the total RV) periodograms could be the following: for nearly pole-on configurations, we see the active structures on one stellar hemisphere only, whereas we see them on two hemispheres for configurations closer to edge-on. This may induce a kind of averaging effect on the long-term signal for the edge-on configuration and explain its increase with decreasing i. Finally, we show in Fig. 13 the simulated TSI for representative configurations. As we found above, the long-term amplitudes of the two time series are nearly the same; on the contrary the dispersion is widely reduced with i.

4.2. RV time series

We now study the impact of stellar inclination i on our simulated RV time series. Our first result is that in contrast with the photometry, both the amplitude and the dispersion of the total RV decrease with decreasing i. This is nicely illustrated in Fig. 14 where we display the RV time series for all activity components and for significant inclinations. The amplitude and the dispersion (RV rms) of the total RV taken over the complete cycle are decreased by a factor of ~6 when going from an edge-on to a nearly pole-on configuration.

We characterize both the long-term and short-term variations of the RV signals to investigate thoroughly the impact of inclination. Our results are illustrated in Fig. 15. The “long-term” variations correspond to the signal taken over the complete stellar cycle. To study the short-term variations of the signal, we perform a running average of the RV time series with a smoothing window of 30 days and subtract it from the original data. Then the long-term variations should primarily be affected by the global activity cycle, whereas the short-term variations will come from rotation-related effects.

1. Peak-to-peak amplitude: For the total RV andwhen going from i = 90° to i = 10°, the amplitude decreases by ~80% over the cycle and by ~85% on the short term. When taking only the photometric contribution (i.e. spots and faculae) into account, the peak-to-peak amplitude is decreased by nearly 90% on the long-term, and by 87% on the short-term.

2. Dispersion: the total RV rms decreases by 70% with inclination on the long term, and by 85% on the short term. In the case of spots and faculae only, the decrease is the same on both timescales and is of nearly 85%. For a star with a solar-like activity pattern seen almost pole-on, a jitter of ~0.5 m s^-1 can be expected if the observation timescale is equal to the activity cycle length, and a jitter of ~0.2 m s^-1 can be expected over a month.

3. Ratio between photometric and convective components: we also study the evolution of the relative contribution of the photometric (i.e., due to spots and faculae) component to the total RV signal with the inclination. Over the cycle, the photometric fraction of the RV is reduced by a factor of ~2 when going from i = 90° to i = 10° (i.e. the convective component is increasingly important with a decreasing inclination). On the contrary, on the short-term the photometric relative contribution remains nearly constant with i, at a level of ~0.3. This value is at least two times larger than the value of the photometric fraction on the long-term (i.e. ~0.16 in the edge-on configuration). This confirms that the spot+facula component has mostly a short-term effect on the RV signal, equal to the rotation period.

We finally study the Lomb-Scargle periodograms of the RV time series for representative inclinations, first in the case of the sum of the flux contributions of dark and bright structures, and then when adding the convective component. The periodograms are shown in Fig. 12. Starting from i = 90°, the RV periodogram for the spot+facula signal is dominated by the stellar rotation. The power is concentrated at the rotation period and its two first harmonics, with about two times and six times less power for the first and second harmonics, respectively. The power is induced by the active structure crossing of the visible hemisphere during about half the rotation period (as the active structures are mainly located near and around the stellar equator). When i decreases, we begin to see the active structures on a time longer than half the rotation period, until we see them permanently for the nearly pole-on configuration. There is thus nearly no remaining component in the spot+facula RV periodogram for i = 10°. This is in good agreement with the 90% decrease in jitter already described.

Starting from the edge-on configuration, the total RV periodogram is dominated i) by a series of peaks at long period (200−1000 days) due to the long-term activity cycle and ii) by peaks at the stellar rotation period and its three first harmonics.

As is expected, the peaks corresponding to the rotation period and its harmonics gradually decrease when going towards smaller inclinations, until they nearly disappear for i = 10°. We also note that the harmonics of the rotation signal decrease significantly faster than the rotation signal itself. On the contrary, the long-term signal corresponding to the activity cycle decreases slowly when i decreases, but remains far above the 1% false-alarm probabilities (FAP) for i = 10°. In the edge-on configuration, the power is equally distributed between the rotation and the long-term signals, but the latter is greater for smaller inclinations. This is in agreement with the decrease of the ratio of the RV photometric component over the total RV signal on the long-term we show above.

We conclude that in the case of a star with a solar-type activity pattern, the stellar inclination strongly affects both the observed stellar irradiance, RV variations, and active regions filling factor (and hence the observed chromospheric activity). Indeed, the short-term variations of these various observables strongly decrease with a decreasing inclination. Then, on short to medium timescales (from a few days to a year, or well under the activity cycle length), it will probably not be possible to distinguish between a solar-like star seen nearly pole-on and an inactive star. Provided that the time baseline is sufficient, on timescales equal to the cycle length, the activity cycle should still be detectable for small inclinations, but with a reduced amplitude. It would then not be possible to distinguish it from a nearly inactive star with a low-amplitude activity cycle.

4.3. Potential impact of stellar inclination on various correlations

4.3.1. Correlation between total RV signal and facula filling factor

In Paper IV, we used the solar calcium (Ca) index (S-index or log$\mathrm{\left(}{\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}\mathrm{\right)}$) as a tool to partially correct the activity-induced RV from its convective component (using the observed solar activity pattern). The Ca index is correlated with the chromospheric plage filling factor and hence with the photospheric facula filling factor (Shapiro et al. 2014, found a linear dependence between the Ca index and the facula filling factor in the case of the Sun).

As already explained, in the case of a solar-like star seen edge-on, the RV signal is closely correlated to the facula filling factor (as the convective component dominates the RV variations, see Sect. 3.3 and Figs. 4 and 8). When going towards smaller inclinations, we note that the correlation gets even stronger, with a Pearson correlation coefficient going from 0.97 for i = 90° to 0.99 for i = 10°. This is in agreement with the fact that the relative contribution of the spot+facula component to the RV signal decreases with i. This means that the Ca index correction method can be used regardless of the stellar inclination.

4.3.2. RV jitter, Ca index and facula filling factor

In recent studies (Wright 2005; Santos et al. 2010; Isaacson & Fischer 2010; Hillenbrand et al. 2015), the average level of RV jitter was commonly taken as a proxy for the stellar magnetic activity level, along with the Ca index. For different samples of FGK stars, the authors found loose correlations between the average RV jitter and the mean Ca index, generally with a significant amount of dispersion. In our simulations, we find that the ratio of the mean facula filling factor over the RV rms (calculated over the complete activity cycle) is not constant over the range of inclinations i we explored. We display in Fig. 16 the facula filling factor to RV amplitude and to RV rms ratios. We find that both ratios increase towards smaller inclinations (in agreement with Sects. 4.1 and 4.2, where the decrease of the RV amplitude and rms with i is more pronounced than the decrease of the facula filling factor with i). Provided that the Ca index evolves in the same way as the facula filling factor with i, this could partly explain the large dispersion in the RV rms-mean Ca index relation found in the studies cited above (as the stellar inclination remains unknown for most of the observed targets). It also means that a clear RV rms-Ca relation will remain hard to establish for future studies, unless the uncertainty on the stellar inclination can be removed.

4.4. Detection limits

In this section, we compute the detection limits for the total RV time series over the range of stellar inclinations explored above.

4.4.1. Approach

We use two different methods: the correlation-based method and the local power amplitude (LPA) method. Both were described in detail and tested on real stellar RV data in Meunier et al. (2012). Briefly, the first method makes the correlation between the periodograms of a generated planetary RV signal (with the same temporal sampling as the data) and of the actual RV data to which the planetary RV signal has been added (i.e. it determines the correlation of the power of the (stellar data + fake planet) periodogram versus the power of the fake planet periodogram). The detection limit corresponds then to the minimal mass and period for which the correlation values are all above a given threshold, for 100 realizations spanning all orbital phases. The threshold (spanning here from 0.003 for i = 10° to 0.05 for i = 90°) corresponds to the maximum of the correlations obtained for a very low-mass planet (here ~0.6 M_Earth). For the LPA method, we compare the periodograms of the actual RV data and of a given generated planetary signal (for the same temporal sampling), but within a localized period range around the given planetary period. The detection limit at this period corresponds then to the minimal mass for which the maximum power of the planetary RV periodogram (within the limited period range) is always above the maximum power of the actual RV data periodogram in the period range.

We compute here our detection limits for only one period of 480.1 days, corresponding to a separation of about 1.2 au, as in Paper II. It is a representative value of the outer boundary of the HZ for a solar-type star. The total time span is always fixed to the complete simulation time span (i.e. 4566 days). As in Papers I and II, we assume the star to be observed eight months a year (meaning that we take into account only 2978 RV points instead of 4566). We assume different temporal samplings: all points (i.e. one-day sampling or 1:1); 1 point every 4 nights (1:4); 1 point every 8 nights (1:8); or 1 point every 20 nights (1:20), which means that the (1:20) series has 20 times fewer points that the (1:1). We also assume different precisions on the RV (no added noise, 0.01, 0.05, and 0.1m s^-1 noise levels), as we did in Paper I.

Finally, we compute the detection limits versus inclination for two different cases. In the first case, we compute the detection limits considering that the hypothetical planet is always seen orbiting edge-on for all stellar inclinations (i.e. there is an increasing spin-orbit misalignment with decreasing inclination). In the second case, we consider that the hypothetic planet always orbits in the stellar equatorial plane (i.e. there is always spin-orbit alignment).

Considering these two distinct possibilities is important; in most cases when searching for planets with RV, we do not know the stellar inclination or the inclination of the planet orbit with respect to the line of sight. The inclination of the planet orbit is taken into account by giving minimal masses m.sin(i) (with i denoting here the inclination of the planet orbit with respect to the line of sight) for detected planets. However, the detection limits are generally computed without knowing the stellar inclination itself, i.e. considering that the star is seen edge-on and that the activity-induced jitter is not reduced by the projection effect. The detection limits we compute in case 1 correspond to this configuration and thus to the best possible detection limits. As we go towards smaller stellar inclinations, they also correspond to more and more unlikely orbital configurations. Most of the exoplanets found so far orbit in or near the stellar equatorial plane⁴ (even if in the specific case of hot Jupiters, a significant fraction of systems show spin-orbit misalignments, e.g. Albrecht et al. 2012; Triaud et al. 2014).

On the contrary, the detection limits that we compute in case 2, i.e. when considering a spin-orbit alignment, correspond to a more conservative but more realistic assumption. The detection limits for the correlation-based and LPA methods are displayed in Fig. 17.

4.4.2. Results

Comparison with previous results –

For the edge-on configuration, we can compare the detection limits (for the 1:1 sampling) to the detection limits computed in Paper IV⁵, i.e. on RV time series derived from observed solar activity patterns, as the latter time series had a similar time span, sampling and number of data points. In the case of the correlation method, they are in fairly good agreement (5.5 M_Earth for the present time series versus 6.8 M_Earth for the time series from observed patterns), whereas there is a certain discrepancy for the LPA method (9.9 M_Earth and 15.7 M_Earth for the simulated and observed RV time series, respectively).

Comparison between the two methods –

The correlation-based method gives lower detection limits than the LPA method. This is also in agreement with our results from Paper IV. We note that we also compared the two methods in Meunier et al. (2012), this time on actual stellar data. Out of a sample of ten targets, we found in most cases that the LPA method gave lower detection limits. However in these cases, the targets were massive A-F main-sequence stars with medium to high vsini (from 7 to ≥175 km s^-1), some of them being young stars. They would then have a very different activity pattern than the solar one.

Impact of stellar inclination –

For the LPA method, the decrease in the detection limit with the stellar inclination i in case 1 (best possible detection limits) is consistent with the decrease in the activity-induced RV rms when going towards a smaller i. In case 2, we observe that the detection limits remain around 10 M_Earth for i ≥ 40° before increasing towards higher masses for smaller inclinations. This should mean that for i ≥ 40°, the decrease in the planetary-induced RV amplitude for a decreasing inclination of the system is counterbalanced by the decrease in the same time of the activity-induced RV jitter. For smaller inclinations, the decrease of the activity RV jitter is less pronounced and hence it becomes more difficult to detect the planet. This is consistent with the fact that the behaviour of the activity-induced RV jitter with respect to i is not sinusoidal.

For the correlation method, the behaviour of the detection limits is quite different. In case 1, the detection limits decrease with i until they reach a plateau at about 1.5 M_Earth for i ~ 40°. In case 2, surprisingly, the detection limits begin to decrease with a decreasing i and reach a minimum and best value of 2−2.5 M_Earth for i = 40−50°. For smaller inclinations, the detection limits increase with a decreasing i. Thus, the optimal configuration for the detection of the planet is an orbital plane inclined of 45° with respect to the line of sight (in the most probable case of spin-orbit alignment).

4.4.3. Impact of the parameters

In this section, we study in greater detail the impact of the main parameters for the computation of the detection limits (i.e. the added RV noise level and the temporal sampling. To this end, we first study their influence on the detection limits computed above (Fig. 17). Then, we extend the study to a wider range of noise level or temporal sampling, but this time concentrating on one case (i = 50°, case 2). This roughly corresponds to the inclination for which we get the better detection limits (in case 2).

Temporal window –

Our detection limit computation is based on the periodograms of the RV time series taken eight months a year so as to better mimic actual observations. The periodograms will thus be modified compared to the periodograms of the total RV time series taken over the full simulation time span as a result of aliasing. We compare the periodograms of the total RV time series (in the i = 90° case) before and after removing four months a year in Fig. 18. Removing four months a year injects a large amount of power in the RV periodogram around one year (~360−400 days) and at lower harmonics (180 and 90 days). It will then affect our detection limit computation and deteriorate our detection limits, but it has to be taken into account to reproduce the observations more accurately.

Temporal sampling –

We found in Paper IV that when taking into account all RV components, the detection limits got worse for the (1:20) sampling (i.e. for a largely degraded sampling), while they did not vary significantly when going from a (1:1) to (1:4) or (1:8) sampling. Figure 17 shows that the LPA detection limits do not vary significantly when going from the (1:1) to the (1:20) sampling. The correlation detection limits do not vary significantly when going from the (1:1) to the (1:8) sampling, but get worse for the (1:20) case.

We now focus on the i = 50° case to study the sampling impact in more details. We compute the detection limits for an extended sampling range (from (1:1) to (1:70)) and show the results in Fig. 19. We find that both the correlation and the LPA detection limits are nearly independent from the temporal sampling up to the (1:10) sampling, and get worse for more degraded samplings (even if they show a large dispersion in the LPA case). The large dispersion at large temporal samplings probably reflects the larger uncertainty on the detection limits for a smaller number of points.

Added noise –

The detection limits computed with the correlation and LPA methods are nearly independent from the added RV noise for a noise level up to 10 cm s^-1 (Fig. 17), i.e. for the best RV accuracy expected on future spectrographs. This is in agreement with the results of Paper IV. We consider that the small variations seen between the detection limits computed for the different noise levels in the case of the correlation method are not significant.

We study the impact of the added noise for an extended range of noise levels, up to 1 m s^-1 (i.e. the current RV accuracy reached on the better spectrographs such as HARPS), in the i = 50° case (Fig. 20). The LPA detection limits are constant up to 1 m s^-1, while the correlation detection limits remain constant up to 0.7 m s^-1 with a slightly larger value at 1 m s^-1.

5. Conclusion and perspectives

We built a fully parametrized model of the activity pattern of a solar-like star, including the dark spots and the bright features (large faculae and network). The model includes about 30 parameters that account for the different activity scales and the active structure behaviour (most of them being well constrained by the literature), and has a daily timescale over a complete activity cycle. Using the same approach as in the previous papers, we deduced the corresponding RV and photometric time series, taking the inhibition of the convective blueshift into account in the case of the RV. The simulated activity pattern, as well as the time series, were compared to the work done in Papers I, II and III with data from solar observations. We found our model to be in remarkably good agreement with the previous data, thus confirming its validity. We then study the case of a solar-like star seen under different configurations so as to estimate the impact of the stellar inclination on the activity-induced jitter. Our results are the following:

• For the stellar irradiance, the stellar inclination i has almost noeffect on the amplitude of the long-term irradiance variations.The decrease of the active structure filling factor with adecreasing i mostly counterbalances the enhanced contrast of thestructures, in agreement with some previous studies (Knaacket al. 2001; Shapiroet al. 2014). On the contrary, theshort-term jitter is reduced by a factor of 6 when going from anedge-on to a nearly pole-on configuration.

• For the RV variations, the peak-to-peak amplitude and the rms are reduced by a factor of 8−10 on both the long-term and the short-term. For a nearly pole-on configuration, the remaining RV jitter is about 0.6 m s^-1 when considering the whole cycle, and is about 0.2 m s^-1 on shorter timescales.

• When computed over the whole activity cycle, the activity-induced periodograms show power mainly in two period ranges: first at the stellar rotation period and its two or three first harmonics (induced by the crossing of active structures on the visible stellar hemisphere), and then at much longer periods, induced by more complex cycle effects. When going towards smaller inclinations, the power at the rotation period is reduced and nearly disappears for a pole-on configuration. The long-term cycle effects increase.

• In the case of a solar-like star, the convective component widely dominates the RV variations. The photometric contribution of spots and bright features account for about 35% of the RV jitter and decreases to about 15% for a nearly pole-on configuration.

• For a solar-like star seen in a pole-on configuration, the photometric and RV (and presumably the chromospheric) activity-induced variations are most probably not distinguishable from those of a less active or almost inactive star.

• Finally, when considering a realistic orbital configuration (i.e., spin-orbit alignment), the optimal configuration for planet detection is a system inclined by about 45°. In this case, the lowest detection limits reach planetary masses of about 2 M_Earth at 480 days, without applying a correction to the RV signal.

According to previous studies, solar-like stars show a great diversity of activity levels and properties (Schröder et al. 2013). The Sun itself is considered to be an average star, not particularly active but not quiet either. Being fully parametrized and validated for the solar case, our activity model now allows us to explore a wide range of activity parameters and stellar properties. In this paper, we have focused on solar-like activity and convection level stars; we expect stars with a lower activity level and/or stars with a lower convection level (such as K-type stars) to be much less affected and then to exhibit lower detection limits in the corresponding HZ. We also emphasize that activity correction methods (such as the calcium index or the H_α line) are commonly applied to RV time series and already allow detection limits to be reached at the M_Earth level in the HZ (Paper IV). Optimized observation and reduction strategies are also promising (such as averaging, see e.g. Dumusque et al. 2011a). Although it is not the focus of our paper, improving such methods and strategies, as well as looking for new ones, is an extremely important question. Future instruments with very high RV accuracy may therefore be critical in implementing very efficient correction tools to extract low-mass planet signals. Our goal is to be able to test and even reproduce the photometric and RV variations corresponding to each activity and stellar configuration. This should allow a better understanding of the magnetic stellar activity and open the way towards a better correction of the activity-induced stellar jitter.

In Paper III, we derived the astrometric time series corresponding to the observed solar activity pattern. Our activity model also allows us to produce the astrometric time series corresponding to the simulated activity pattern (apart of the photometric and RV ones). We will present our results in astrometry in a separate paper. This is justified as upcoming instruments such as the Nearby Earth Astrometric Telescope (NEAT, see e.g. Malbet et al. 2012) should allow for the first time the search for low-mass planets around nearby FGK stars with astrometry.

---

Acknowledgments

We would like to thank our referee (X. Dumusque) for his very useful comments on the manuscript.

References

All Tables

All Figures