© N. Meunier et al. 2019

1 Introduction

It is now well recognized that stellar activity strongly affects the detectability of exoplanets. First rough attempts to model the amplitude of this effect through radial velocity (RV) have been made with simple models that related simple activity coverage with jitter (Saar & Donahue 1997; Hatzes 2002; Saar et al. 2003; Wright 2005). Desort et al. (2007) modeled the RV that is caused by single spots. More sophisticated models describing the full behavior of the activity that causes the RV variations are needed, however, to estimate the effect of stellar activity more quantitatively and to test analysis and correction methods. Such models have been made for the Sun (Borgniet et al. 2015) and for a few configurations of other stars (Dumusque 2016; Dumusque et al. 2017). Other models have been proposed by Herrero et al. (2016); they reproduce contributions of spots and plages. Santos et al. (2015) modeled the contributions of spots alone.

We made a significant step when we modeled the solar RV and photometry using observed solar spots, plages, and network structures (Lagrange et al. 2010; Meunier et al. 2010a). This allowed us to show that the inhibition of the convective blueshift in plages dominates the long-term variability, which we validated by reconstructing the solar RV variation from MDI/SOHO (Michelson Doppler Imaging/SOlar and Heliospheric Observatory) dopplergrams (Meunier et al. 2010b). Direct or indirect (Moon, asteroids, Jupiter satellites) observations of the Sun later confirmed these results (Dumusque et al. 2015; Lanza et al. 2016; Haywood et al. 2016). We also studied the effect of activity on future astrometric measurements (Lagrange et al. 2011), which are important in the context of the current Gaia mission. Our second step was to generate similar time series based on randomly generated solar spots and plages, for which we used realistic properties over the solar cycle (Borgniet et al. 2015): this allowed us to study the effect of inclination, and to open the way to model stars other than the Sun. This is the objective of the present paper. Our approach focuses on using the proper spatio-temporal distribution of spots and plages, and on a physical relationship between spots and plages together with realistic physical properties. This is complementary to other approaches that focus on estimating finer details in contrast variations, for example (e.g., Cegla et al. 2018).

In this paper, we therefore extend the solar model described in Borgniet et al. (2015) to other stars. We propose consistent parameter sets to build RV, photometric, and astrometric time series. We also implement a model to describe the chromospheric emission as a function of time. The goal of such simulations is threefold: (1) to compare our model with observations for these different observables; (2) to help with the interpretation of these observations, and in particular to understand the degeneracies and biases well, as well as the effect of the various parameters (including “hidden” parameters such as inclination); and (3) to test correction methods and estimate the effect on exoplanet detectability through various techniques (RV, transits, and astrometry). Our objective in building the whole parameter set is to be as consistent as possible in the various choices so that we retain a large amount of the complexity of stellar variability while keeping the parameters to a reasonable number for this first set: all parameters that correspond to a given time series are compatible with each other.

The outline of the paper is the following. In Sect. 2 we explain how we adapt the solar parameters to other stars to generate spots and plages: this section is devoted to the procedures and laws we used to produce lists of spots and plages as well as their properties at each time step. In Sect. 3 we describe the required contrast and RV properties for producing the observables. We then present the chromospheric emission model in Sect. 4 as well as the calibration that must be made to produce realistic time series. Finally, we compare in Sect. 5 the obtained Log $R_{\text{HK}}^{\prime }$ values with what is expected from the input parameters, followed by a conclusion and a description of future works in Sect. 6.

Fig. 1 Variable parameters in our grid. Parameters with specific indications in orange depend on B–V and/or Log $R_{\text{HK}}^{\prime }$, and we also consider two or three laws, as indicated. A few parameters depend on specific conditions (plage contrast). The 2D grid in B–V and Log $R_{\text{HK}}^{\prime }$ is detailed in Fig. 3.

2 Generation of spots and plages

2.1 General principles

The model we used to produce spots and plages at each time step and to follow their evolution is described in detail in Borgniet et al. (2015). We summarize here the main parameter categories:

• spatio-temporal distribution of spots and plages: butterfly diagram and active longitudes;

• long-term variability: cycle length, amplitude, and shape;

• individual properties: size distributions, decay distributions, plage-to-spot size ratio, and plage-to-network decay;

• dynamics: rotation period, differential rotation, meridional circulation, and diffusion.

Plages are created at the same time as spots, then each type of structure follows its own evolution. Part of the plage decay creates network features. This leads to an entirely coherent model that describes spots, plages, and the network. In the following, unless otherwise mentioned (Sect. 2.7), the term plage refers to large plages and network structures, that is, all bright features, in order to simplify the presentation. We also recall that Meunier et al. (2010a) and Borgniet et al. (2015) considered plage sizes that were obtained from MDI data with a threshold of 100 G, leading to an adjustment of the various contrasts to match the observed solar photometric variations (while associated with the spot distribution). In this paper, we keep size distributions and contrasts that are consistent with that definition. The detailed parameters are described in Table E.1. The main differences with the model of Borgniet et al. (2015) are that we simplified the input spot number (see Sect. 4) and adapted the dispersion that was added to the shape of the reference cycle (see Sect. 2.6.3).

A large number of parameters were involved in our solar simulation. When we adapted these parameters to other stars, we did not have to explore the full space of possible parameters, as some parameters may depend on others. For example, the rotation period depends on the spectral type and on the activity level, so that for a given spectral type and activity level, the range of possible periods is limited. Empirical laws, sometimes with large dispersion that represents an actual variability between stars, have been established in the literature, allowing us to establish a correspondence between certain variables. These relations are sometimes multivariable. We use these laws in the following to build the parameter sets.

Our objective is to study the effect of important parameters on time series, for example, to establish which types of stars are most suitable for exoplanet detection, or to estimate the performance of correction methods in various conditions. Some parameters are not constrained at all for stars other than the Sun, therefore we keep some to the solar values in this work. The list of parameters that are different from the solar values is described in Fig. 1. A summary of the laws described in the next sections is also given in Table E.3. For five of these parameters, we used two (upper and lower bound) or three laws (median law as well) to cover a range of realistic values because we estimate that the observed variability is real. The remainder of this section is devoted to describing the way in which we derived all these parameters based on our current knowledge of stellar activity.

Fig. 2 First panel: Teff vs. spectral type. Second panel: B–V vs. spectral type. Third panel: stellar radius (in R⊙) vs. spectral type. Fourth panel: stellar mass (in M⊙) vs. spectral type.

2.2 Fundamental stellar parameters

The spectral type constitutes our first axis in the grid (the second is the average activity level, see the next section). We translate it into B–V values because many empirical laws used in the following are available as a function of B–V in the literature. We consider a wide range of stellar types, F6 to K4, that is, stars whose activity patterns are not very different from that of the Sun. For example, the convective blueshift, which has a critical effect on the RV amplitudes, has been estimated with a good precision by Meunier et al. (2017b) for this range of spectral types, but it is not well constrained beyond this range. The four laws we use are illustrated in Fig. 2 and the values are listed in Table 1.

• T_eff. Effectivetemperatures are derived from the spectral type using a fourth-degree polynomial from the observations of Gray et al. (2003). The validity domain is A2–K3, therefore we extrapolate this function over a small range for K4.

• B–V. are derived from T_eff (see above) with the law provided by Gray (2005).

• Radius. Many stellar radii have been measured using interferometry (Boyajian et al. 2012, 2013). We use Eq. (8) from Boyajian et al. (2012) to relate the radius to T_eff. This formula is valid for T_eff up to 5500 K, and we verified that the formula can be extrapolated from the measurements of Boyajian et al. (2013): the extrapolation is appropriate up to 6400 K, but with a larger dispersion.

• Masses. Stellar masses are derived from the radius in Boyajian et al. (2012, 2013) using a third-degree polynomial fit from their Tables 6 and 3, respectively.

2.3 Average Log R${}_{HK}^{\prime }$–(B–V) relationship

The activity level we consider here is the average Log $R_{\text{HK}}^{\prime }$ over time for a given star over time-scales of a few years because a star at a given age does not have a single Log $R_{\text{HK}}^{\prime }$. The average activity level constitutes our second axis.

2.3.1 Lower limit in Log R${}_{HK}^{\prime }$

We first estimate the lower limit for the Log $R_{\text{HK}}^{\prime }$. Several papers have estimated the average Log $R_{\text{HK}}^{\prime }$ versus B–V for large samples of stars and dates (Henry et al. 1996; Gray et al. 2003, 2006; Jenkins et al. 2008; Isaacson & Fischer 2010; Jenkins et al. 2011; Arriagada 2011; Schröder et al. 2012; Mittag et al. 2013). We first consider the B–V range from 0.45 to 0.94. Mittag et al. (2013) obtained a flat minimum S-index in this domain. This limit is not entirely strict because occasionally,a few stars lie below it, but this lowest flux derived from Mittag et al. (2013) is consistent with previous publications. Therefore we consider their value (corresponding to a S-index of 0.144) in the following to be the lower limit for activity in this B–V domain. This S-index can then be converted into a Log $R_{\text{HK}}^{\prime }$ value using the commonly used formula from Noyes et al. (1984a), as we do here.

Finally, we consider B–V above 0.94. A strong increase of the lower limit in activity (Isaacson & Fischer 2010; Mittag et al. 2013) corresponds to stars with a significant degree of activity, implying that low-activity stars are not observed in this B–V range. We have checked the HARPS spectra of stars close to this apparent limit, and they indeed show strong calcium emission. This lower limit can therefore be used to identify where stars are located in the 2D space (B–V, Log $R_{\text{HK}}^{\prime }$). In conclusion, the considered Log $R_{\text{HK}}^{\prime }$ values lie above the dashed line in Fig. 3 (solid for B–V below 0.94 as the two coincide in that domain).

Fig. 3 Upper panel: average Log $R_{\text{HK}}^{\prime }$ versus B–V representing our 2D grid (orange stars). The lower solid line shows the basal flux used in this paper, from Mittag et al. (2013) and Schröder et al. (2012) and the upper solid line shows the upper bound for the activity level versus B–V, from Lockwood et al. (2007). The dashed line corresponds to the minimum Log $R_{\text{HK}}^{\prime }$ we consider for stars with B–V above 0.94. The vertical dotted lines show the range in B–V covered by our simulations, the vertical dashed line vizualises B–V = 0.94, and the red horizontal line approximately corresponds to the Vaughan-Preston gap (see text). Lower panel: same for the S-index versus B–V.

2.3.2 Upper limit in Log R${}_{HK}^{\prime }$

We now consider the upper limit in Log $R_{\text{HK}}^{\prime }$. A first simple choice would be to consider a threshold from the Vaughan-Preston gap in the usually bimodal distribution of Log $R_{\text{HK}}^{\prime }$ values. Depending on the publication, the position of the gap ranges from −4.80 (Noyes et al. 1984a; Jenkins et al. 2011) to −4.6 or −4.7 (Mamajek & Hillenbrand 2008) with intermediate values (Wright et al. 2004; Jenkins et al. 2008; Henry et al. 1996; Gray et al. 2003, 2006). However, our purpose is to model stars with properties similar to solar properties in terms of plage-to-spot ratio, for example, at least given our current knowledge. We therefore used the results from Lockwood et al. (2007), which show this type of correlation versus B–V and Log $R_{\text{HK}}^{\prime }$. The interface between the spot-dominated regime (younger stars) and the plage-dominated regime (the older stars we are interested in) varies with B–V, and is about −4.5 for the most massive and −4.85 for the less massive stars. We use this as an upper limit for our Log $R_{\text{HK}}^{\prime }$ values (shown as the upper solid line in Fig. 3).

We could also have derived this upper limit from age isochrones (Mamajek & Hillenbrand 2008), but because we are more interested in the plage-to-spot ratio in our input parameters, this choice would be less pertinent. The age range covered by our simulations may therefore vary with B–V (see Sect. 2.4).

2.3.3 Log R${}_{HK}^{\prime }$ values between the lower and upper limits

Stars are observed with Log $R_{\text{HK}}^{\prime }$ between the lower and upper limits that we defined in the previous sections. The distribution of stars within that domain is not necessarily homogeneous, but this was not taken into account when we built the grid.

We considered Log $R_{\text{HK}}^{\prime }$ values higher than the lower level by 0.07 dex and then with a step of 0.05 dex up to the upper bound. Theses values are shown as orange stars in Fig. 3. They lead to 141 points in 2D space and correspond to the average Log $R_{\text{HK}}^{\prime }$ over time. For each of these positions, parameters were defined according to Fig. 1 and several time series were built. These parameters are described in detail in the remainder of this section.

2.4 Rotation period versus B–V and Log R${}_{HK}^{\prime }$

In this section we wish to determine which rotation rate (or range of rotation rates) to use to simulate a star of a given spectral type and average activity level. Several estimates of the rotation period as a function of the average activity level have been published using large samples of stars, either directly or through an estimate of the Rossby number. A comparison of these different laws is provided in Appendix A. We have then chosen to use the law relating the Rossby number and the average Log $R_{\text{HK}}^{\prime }$ from Mamajek & Hillenbrand (2008), with the estimated turnover time from Noyes et al. (1984a) to relate the rotation period and the Rossby number. Long periods are difficult to estimate, and samples are usually biased toward short periods. Laws are therefore uncertain for long periods, which correspond to our lower mass stars. We have then taken into account the observed dispersion around the law provided by Mamajek & Hillenbrand (2008), which we estimated from their data to be of about ±0.2 in Rossby number: from these upper and lower bound laws we derived the rotation period versus B–V for each Log $R_{\text{HK}}^{\prime }$, which we show in Fig. 4 as red and green curves, respectively. Although we took the observed dispersion into account, we might still underestimate the longest rotation periods: if this were the case, the effect on the final RV or photometric jitter is expected to be very small. However, it is expected to affect the frequency analysis in two ways: the power due to rotation will naturally be localized at a different period, and this may affect the morphology of the curves because the ratio between the rotation rate and the typical lifetime of the magnetic features will be different. The fact that simulations are always made for three different rotation rates will help to analyze the effect of our assumptions in future works.

Age is not a parameter in our simulation. However, we know that there is a relationship between rotation, activity level, and age (e.g., Wilson 1963; Skumanich 1972). For instance, when we use the laws of Mamajek & Hillenbrand (2008) for the most massive stars in our simulations, our range in Log $R_{\text{HK}}^{\prime }$ corresponds to ages between 0.5 and 3 Gyr. Lower mass stars in our simulations correspond to older stars, typically between 4 and more than 10 Gyr, depending on their average activity level.

Fig. 4 Chosen rotation periods vs. B–V for eight different Log $R_{\text{HK}}^{\prime }$ values between −5.1 and −4.75 (from top to bottom) from Mamajek & Hillenbrand (2008). The median law is shown in black, and extreme laws are shown in green and red.

2.5 Differential rotation and latitude coverage

The implementation of differential rotation is strongly related to the latitudinal extension over which magnetic activity is present because measurements of the differential rotation, based on the presence on active structures, only provide the differential rotation over that range in latitude and not the value corresponding to a full range of 0–90°. In this section we therefore discuss these two parameters together.

2.5.1 Differential rotation versus temperature

To derive a practical relation using our other input parameters, we used the differential rotation measured by Reinhold & Gizon (2015) from Kepler data for a very large sample of stars. As for the rotation period, we should keep in mind that observations are biased toward active stars, that is, fast rotators. Our objective is to define a law Ω(θ), where Ω is the rotation rate and θ is the latitude that can be used in our simulations. A parameter α is commonly defined from the minimum and maximum rotation periods given in Reinhold & Gizon (2015), P_min and P_max, $$\alpha =\frac{P_{\text{max}}-P_{\text{min}}}{P_{\text{max}}},$$(1)

which is a relative differential rotation. α is then available as a function of T_eff. We note that the Ω(θ) function for the Sun is usually described with three parameters (e.g., Snodgrass & Ulrich 1990) as Ω(θ) = Ω₀ + Ω₁sin²(θ) + Ω₂sin⁴(θ). Because the differential rotation for stars is much less well defined, we used only the first two coefficients in the following (Ω₂ = 0). When the following results are compared with solar differential rotation, caution is therefore advised.

We considered the stars in the sample from Reinhold & Gizon (2015) with log(g) between −3.94 and −4.94 as in some previous analyses of solar type stars (Das Chagas et al. 2016). For ten bins in T_eff in our range in temperature, we performed a linear fit between Log(α) and Log(P_rot), where P_rot is defined as $$P_{\text{rot}}=\frac{P_{\text{min}}+P_{\text{max}}}{2}.$$(2)

The two coefficients of these linear fits, p₀ and p₁, are shown in Fig. 5 as a function of T_eff. We then modeled p₀ and p₁ as a linear function of T_eff, which gives $$p_0=-3.485+{2.47810}^{-4}\times T_{\text{eff}}$$(3)

and $$p_1=1.597-{1.3510}^{-4}\times T_{\text{eff}}.$$(4)

For each T_eff and P_rot (previous section) in our grid, we can then derive α. We discuss differential rotation in more detail in Appendix B.

The effect of our choice of Ω(θ) is not expected to be critical for our simulations. When differential rotation is present, periodograms of the time series are thought to exhibit multiple peaks around the rotation period, which complicates estimating the rotation period, for example (in addition to the limited lifetimes of structures). The precise choice of Ω(θ) therefore mostly affects the complexity of the peak structures in the periodograms, but it does not affect the signal amplitude, for example.

Fig. 5 Upper panel: coefficient p0 vs. Teff from the fit of Log(α) = p0+p1×Log(Prot) for stars with log(g) between 3.94 and 4.94. Computations are made from the data published by Reinhold & Gizon (2015). Middle panel: same for p1. Lower panel:number of stars in each Teff bin.

2.5.2 Maximum latitude

We assumed(1) that structures are always present at low latitude at the end of the cycle (when the maximum latitude is higher, there might be no activity close to the equator either, as shown, e.g., by Işık et al. 2011, but these effects are not full understood so far), and (2) that the latitude coverage is directly related to the maximum latitude of the butterfly diagram (in the case of the Sun, we used an average latitude at the beginning of the cycle of 22°, with a possible extension of activity to 42°), hereafter θ_max. How θ_max varies with T_eff is not constrained. For lower mass stars, where the convective zone is thicker, we expect higher values of θ_max (e.g., Işık et al. 2011). On the other hand, for a shorter P_rot (in our case, for higher mass stars), we also expect larger θ_max (e.g., Schuessler & Solanki 1992). Because these two effects compete with each other, we do not know the proper trend from observationsor numerical simulations. Simulations such as the one made by Işık et al. (2011) are expected lead to some results in the future, but so far, the coverage in parameters is too sparse to conclude. As for the observation, the analysis of the data results of Reinhold & Gizon (2015) does not allow us to conclude either. Recent results using planetary transit across spot at the stellar surface allowed determining the latitudinal distribution for a specific star over a short period of time, see Morris et al. (2017), but the statistics is not yet sufficient to ascertain any trend.

Because the coverage in latitude of magnetic structures is not well constrained by either observations or observations, we considered in our simulations three possible levels for θ_max, remembering that we do not know how other stars differ from the Sun in that respect: the solar value itself θ_{max, ⊙}, θ_max,⊙+10°, and θ_max,⊙+20°. For each of these values of θ_max, P_rot (derived from the previous section and as defined in Eq. (2)) and α (derived as described above and defined in Eq. (1)) lead to Ω₀ and Ω₁, assuming that P_max corresponds to θ_max and P_min corresponds to the equator. This allows us to estimate the effect of these parameters on the time series.

2.5.3 Antisolar rotation

Numerical simulations have indicated that some stars probably present antisolar differential rotation (rotation is slower at the equator than at the poles). This would occur for stars with large Rossby number, that is, long P_rot and high masses (e.g., Brun et al. 2017). This is very difficult to observe from light curves (Santos et al. 2017), however, although there have been a few indications that this could be present: Reinhold & Arlt (2015) made tests on synthetic time series, applied their method to a small sample of 50 Kepler stars, and estimated that there was a possibility that the rotation in 10–20% of the stars is antisolar. When we compare our parameter grid in the (P_rot, stellar mass) space with the results from the magnetohydrodynamics (MHD) simulation of Brun et al. (2017), we find that fewer than 6% of our simulation stars (all are less massive than the Sun and are very quiet) may have such an antisolar differential rotation, although the threshold between the solar and antisolar regimes is not well defined. It is therefore still very uncertain, especially for stars with our parameter range, and probably does not concern many stars. Because this effect would not significantly affect our results, we consider only solar differential rotation here.

2.6 Cycle properties

The cycles of the stars we simulated are similar to the solar cycle in shape, although the amplitudes and ratio between maximum and minimum may be different. Stars with no variability will not be reproduced adequately, although the simulation with a very low cycle amplitude will present some similarities with such stars.

How many stars have a cycle is subject of debate. However, the existence of long-term variations is crucial for us here because these variations are critical for studying the effect on exoplanet detectability. Baliunas et al. (1998) analyzed Mount Wilson data and found that 15% of the stars had a constant activity level, 25% had a variability without any obvious periodicity (they did not show any smooth cycle like the Sun, but rather some erratic variations, and they correspond to young fast-rotating stars), and 60% had solar-like cycles. Lovis et al. (2011) were unable to find any period (defined as the period derived from the fit of a sine function on the data, even if a single cycle was observed) for 66% of a large sample of variable stars observed with HARPS, which is likely due to the sampling: stars without an identified period are very strongly biased toward a very poor sampling compared to the list of stars with an identified “cycle” (“cycle” here means that a proper fit with a sinusoidal was possible, not that it was repeating itself). When this bias is taken into account, the percentage of stars without a long-term variation similar to a “cycle” decrease to only 15–20%. More recently, results obtained from the analysis of the long-term Mount Wilson survey together with the Lowell survey (Hall et al., in prep.) show that 40% of the stars may have a relatively flat chromospheric emission over decades, although some of them show high chromospheric emission.

The statistics of the various stellar categories is therefore still uncertain. In practice, the lower cycle amplitude in our grid will allow us to cover almost no variability reasonably well, at least at low average activity level, because the ratio between the number of spots at cycle maximum and at cycle minimum will be able to reach values close to 1 in some cases. We do not attempt here to add more complexity to our time series, as this would represent additional parameters that are not at all constrained (e.g., what does the butterfly diagram look like when two cycle periods are present?), but future work will have to consider these configurations more precisely. Therefore our simulations are quite representative of stars that have some significant variability, except for the more complex stars, as well as of stars with very low variability.

2.6.1 Cycle period

We compared the cycle period versus rotation rate from various sources (Baliunas et al. 1996; Saar & Brandenburg 1999; Böhm-Vitense 2007; Oláh et al. 2009, 2016; Suárez Mascareño et al. 2016) for the so-called inactive branch when relevant. Stars on this branch are old stars similar to the Sun as in our model, and not young active stars. Except for Böhm-Vitense (2007), which lies apart, these stars provide a coherent picture. The slope of Log(P_cyc/P_rot) versus Log(1/P_rot) is in the range 0.74–1.09, and we consider here an average between these different sources, that is, a coefficient of 0.84: this gives the cycle period we show in Fig. 6. The curve is relatively flat. However, the dispersion in the observations is likely to be real, and the dashed lines shows the two extreme laws that we also considered to account for the observed variability. We therefore explored a wide range of cycle periods.

Fig. 6 Pcyc (in years) vs.Prot used in our grid (upper bound as solid line, lower bound and intermediate range as dashed lines).

2.6.2 Cycle amplitude

Several studies produced amplitudes for the cycle period, especially as a function of Log $R_{\text{HK}}^{\prime }$. We have compared the laws from various sources: Radick et al. (1998), Saar & Brandenburg (2002), Lockwood et al. (2007), Hall et al. (2009), and Lovis et al. (2011). When the observed dispersions are taken into account, the agreement between them is good overall, except for the existence of very large amplitudes in Lovis et al. (2011); we discuss this below. The trend for large amplitudes for larger Log $R_{\text{HK}}^{\prime }$ is globally weak. The values obtained by Hall et al. (2009) seem to be slightly lower than in other studies. It is also important to note that in general, these samples contain very few quiet stars with Log $R_{\text{HK}}^{\prime }$ below −5.0, so that the amplitudes in this domain are not well constrained (and are also most likely to be affected by noise).

Because the very large upper bound derived by Lovis et al. (2011) is very puzzling, we verified the temporal variability of all stars whose A_cyc (half-amplitude in ${R^{\prime }}_{\text{HK}}\times {10}^5$) was larger than 0.3 in their sample. To do this, we used HARPS archive data. We find that with extended observations since 2011, all of them fall below 0.33, which agrees very well with the other publications.

Figure 7 shows A_cyc from Lovis et al. (2011) versus B–V and average Log $R_{\text{HK}}^{\prime }$. Stars of different spectral types have different A_cyc. We used the boundaries indicated in the figure to derived three laws: an upper value, a lower value, and an intermediate between the two for each point of the grid in (B–V, Log $R_{\text{HK}}^{\prime }$). We chose a minimum A_cyc of 0.005: this corresponds to stars with very low variability.

2.6.3 Cycle shape

As discussed at the beginning of this section, we considered cycles similar to the solar cycles. The chosen shape of the cycle is the shape of thelast solar cycle, as shown in Fig. 8. The ratio between maximum and minimum can be different, however. At each time step, some random variability was added to that curve (25% of the amplitude at that time step) to represent the stochastic variability that can be introduced by the dynamo in terms of flux emergence. This amplitude is somewhat arbitrary, but it gives a final realistic dispersion for the Sun. The input parameters, in addition to this shape and dispersion (constant), are therefore the minimum and maximum level in Log $R_{\text{HK}}^{\prime }$, which must correspond to the average Log $R_{\text{HK}}^{\prime }$ we wish to obtain.

Fig. 7 Upper panel: half full amplitude of stellar cycles vs. B–V, derived from Lovis et al. (2011) after revision of the largest amplitudes (see text) for different types of stars: B–V < 0.7 (black stars), 0.7 < B–V < 0.9 (red squares), and B–V > 0.9 (green triangles). The black lines correspond to the lower and upper boundaries that were taken into account in building the grid. Upper panel: same vs. average Log $R_{\text{HK}}^{\prime }$. The solid line represents the lower limit that was taken into account in building the grid.

2.7 Small-scale convection level and convective blueshift

The inhibition of the convective blueshift in plages is an important contribution to the final RV. This parameter has no effect on the other observables. We used the study of the activity effect on convective blueshift based on a large sample of F-G-K stars that was presented in Meunier et al. (2017b). We found that the convective blueshift depended not only on B–V, but also on Log $R_{\text{HK}}^{\prime }$. For several B–V bins, we extrapolated these convective blueshifts to the basal Log $R_{\text{HK}}^{\prime }$ (i.e., the level of convective blueshift we would have if no activity were present). This gives the dashed line in Fig. 9. The convective blueshift was estimated as in Meunier et al. (2017a), that is, it is based on the Sun from Reiners et al. (2016), with 355 m s⁻¹ for the solarconvective blueshift.

After we derived the convective blueshift, we applied an attenuation factor (which provides the amplitude of the RV in plages and networkstructures) and a correction factor for projection effects (considering effects perpendicular to the surface, as in our previous work). In previous work, we used an attenuation factor of two-thirds based on Brandt & Solanki (1990), but also a smaller solar convective blueshift. Because our amplitude led to good results when we compared out results with a solar reconstruction of the long-term RV variation (Meunier et al. 2010b), we would need to use a smaller attenuation factor (0.38) to obtain the same results, given our new convective blueshift: this is what was used in our simulations, which gives the local ΔV applied to each structure as a function of B–V, shown as the solid line.

Another result obtained by Meunier et al. (2017b) is a possible trend versus T_eff for the attenuation factor, which would imply a correction factor of −2.077+5.324 10⁻⁴ T_eff. The effect of thistrend is shown in Fig. 9 as the dotted line. However, this trend is poorly constrained below 5300 K, and we therefore chose to make our simulation without this factor. The effect on the resulting RV can be estimated during analysis, as this correction can be applied afterward to the time series.

Finally, the convective blueshift is higher in larger structures. We implemented the dependence between the velocity and the sizederived in Meunier et al. (2010b). There is typically a ratio of 6 between the largest and smallest structures.

Fig. 8 Wolf number vs. phase smoothed over the last solar cycle that was used as a reference.

Fig. 9 Convective blueshift derived for the basal Log $R_{\text{HK}}^{\prime }$ (dashed line), adapted from Meunier et al. (2017b), and scaled with the solar values derived in Meunier et al. (2017a). The solid line shows the local convective blueshift we used in our model, after correction for projection effects for a constant attenuation of the convective blueshift in plages of 0.38 (see text). The dotted line shows the same parameter when the trend of the attenuation factor vs. Teff is taken intoaccount.

Fig. 10 Upper and lower bounds for the difference between spot and photosphere temperature vs. Teff from Berdyugina (2005). The lower bound is the solar value derived in Borgniet et al. (2015).

2.8 Spot temperature

The temperature of stellar spots remains poorly constrained because it is very difficult to measure (and strongly degenerated with spot size in photometric light curves). We used the results of Berdyugina (2005), which show a trend with T_eff (lower spot contrast for lower stellar T_eff), and a large dispersion because stars (including the Sun) exist at a much lower temperature contrast. The trend and order of magnitude have recently been confirmed in numerical spot simulations (Panja et al., in prep.). We therefore used two laws that represent two extreme configurations, assuming that stars have temperatures within that range. This is shown in Fig. 10: the lower boundary law corresponds to the solar contrast we used in Borgniet et al. (2015), and the upper boundary is derived from Berdyugina (2005). The computations were made at 6000 Å, as in Borgniet et al. (2015), and the contrasts were adjusted to correspond to the bolometric photometric variability (see Borgniet et al. 2015, for details).

2.9 Plage contrasts

Contrasts of stellar plages are poorly constrained as well. In our previous work for the Sun (Meunier et al. 2010a), we used a law that described a temperature contrast versus μ (cosine of the angle between the line of sight and the local vertical at the solar surface) similar to the law described in Unruh et al. (1999). We then adjusted this slightly (together with the contrast of the spot temperature that we described in the previous section) to fit the observed solar irradiance. Borgniet et al. (2015) used a description in terms of intensity contrast (plage intensity divided by quiet-Sun intensity, minus one), which was described as a second-degree polynomial in μ. A similar adjustment was made to fit the photometry (which was necessary because we used a slightly different center-to-limb darkening function).

In the present paper, we not only need to describe the contrast as a function of μ for the Sun,but also for other spectral types. We used the results from the MHD simulations performed with the MURAM code by C. Norris (Norris et al. 2016, 2017; Norris 2018). She provided coefficients describing the plage intensity for different magnetic field levels (0, 100, and 500 G) versus μ (in the 0.2–1 range, not available for μ below 0.2) for G2, K0, and M0 stars as described in Norris (2018). These intensities where computed for the HARPS wavelength range. The contrasts take slightly different shapes depending on the parameters, but the global trend is a higher contrast for stronger magnetic fields and higher T_eff.

Different functions can be used to describe the intensity variations versus μ (e.g., Yeo et al. 2013), and they provide different values in the 0–0.2 μ range: however, their effect on the final RV is very low (lower than 0.1%), therefore we use a quadratic form in μ in the following. We then compute the contrasts as quadratic functions of μ for these three spectral types and two magnetic field levels and interpolate (or extrapolate for stars in the F6–G1 range) for other spectral types and different magnetic fields.

Our simulations provide sizes. We therefore established a law relating the size (A) and the magnetic field flux (B) for the plages and magnetic features we are interested using MDI/SOHO (Scherrer et al. 1995) magnetograms that cover a full solar cycle, which gives $$\text{Log}(B)=2.1134+0.1355\times \text{\hspace{0.17em}Log}(A),$$(5)

where B is in G and A in ppm of the hemisphere. For each structure in the simulation, we therefore computed its associated magnetic field according to this law, to be able to interpolate (all values are between 100 and 500 G).

When we apply this procedure to a G2 star with an activity level similar to that of the Sun, we find that the contrasts are slightly higher on average than those of Borgniet et al. (2015), they are higher by a factor 1.5. We therefore divided the contrasts by this value for all stars for consistency with our definition of structures sizes and the good agreement with the solar irradiance variability (and the corresponding definition of spot size and spot contrast).

2.10 Other parameters

• Spatio-temporal distribution. Latitude and longitude distributions as well as north-south asymmetry and active longitude parameters were kept to the solar values from Borgniet et al. (2015). The latitude at the beginning of the cycle (related to the maximum possible latitude) was also tested with different values (see Sect. 2.5.2). Migration was considered to be equatorward as for the Sun, although there have been indications that poleward migration could exist (Messina & Guinan 2003; Moss et al. 2011), but this is not well constrained.

• Large-scale dynamics. The differential rotation discussed in Sect. 2.5 that we adapted to each grid point is described with only two coefficients (instead of three for the Sun). The same law was used for all structures (spots, plages, and network). The meridional circulation was kept to the solar value used in Borgniet et al. (2015) based on Komm et al. (1993), and was also the same for all structures.

• Spot properties. The distributions of spot size and decay rate were kept to the solar values used in Borgniet et al. (2015), which were adapted from Martinez Pillet et al. (1993), Baumann & Solanki (2005), Lagrange et al. (2010), and Meunier et al. (2010a).

• Faculae properties. The faculae properties were similar to those used in Borgniet et al. (2015), which included the ratio distribution between plage and spot sizes, and decay rates (here we considered the plages that were produced each time we generated a spot).

• Network properties. Network properties were similar to those used in Borgniet et al. (2015), which include the diffusion coeficient from Schrijver (2001), the fraction of plage flux that was used to build the network, and the decay rate. The diffusion coefficient was scaled with the amplitude of the convective blueshift, as discussed in Sect. 2.7.

3 From structures to observables

3.1 Filling factor and photometric and radial velocity time series

Because of the huge number of simulations, it is not possible to compute the observables as done in our previous work (Lagrange et al. 2010; Meunier et al. 2010a; Borgniet et al. 2015), in which we computed maps, then spectra, and finally RVs (and other observables) in the same way as for stellar observations. We therefore simplified the computations as follows: we directly summed the contribution from each structure to the RV, photometry, and astrometry. The filling factors were also computed. This implies that we assumed that the structures are point-like, which means that we did not need to check for superimpositions. This is different to what was done in Borgniet et al. (2015). This also means that we neglected any geometrical effect that would be due to a large area covered by a structure, which is a good assumption because we modeled relatively quiet stars with moderate structure sizes. The assumption is very good for spots because the maximum size of a spot would correspond to a radius of about 3° (which would be a very rare case; typical spots are much smaller, see Table E.1). Plages may cover a larger area, which is expected to add a second-order distortion to an ideal time series: in most cases, the extreme sections of a structure (east and west) are expected to produce a signal that is very similar to the central part of the structures on average; in addition, we considered at each time step an irregular decay of the structures that has a random amplitude, so that the distortion produced by the assumption would not be identifiable given this other source of irregularity. The same structures were used for all stellar inclinations. The formulae are provided in Appendix C.

We recall that we chose to study different laws for several of our parameters (as summarized in Fig. 1), so that several time series of spots, plages, and network features were produced for each of the 141 points of the 2D grid in (B–V, Log $R_{\text{HK}}^{\prime }$), each corresponding to a different parameter set (Table E.1). This leads to 11 421 time series, or 22 842 when the two levels for Tspot are considered for each inclination and observable, hence a total of 228 420 realizations for each observable. Inclinations take values between 0° (pole-on) and 90° (edge-on), with a step of 10°. Figure E.1 shows a summary of all parameters.

3.2 Temporal sampling and duration

Time series must have a sufficiently long duration to allow us to test analysis methods, and they must cover at least a cycle period. To keep it reasonable, however, we imposed a maximum of 15 yr (which is just above the maximum cycle period we considered). We then simulated an integer number of cycles, choosing the maximum number that would lead to a duration shorter than 15 yr.

The time step was one day on average, but as in Borgniet et al. (2015), we added a small random departure (within ±4 h) from the regular sampling to mimic a realistic sampling.

3.3 Addition of short-term variablity in RV

3.3.1 Principle

To produce realistic RV time series that include all contributions at various timescales, we added the contribution of oscillations, granulation, and supergranulation, as was done in Dumusque (2016). We call these three contributions to RV the OGS signal hereafter. The principle is the following. For each spectral type and each variable, we computed one time series that covered 15 yr and had a time step of 30 s. From this, the time series corresponding to a given sampling can be extracted. To produce a time series like this, we computed the inverse Fourier transform of the power spectrum as a function of the frequency ν, P(ν), for each of these contributions. The parameters describing the power depend on the spectral type. We also computed smoothed time series (with a bin of one hour) to simulate the effect of long exposure times. In the following, we mostly use such long-time exposure time series, assuming some good observing conditions. The series with no smoothing may be used for a comparison with observations that were made with short exposure times, however. In practice, a long-duration (15 yr) time series was produced for each spectral type, with a time step of 30 s; for a given time series in our grid, we extracted either the instantaneous values or the one-hour average corresponding to the same sampling.

We can also add a white noise of 0.6 m s⁻¹ to simulate instrumental noise (this value corresponds to typical uncertainties on individual measurement from HARPS data for G stars). We considered four types of time series in our analysis:

• 1: original RV time series caused by magnetic activity. These are useful to study only the activity contribution.

• 2: original RV time series plus oscillation, granulation, or supergranulation signal (no smoothing). This is useful to compare with observations under ideal conditions (assuming the instrumental noise was totally corrected for).

• 3: original RV time series plus oscillation, granulation, or supergranulation signal (no smoothing) plus instrumental noise.

• 4: original RV time series plus oscillation, granulation, or supergranulation signal (one-hour smoothing) plus instrumental noise.

3.3.2 Oscillations

We used thefollowing power function for the oscillations: $$P(\nu )=A\times {\text{e}}^{-{(\nu -{\nu }_0)}^2/2/{\Gamma }^2},$$(6)

which describes the mode envelopes and not the individual modes themselves (e.g., Kallinger et al. 2014). We adopted the following scaling laws for different types of stars: $$A={(T_{\text{eff}}/T_{\text{eff},\odot })}^4\times {(R/R_{\odot })}^2/{(M/M_{\odot })}^{0.7}$$(7)

from Kjeldsen & Bedding (1995), with the adaptation of Samadi et al. (2007) for the exponent, where A is relative to the solar value, $${\nu }_0=(M/M_{\odot })/{(R/R_{\odot })}^2/\sqrt{T_{\text{eff}}/T_{\text{eff},\odot }},$$(8)

from Bedding & Kjeldsen (2003), where ν₀ is relative to the solar value, and $$\Gamma =\sqrt{(M/M_{\odot )}}/{(R/R_{\odot })}^{1.5}$$(9)

from Kippenhahn & Weigert (1990) and Belkacem et al. (2013), where Γ is relative to the solar value.

These laws were scaled with the following solar values: A_⊙ = 200 (m s)⁻² Hz⁻¹ (which provides an amplitude of the power that agrees well with the observed power, e.g., from Davies et al. 2014), ν_0,⊙ = 3140 × 10⁻⁶ Hz, and Γ_⊙ = 361 × 10⁻⁶ Hz (both from Kallinger et al. 2014).

3.3.3 Granulation

We used thefollowing power function for the granulation signal: $$P(\nu )=A/(1+{(\tau \nu )}^{\beta })$$(10)

from Harvey (1984). The power spectrum in RV of a few stars was analyzed (Dumusque et al. 2011), but the sample is not large enough to derive a proper trend. We therefore used a scaling derived from the numerical simulation of Beeck et al. (2013) to obtain the scaling of A and τ by fitting a linear law on their results: $$A=(0.3+6.323\times {10}^{-4}\times (T_{\text{eff}}-4594)),$$(11)

and $$\tau =-2.831+1.574\times {10}^{-3}\times T_{\text{eff}}$$(12)

are normalized by the same amplitude and τ, respectively, computed for 5784 K (G2). β was kept to the solar value. The solar values were derived from a fit on the simulated time series produced in Meunier et al. (2015): A_⊙ = 154 (m s)⁻² Hz⁻², τ_⊙ = 2781 s, and β_⊙ = 1.97.

3.3.4 Supergranulation

The formula for the supergranulation power is similar to Eq. (10) for granulation. Supergranulation seems to be present in stars other than the Sun (Dumusque et al. 2011), but the statistics is not sufficient to describe theparameters as a function of spectral type. Given the lack of knowledge on stellar supergranulation and becauseit is likely to be related to the granulation pattern and amplitude (e.g., Roudier et al. 2016), we used the granulation scaling relation. We considered the solar time series simulated in Meunier et al. (2015), which correspond to intermediate parameters between the two extremes the authors evaluated (supergranulation is less strongly constrained than granulation), and then fitted them, which gives the following parameters: A = 43 000 (m s)⁻² Hz⁻¹ and τ = 1.1 × 10_⁶ s. β was kept fixed to the granulation value.

3.3.5 RVjitter caused by the OGS signal

The rms RV produced by the OGS signal alone is shown in Fig. 11 as a function of B–V. It ranges from 1.8 to 1.1 m s⁻¹ for stars between F6 and K4. After averaging over one hour, the values are lower than 1 m s⁻¹. They lie between 0.9 and 0.5 m s⁻¹ from F6 to K4.

Fig. 11 RV jitter due to the oscillation, granulation, and supergranulation vs. B–V for instantaneous values (black) and averaged over one hour (red).

4 Chromospheric emission – calibrating the spot number

4.1 Objectives and principle

The laws described in Sects. 2 and 3 depend on Log $R_{\text{HK}}^{\prime }$, which is an observable. However, the input parameters of our simulations are the number of spots (Sect. 2.1). We therefore need to know how many spots to inject if we wish to reach a given activity level that is described by its average Log $R_{\text{HK}}^{\prime }$. For this purpose, we need two elements: a chromospheric emission model that uses the input parameters of our simulations, and a calibration law relating Log $R_{\text{HK}}^{\prime }$ and spot number. This model will also be very useful after the simulations are performed to determine to which Log $R_{\text{HK}}^{\prime }$ they correspond as well: the exact average Log $R_{\text{HK}}^{\prime }$ of simulation may be slightly different from the one in Fig. 3, but our objective is that it should be close (so that the input parameters we have used for that particular simulation are valid).

Assuming a solar chromospheric emission model, we therefore performed a series of simulations with constant activity levels, which will give the typical contribution per injected spot to the S-index. This was made for a constant number of spots¹ on eight-year time series. The average was then made on inclination (because the same structures were used for all inclinations, they correspond to the same inputs in terms of spot number), θ_max (very small variation), and spectral type (no trend observed, which is expected because we used the same law to compute the plage contribution). We note that this calibration depends on the plage-to-spot size ratio: we kept it constant in this paper, but if it were to vary, new calibrations are required. The same is true if the size distribution of spots changes.

4.2 Principle of the chromospheric emission model

A necessary step is therefore to implement a model providing an S-index (and then Log $R_{\text{HK}}^{\prime }$) based on a list of plages and network features at each time step. The full model is described in Appendix D. We provide here the general principles. The model is based on the work of Meunier (2018) for the Sun, and includes three contributions: (1) the basal flux (when no activity is present, determined in Sect. 2.3.1); (2) the contribution of plages and network structures, with a law that depends on their size, following Harvey & White (1999); (3) the contribution of the quiet star (“quiet” here means outside active regions and network, hereafter QS). This last contribution is important because it must vary from one star to the other: if it were kept constant, it would be impossible to observe stars with variability while having an average activity level below that of the quiet Sun because the level at solar minimum (corresponding to no structures) is significantly above the basal level. This is discussed in more detail in Appendix D. We propose that this contribution depends on the average activity level of the star: the more active the star, the stronger the (weak) magnetic field in the quiet star, and the larger the contribution of the quiet star to the chromospheric emission. The exact choice of the QS contribution will affect the number of spots, and most especially, the number at cycle minimum.

For a given Log $R_{\text{HK}}^{\prime }$, we computed the S-index, from which the basal flux was removed, as well as the QS contribution for that Log $R_{\text{HK}}^{\prime }$ (see Appendix D.1). The resulting flux was divided by the typical flux per structure, providing the number of spots to inject to obtain the Log $R_{\text{HK}}^{\prime }$ we wish to reach. The resulting calibration works very well, as we show in the next section.

4.3 List of the different time series

Table 3 shows the list of time series that we produced during the simulations. We recall that except for the number of structures, plage refers to all bright structures, from large structures in active regions to the smallest structures in the network.

5 Log R${}_{HK}^{\prime }$ behavior

5.1 Comparison between objective and realization

As explained in Sects. 2 and 4, we wish to simulate time series with a given average activity level, as well as a certain cycle amplitude (in Log $R_{\text{HK}}^{\prime }$), that is, Log ${R^{\prime }}_{\text{HK},\text{obj}}$ and ΔLog ${R^{\prime }}_{\text{HK},\text{obj}}$. A calibration was necessary to achieve this goal (Sect. 4.2), therefore we must check that the simulations behave as planned. We computed the average Log ${R^{\prime }}_{\text{HK},\text{out}}$ for each time series, where “out” stands for the output Log ${R^{\prime }}_{\text{HK}}$ from the simulations. The properties of this realized Log ${R^{\prime }}_{\text{HK}}$ (“out”) time series was compared tothe targeted Log ${R^{\prime }}_{\text{HK}}$ (“obj”) fromthe grid of parameters from Sect. 2.3. The amplitude, ΔLog ${R^{\prime }}_{\text{HK},\text{out}}$, was computed using a sinusoidal fit of each smoothed Log $R_{\text{HK}}^{\prime }$ time series.

The two upper panels (A and B) in Fig. 12 compare the Log ${R^{\prime }}_{\text{HK},\text{obj}}$ and Log ${R^{\prime }}_{\text{HK},\text{out}}$. We observe a strong inclination effect on Log $R_{\text{HK}}^{\prime }$; departures from the average are within 0.05 dex. However, on average, the difference between the expected value and the final value is much smaller when this inclination dependence is removed (below 0.01 typically). The differences are smaller than the typical uncertainties on Log $R_{\text{HK}}^{\prime }$ values as estimated by Radick et al. (2018), of about 0.06 dex, and our Log $R_{\text{HK}}^{\prime }$ should not be considered to be more precise than this in absolute value (although for a given simulation, the relative variability will be much more precise, of course).

In Appendix D, we mention the possibility of a trend in chromospheric emission versus T_eff. We did not include this trend because it is still very uncertain. When it is applied, the difference is small; the largest difference for our lowest mass stars is about 0.03 dex. If it is real, we would need fewer spots for the stars with the lowest mass than are included in the present simulations to reach a given objective because the emission for a given plage would be higher.

A fit of the time series with a sinusoidal provides an estimate of the amplitude of the cycle and of its period, which can be compared to expectations. Panel C in Fig. 12 shows the inclination effect on the cycle amplitude, while panels D and E in Fig. 12 compare ΔLog ${R^{\prime }}_{\text{HK},\text{obj}}$ and ΔLog ${R^{\prime }}_{\text{HK},\text{out}}$ for average inclinations. The cycle amplitudes from the simulations are also very close to the expected ones. Finally, the last plot (panel F) compares cycle periods, which for most simulations are in good agreement. There are a few outliers, but these are mostly due to low-amplitude simulations, for which the measurement itself is not reliable.

We conclude that after averaging, the average Log $R_{\text{HK}}^{\prime }$, the amplitude, and period of the cycle agree with the input parameters. The inclination effect is discussed in the next section.

Fig. 12 Panel A: Log ${R^{\prime }}_{\text{HK},\text{out}}$–Log ${R^{\prime }}_{\text{HK},\text{obj}}$ vs. Log ${R^{\prime }}_{\text{HK},\text{obj}}$ for all simulations. The color indicates stellar inclination from 0° (pole-on, yellow) to 90° (edge-on, blue). Panel B: Log ${R^{\prime }}_{\text{HK},\text{out}}$ averaged over all inclinations vs. Log ${R^{\prime }}_{\text{HK},\text{obj}}$. The solid red line indicates the y = x linear function. Panel C: ΔLog ${R^{\prime }}_{\text{HK},\text{out}}$–ΔLog ${R^{\prime }}_{\text{HK},\text{obj}}$ vs. Log ${R^{\prime }}_{\text{HK},\text{obj}}$, where Δ represents the amplitude of the cycle, for different inclinations (the color code is similar to panel A). Panel D: ΔLog ${R^{\prime }}_{\text{HK},\text{out}}$–ΔLog ${R^{\prime }}_{\text{HK},\text{obj}}$ vs. Log ${R^{\prime }}_{\text{HK},\text{obj}}$, where Δ representsthe amplitude of the cycle, after averaging the simulations made for inclination of 40° and 50°. Panel E: ΔLog ${R^{\prime }}_{\text{HK},\text{out}}$ vs. ΔLog ${R^{\prime }}_{\text{HK},\text{obj}}$, after averaging for inclination of 40° and 50°. The solid red line indicates the y = x linear function. Panel F: cycle period (in years) derived from a fit on smoothed Log $R_{\text{HK}}^{\prime }$ time series vs. the prescribed cycle period.

5.2 Dependence on inclination

We observe a strong inclination effect on the average Log $R_{\text{HK}}^{\prime }$, with a stronger value for an edge-on than for a pole-on configuration. This is in agreement with the results of Shapiro et al. (2014), which were based on simulations with a simpler model of the chromospheric emission (no structures or size dependence). Knaack et al. (2001) obtained a much weaker dependence, probably because they used a model that did not take all parameters into account. The inclination effect is also strong on the long-term amplitude in Log $R_{\text{HK}}^{\prime }$, although it presents a large dispersion. For the solar θ_max, the amplitude is larger for the edge-on configuration, with a difference of about 20–40% depending on the simulation. This has been observed by Knaack et al. (2001). For larger θ_max, the difference is smaller, with a slight predominance of larger amplitude when edge-on for θ_max,⊙+10° and a reversal for θ_max,⊙+20° (with a large dispersion and difference occasionally up to 20%).

5.3 Example of RV and Log R${}_{HK}^{\prime }$ time series

Figure 13 shows a few examples of time series for a small sample of spectral types, activity levels, and inclinations in chromospheric emission and radial velocity. The different activity levels correspond to different ages. When the law derived by Mamajek & Hillenbrand (2008) to relate rotation and age is assumed, the first two panels would correspond to an age of 3.3 Gyr, thethird panel to 4.3 Gyr, and the fourth to 8.2 Gyr. This shows a good similarity between RV and Log $R_{\text{HK}}^{\prime }$ long-term variations. It is also possible to obtain relatively flat curves for the lowest cycle amplitudes.

6 Conclusion

We have proposed a model to produce realistic time series of different variables (RV, photometry, astrometry, and chromospheric emission) that represent complex activity patterns for a wide range of stars. We have described the model in detail: a specificity of our simulations is that we use consistent parameter sets for a wide range of stars, that is, old F6-K4 star with different activity levels.

Our very large set of time series will be analyzed in detail in subsequent papers. We will compare the RV jitter between simulations and observations, analyze the effect of parameters on the RV jitter, and use this RV jitter to predict the detectability of exoplanets as a function of B–V and Log R_HK (Meunier &Lagrange 2019a). The detailed relationship between RV and Log $R_{\text{HK}}^{\prime }$ will be studied in order to understand why the corrections of RV time series using a linear function of Log $R_{\text{HK}}^{\prime }$ are limited (Meunier & Lagrange 2019b). RV times series will be further analyzed to produce detection limits by taking the frequency behavior of the stellar variability into account. Finally, a similar analysis will be made for astrometry. The effect of oscillation, granulation, and supergranulation will also be studied in more detail as a function of spectral type and activity level.

These time series are a good tool to provide clues to help interpret stellar variability from brightness time series: this is crucial because there are many degeneracies and biases, and synthetic time series are useful to determine the effect of the different parameters. They can also be used to test new methods, not only a correcting method for purposes of exoplanet detection, but of stellar activity analysis.

When these globally consistent parameter sets are built, the main limitation in our opinions is the poorly constrained QS contributionto the chromospheric emission. We have made two strong assumptions because our knowledge is incomplete. First, we have neglected the variation of this contribution with time, although we expect a small variation with a complex pattern from our solar study Meunier (2018; competition between stronger magnetic field at cycle maximum, but also a lower surface coverage). Second, we imposed that the number of spots at cycle minimum varies within a small range (and is small). We cannot exclude that for very active stars, for example, there could be a trend of a lower QS contribution and larger spot number at cycle minimum. This is not constrained, however, although we know that the Sun, which lies in the middle of our grid, has very few spots at cycle minimum and therefore conforms to our assumption. The only way to go beyond this limitation would be to better understand this QS contribution over the cycle and for different activity levels, most likely from dedicated simulations (MHD or using flux tubes down to very small spatial scales).

Fig. 13 Firstline: Log $R_{\text{HK}}^{\prime }$ (left panel) and RV (right panel) vs. time for a moderately active G2 star and an inclination of 90°. Second line: same for an inclination of 0°. Third line: same for a quiet G2 star and an inclination of 90°. Fourth line: same for a moderately active K2 star and an inclination of 90°.

Acknowledgements

This work has been funded by the Université Grenoble Alpes project called “Alpes Grenoble Innovation Recherche (AGIR)” and the ANR GIPSE ANR-14-CE33-0018. We are very grateful to Charlotte Norris, who has provided us the plage contrasts used in this work prior the publication of her thesis. We thank Frédéric Baudin and Sasha Brun for useful discussions about some aspects of this work. We thank Pascal Rubini for his work of convert the simulation code into C++. We thank the anonymous referee for his/her useful comments which helped to improve the paper.

Appendix A Rotation period

Figure A.1 shows a comparison of the rotation period derived from observations presented in different sources. The first two plots show the rotation period versus B–V for two values of Log $R_{\text{HK}}^{\prime }$ (−5.10 and −4.75, respectively) from Mamajek & Hillenbrand (2008), Noyes et al. (1984a), and Saar & Brandenburg (1999). Mamajek & Hillenbrand (2008) and Noyes et al. (1984a) are quite close to each other. Saar & Brandenburg (1999) give rotation periods longer than all the others for the most quiet stars, as illustrated in the last plot, but these are always poorly constrained.

Fig. A.1 First panel: Comparison of the rotation period vs. B–V for Log $R_{\text{HK}}^{\prime }$ = −5.1 from various sources: the solid black line is from Mamajek & Hillenbrand (2008), the green dot-dashed line from Noyes et al. (1984a), and the dashed red line from Saar & Brandenburg (1999). Second panel: Same for Log $R_{\text{HK}}^{\prime }$ = −4.75. Third panel: Similar comparison vs. Log $R_{\text{HK}}^{\prime }$ for eight levels in B–V between 0.5 and 1.1. The vertical line approximately shows the lower limit for the validity of the laws (the shaded area indicates the zones where it is not valid).

Other papers also provide rotation periods for very large samples, for example, from Kepler light curves (Nielsen et al. 2013; McQuillan et al. 2014; García et al. 2014; Reinhold & Gizon 2015), but they depend on the photometric variability, which cannot be translated directly into an average Log $R_{\text{HK}}^{\prime }$ level. Others are only given as a function of magnetic fields (Vidotto et al. 2014) or without any indication of activity level (Strassmeier et al. 2012). They can therefore not be easily included in our set of parameters as such.

Appendix B Differential rotation and maximum latitude

Reinhold & Gizon (2015) compared their measured differential rotations with previous results. They obtained a good agreement with laws obtained by Hall (1991) and Donahue et al. (1996). Barnes et al. (2005) and Collier Cameron (2007) obtained much steeper laws versus T_eff, however. We note that Das Chagas et al. (2016) obtained a differential rotation that is very similar to the solar one from the analysis of 17 Kepler solar-type stars. The results obtained by Balona & Abedigamba (2016) are more difficult to compare because of their normalization. Most numerical simulations of stellar differential rotation cover a small range in parameters, either have very fast rotation periods (e.g., Küker & Rüdiger 2011), or they are too close to the solar case (e.g., Küker et al. 2011). The simulations of Brun et al. (2017) cover a wider range (G and K stars): they derive a scaling law of ΔΩ versus mass and Ω. For their solar-type differential rotation, they obtained a scaling as M^0.73Ω^0.66, for ΔΩ between equator and 60°. They attempted a comparison with the scaling laws in the literature, which they found to be not very conclusive, but they did not discuss the effect of θ_max on the observations. It is therefore difficult to use these laws to build our simulations.

Appendix C Computation of the observables

We detail here how the observables were computed at each time step to produce the time series. The sum in all formulae is made on all structures of a given type at the corresponding time step, and all observables below are functions of time. We use the following notations: A_j is the size of the structure j (in ppm of the hemisphere); θ_j and ϕ_j are their latitude (between −90° and 90°) and longitude (between 0° and 360°), respectively; μ_j is their position on the disk (cosine between the local surface and the line of sight, it takes a value of 1 at disk center and zero at the limb, and depends on θ_j, ϕ_j, and inclination).

Other variables are as follows: i is the star inclination; Pcb is the center-to-limb darkening function at a given temperature, as in Borgniet et al. (2015), from Claret & Hauschildt (2003), using a log(g) of 4.5 and solar metallicity; C_pl is the relative contrast of the plages and is a function of μ; subscripts “phot”, “pl”, and “sp” are for photosphere, plages (including network), and spots, respectively; and T_eff is the photospheric temperature.

C.1 Filling factors of the spots and plages

The filling factor of either spots or plages is defined as $$\text{ff}=\Sigma A_j\times {\mu }_j\times 2$$(C.1)

and is in ppm of the stellar disk.

C.2 Photometry of the spots and plages

The plage contribution to the photometry is defined as $$\begin{array}{lll}{I}_{\text{pl}}\hfill & =\hfill & \Sigma {\text{A}}_j\times {\mu }_j\times C_{\text{pl}}({\mu }_j,{\text{A}}_j,T_{\text{eff}})\times 2\hfill \\ \hfill & \hfill & \times f_{\text{phot}}(T_{\text{eff}})\times \text{Pcb}({\mu }_j,T_{\text{eff}})/{\Phi }_{\text{tot}}(T_{\text{eff}}),\hfill \end{array}$$(C.2)

where C_pl is a relative contrast DeltaI/I = (Ipl − Iphot)/Iphot (local (= Iplage/Iphot − 1)) and Φ_tot is equal to f_phot × Pcb then integrated over the disk, where f is the Planck function for the indicated temperature.

The spot contribution to the photometry is defined as $$I_{\text{sp}}=\Sigma {\text{A}}_j\times {\mu }_j\times C_{\text{sp}}({\mu }_j,T_{\text{sp}},T_{\text{eff}})\times 2,$$(C.3)

where $$\begin{array}{lll}{C}_{\text{sp}}\hfill & =\hfill & (f_{\text{sp}}(T_{\text{sp}})\times \text{Pcb}({\mu }_j,T_{\text{sp}})-f_{\text{ph}}(T_{\text{eff}})\hfill \\ \hfill & \hfill & \times \text{Pcb}({\mu }_j,T_{\text{eff}}))/{\Phi }_{\text{tot}}(T_{\text{eff}}).\hfill \end{array}$$(C.4)

f and Φ_tot are defined as above. The photometric contributions are in ppm of the quiet-star brightness. The sum of the two provides the total variability in ppm.

C.3 RV of spots, plages, and convection inhibition

The plage contribution to the RV is defined as $$\begin{array}{lll}{\text{RV}}_{\text{pl}}\hfill & =\hfill & \Sigma A_{\text{j}}\times {\mu }_j\times 2\times {10}^{-6}\times \Omega ({\theta }_j)\times \sin ({\varphi }_j)\times \sin (i)\hfill \\ \hfill & \hfill & \times C_{\text{pl}}({\mu }_j,A_{\text{j}},T_{\text{eff}})\times \text{Pcb}({\mu }_j,T_{\text{eff}})\times {\text{f}}_{\text{ph}}(T_{\text{eff}})/{\Phi }_{\text{tot}}(T_{\text{eff}}),\hfill \end{array}$$(C.5)

where Ω is the rotation rate converted in m s⁻¹.

The inhibition of the convection contribution to the RV is defined as $$\begin{array}{lll}{\text{RV}}_{\text{conv}}\hfill & =\hfill & \Sigma {\text{A}}_j\times {\mu }_j\times 2\times {10}^{-6}\times \Delta \text{V}({\mu }_j)\times \text{Pcb}({\mu }_j,T_{\text{eff}})\hfill \\ \hfill & \hfill & \times (1+C_{\text{pl}}({\mu }_j,A_{\text{j}},T_{\text{eff}}))\times f_{\text{ph}}(T_{\text{eff}})/{\Phi }_{\text{tot}}(T_{\text{eff}}),\hfill \end{array}$$(C.6)

where ΔV is the attenuation of the convective blueshift (vertical contribution, see Sect. 3.2 for a discussion).

The spot contribution to the RV is defined as $$\begin{array}{lll}{\text{RV}}_{\text{sp}}\hfill & =\hfill & \Sigma A_{\text{j}}\times {\mu }_j\times 2\times {10}^{-6}\times \Omega ({\theta }_j)\times \sin ({\varphi }_j)\times \sin (i)\hfill \\ \hfill & \hfill & \times C_{\text{sp}}({\mu }_j,T_{\text{sp}},T_{\text{eff}}),\hfill \end{array}$$(C.7)

where C_sp is defined as above. All RV are in m s⁻¹. The Zeemaneffect (Reiners et al. 2013) is not taken into account. The sum of the three components then provides the total RV variation in m s⁻¹.

C.4 Astrometric displacements

The astrometric contribution in the x direction (x taken along the equator) is the same for plages and spots: $$\Delta x=\Sigma I_j\times x_j\times R(R_{\odot })\times 180\times 3600\times /D_{\text{star}}/\pi \times N,$$(C.8)

where I_j is the individual contribution of a structure to the photometry (see Eqs (C.2) and (C.3)). It is normalized to the quiet star. The formula is similar in y (along the rotation axis): x_j and y_j are positions in this referential system, relative to the stellar radius. D_star is arbitrarily chosen to be 10 pc (for a star at a different distance, the corresponding factor must be applied). N is a normalization factor equal to 1/(1+I_pl+I_sp) to normalize with the actual flux of the star. Astrometric displacements are in μas.

Appendix D Chromospheric emission model: practical recipe

D.1 Model for a G2 star

We followed the approach described in Meunier (2018), with some simplifications because magnetic field maps are not available. The S-indexmodel consists of three components and is described as $$S(t)=S_{\text{basal}}+[S_{\text{act}}(t)+S_{\text{qs}}(t)]\times f/N,$$(D.1)

where N is equal here to 10⁶ (areas are in ppm) and the three contributions are therefore the following:

• The basal flux corresponds to stars with no activity (see discussion in Meunier 2018), and must not be confused with the lowerlimit in Log $R_{\text{HK}}^{\prime }$ discussed in Sect. 2.3.1. For B–V lower than 0.94, we used the basal S-index of 0.144 derived by Mittag et al. (2013). For B–V above 0.94, the true basal flux is significantly below the lower limit observed in Log $R_{\text{HK}}^{\prime }$: we used the basal flux built from models and observations made by Schröder et al. (2012). The corresponding basal flux used in this paper is shown as the lower solid line in Fig. 3.

• The active component is due to the plages and network features that are simulated in Sect. 2. We attributed a magnetic field to each size. The sizes used in this work correspond to a threshold of 100 G on MDI/SOHO Scherrer et al. (1995) magnetograms (see Sect. 2.1), as defined in Meunier et al. (2010a). These sizes were then converted into sizes corresponding to a 40 G threshold as in Meunier (2018) before we used the scaling law of the magnetic field versus size from MDI/SOHO. Then we computed the contribution of each plage as a multiplying factor times this magnetic field to a certain exponent. The parameters were kept to the solar values (see Appendix D.3 for a discussion) and depend on the size of the region, as detailed in Meunier (2018).

• The quiet-star component (QS) is due to the weak magnetic field everywhere else on the surface. This contribution is critical and necessary (see Meunier 2018, for a discussion): if we had kept the solar relationship between filling factor ff and S-index as obtained by Shapiro et al. (2014) for example, it would be impossible to model the S-index for stars with an activity level much below the minimum solar cycle (since ff = 0 then). We assumed that the magnetic flux in the quiet star is related to the average activity level of the star: if it is more active, it provides more flux to the quiet star and therefore a higher chromospheric emission. The procedure is detailed below.

We tested a description of the QS component similar to the one used in our solar model (Meunier 2018), with the form f′ ×B${}_{\text{QS}}^{\beta }$, where B_QS is the average magnetic field in the quiet Sun. It is difficult to apply directly, however, if f′ is kept to the solar value, this contribution remains strong even for very low values of B_QS. If we adapt f′ to B_QS by allowing it to decrease sharply toward zero at low activity levels, it leads to very strong variations of the S-index with time due to the quiet star, which does not seems realistic. A QS contribution that is too large for quiet stars prevents us from modeling stars with a large cycle amplitude because the expected level at cycle minimum would then fall below the basal+QS level. Instead, we chose an empirical approach, given the lack of constraints, which allowed us to model stars with a large amplitude as observed. For a given cycle amplitude in our (B–V, Log $R_{\text{HK}}^{\prime }$) grid, given the average Log $R_{\text{HK}}^{\prime }$ and the minimum Log $R_{\text{HK}}^{\prime }$ corresponding to the amplitude, it is possible to compute the maximum QS contribution that we should allow to model such variability. Then we chose a QS level slighly below the minimum value: this imposed a small number of spots at cycle minimum, typically about 2–4. A more sophisticated model would require numerical simulations of magnetic flux tubes from active regions down to very small scales and is beyond the scope of this paper.

D.2 S-index for other stars

This S-index, built as in Appendix D.1, is valid for G2 stars (because it has been validated for the Sun) but not for other stars: the pertinent value to compare chromospheric emissions between stars is the Log $R_{\text{HK}}^{\prime }$, however. For a given simulation at B–V, we transformed the obtained S-index into a Log $R_{\text{HK}}^{\prime }$ using the solar value (0.65), then back into an S-index with the stellar Log $R_{\text{HK}}^{\prime }$. In practice, we mostly used the Log $R_{\text{HK}}^{\prime }$.

D.3 Possible intrisinc chromospheric variability versus B–V?

A dependence of the chromospheric emission (for similar plages) on the spectral type is not taken into account in the previous model. There have been very few studies of the chromospheric emission that would allow us to answer this question precisely. Numerical simulations of magnetic structures (Steiner et al. 2014; Beeck et al. 2015) suggest that the magnetic field remains similar in their magnetic regions (and if the magnetic flux increases, it is mostly due to the size, not the flux density). Therefore, we do not expect large variations. Cuntz et al. (1999) expected the rotation period to affect the chromospheric emission because a low rotation would lead to a smaller spreading of the flux tubes, which in turn would cause a slightly higher chromosperic emission, but it is difficult to directly apply to our parameters. Fawzy et al. (2002c), based on models of magnetic waves propagating in flux tubes to explain the chromospheric heating (Fawzy et al. 2002a,b), established a relation between chromospheric emission and magnetic coverage. They obtained an increasing emission for larger B–V. We can derive an average trend from their results, which would give a corrective factor equal to 1.12–0.18 ×(B–V). However, this is very uncertain because they specified that they might have uncertainties of up to a factor two on the flux emissions. Therefore, we chose not to include this trend in our model.

Appendix E Input parameters

In this appendix we summarize the parameters. Table E.1 shows the list of parameters and their values (or ranges). These were used to produce the list of spots, plages, and network features at each time step (size and position). Table E.2 shows the same summary for the parameters that allowed us to compute the observables. Finally, Fig. E.1 illustrates the range of values covered by the parameters, which are adapted to each spectral type as a function of B–V.

Fig. E.1 Summary of the input parameters vs. B–V, illustrating the range they cover for each spectral type, for the 11 421 simulations. All activity levels are superimposed. The color code corresponds to different average activity levels (well illustrated in the first panel). Temperatures correspond to Teff (triangles), solar spot contrast (stars), and the largest spot contrast (diamonds).

References

All Tables

All Figures

	Fig. 1 Variable parameters in our grid. Parameters with specific indications in orange depend on B–V and/or Log $R_{\text{HK}}^{\prime }$, and we also consider two or three laws, as indicated. A few parameters depend on specific conditions (plage contrast). The 2D grid in B–V and Log $R_{\text{HK}}^{\prime }$ is detailed in Fig. 3.
In the text

	Fig. 2 First panel: Teff vs. spectral type. Second panel: B–V vs. spectral type. Third panel: stellar radius (in R⊙) vs. spectral type. Fourth panel: stellar mass (in M⊙) vs. spectral type.
In the text

Fig. 3 Upper panel: average Log $R_{\text{HK}}^{\prime }$ versus B–V representing our 2D grid (orange stars). The lower solid line shows the basal flux used in this paper, from Mittag et al. (2013) and Schröder et al. (2012) and the upper solid line shows the upper bound for the activity level versus B–V, from Lockwood et al. (2007). The dashed line corresponds to the minimum Log $R_{\text{HK}}^{\prime }$ we consider for stars with B–V above 0.94. The vertical dotted lines show the range in B–V covered by our simulations, the vertical dashed line vizualises B–V = 0.94, and the red horizontal line approximately corresponds to the Vaughan-Preston gap (see text). Lower panel: same for the S-index versus B–V.

In the text

	Fig. 4 Chosen rotation periods vs. B–V for eight different Log $R_{\text{HK}}^{\prime }$ values between −5.1 and −4.75 (from top to bottom) from Mamajek & Hillenbrand (2008). The median law is shown in black, and extreme laws are shown in green and red.
In the text

	Fig. 5 Upper panel: coefficient p0 vs. Teff from the fit of Log(α) = p0+p1×Log(Prot) for stars with log(g) between 3.94 and 4.94. Computations are made from the data published by Reinhold & Gizon (2015). Middle panel: same for p1. Lower panel:number of stars in each Teff bin.
In the text

	Fig. 6 Pcyc (in years) vs.Prot used in our grid (upper bound as solid line, lower bound and intermediate range as dashed lines).
In the text

Fig. 7 Upper panel: half full amplitude of stellar cycles vs. B–V, derived from Lovis et al. (2011) after revision of the largest amplitudes (see text) for different types of stars: B–V < 0.7 (black stars), 0.7 < B–V < 0.9 (red squares), and B–V > 0.9 (green triangles). The black lines correspond to the lower and upper boundaries that were taken into account in building the grid. Upper panel: same vs. average Log $R_{\text{HK}}^{\prime }$. The solid line represents the lower limit that was taken into account in building the grid.

In the text

	Fig. 8 Wolf number vs. phase smoothed over the last solar cycle that was used as a reference.
In the text

Fig. 9 Convective blueshift derived for the basal Log $R_{\text{HK}}^{\prime }$ (dashed line), adapted from Meunier et al. (2017b), and scaled with the solar values derived in Meunier et al. (2017a). The solid line shows the local convective blueshift we used in our model, after correction for projection effects for a constant attenuation of the convective blueshift in plages of 0.38 (see text). The dotted line shows the same parameter when the trend of the attenuation factor vs. Teff is taken intoaccount.

In the text

	Fig. 10 Upper and lower bounds for the difference between spot and photosphere temperature vs. Teff from Berdyugina (2005). The lower bound is the solar value derived in Borgniet et al. (2015).
In the text

	Fig. 11 RV jitter due to the oscillation, granulation, and supergranulation vs. B–V for instantaneous values (black) and averaged over one hour (red).
In the text

Fig. 12 Panel A: Log ${R^{\prime }}_{\text{HK},\text{out}}$–Log ${R^{\prime }}_{\text{HK},\text{obj}}$ vs. Log ${R^{\prime }}_{\text{HK},\text{obj}}$ for all simulations. The color indicates stellar inclination from 0° (pole-on, yellow) to 90° (edge-on, blue). Panel B: Log ${R^{\prime }}_{\text{HK},\text{out}}$ averaged over all inclinations vs. Log ${R^{\prime }}_{\text{HK},\text{obj}}$. The solid red line indicates the y = x linear function. Panel C: ΔLog ${R^{\prime }}_{\text{HK},\text{out}}$–ΔLog ${R^{\prime }}_{\text{HK},\text{obj}}$ vs. Log ${R^{\prime }}_{\text{HK},\text{obj}}$, where Δ represents the amplitude of the cycle, for different inclinations (the color code is similar to panel A). Panel D: ΔLog ${R^{\prime }}_{\text{HK},\text{out}}$–ΔLog ${R^{\prime }}_{\text{HK},\text{obj}}$ vs. Log ${R^{\prime }}_{\text{HK},\text{obj}}$, where Δ representsthe amplitude of the cycle, after averaging the simulations made for inclination of 40° and 50°. Panel E: ΔLog ${R^{\prime }}_{\text{HK},\text{out}}$ vs. ΔLog ${R^{\prime }}_{\text{HK},\text{obj}}$, after averaging for inclination of 40° and 50°. The solid red line indicates the y = x linear function. Panel F: cycle period (in years) derived from a fit on smoothed Log $R_{\text{HK}}^{\prime }$ time series vs. the prescribed cycle period.

In the text

	Fig. 13 Firstline: Log $R_{\text{HK}}^{\prime }$ (left panel) and RV (right panel) vs. time for a moderately active G2 star and an inclination of 90°. Second line: same for an inclination of 0°. Third line: same for a quiet G2 star and an inclination of 90°. Fourth line: same for a moderately active K2 star and an inclination of 90°.
In the text

Fig. A.1 First panel: Comparison of the rotation period vs. B–V for Log $R_{\text{HK}}^{\prime }$ = −5.1 from various sources: the solid black line is from Mamajek & Hillenbrand (2008), the green dot-dashed line from Noyes et al. (1984a), and the dashed red line from Saar & Brandenburg (1999). Second panel: Same for Log $R_{\text{HK}}^{\prime }$ = −4.75. Third panel: Similar comparison vs. Log $R_{\text{HK}}^{\prime }$ for eight levels in B–V between 0.5 and 1.1. The vertical line approximately shows the lower limit for the validity of the laws (the shaded area indicates the zones where it is not valid).

In the text

	Fig. E.1 Summary of the input parameters vs. B–V, illustrating the range they cover for each spectral type, for the 11 421 simulations. All activity levels are superimposed. The color code corresponds to different average activity levels (well illustrated in the first panel). Temperatures correspond to Teff (triangles), solar spot contrast (stars), and the largest spot contrast (diamonds).
In the text