© ESO, 2011

1. Introduction

Several studies have recently addressed the perturbation of the RV induced by stellar activity. This jitter may severely hamper the detection of low-mass planets or the measurement of the mass of those found by the method of transits. Lagrange et al. (2010), Meunier et al. (2010a), and Meunier et al. (2010b) have considered the Sun as a star finding that the reduction of convective blueshifts in active regions, due to the quenching of turbulent motions by the photospheric magnetic fields, has a remarkable effect inducing variations with an amplitude up to 5–10 m s^-1. Dumusque et al. (2011) have extended those studies introducing a model for the impact of dark spots on the RV variations of solar-like stars based on the statistics of sunspot groups.

Several techniques have been proposed to reduce the impact of the RV perturbation induced by stellar activity. They can be classified into two broad groups. Those in the first group exploit the correlation between the RV perturbation and other properties of the cross-correlation function used to measure the RV itself, such as the mean line bisector shift and/or its FWHM (e.g., Melo et al. 2007; Boisse et al. 2009; Pont et al. 2011). The main limitation of such techniques is the lack of a tight correlation between the RV perturbation and those proxies for the average spectral line distortion, as demonstrated by, e.g., Boisse et al. (2011). The other group consists of the techniques that exploit the modulation of the activity-induced perturbation by the stellar rotation. Most of the power due to the activity-induced perturbation is indeed concentrated at the rotation frequency and its first two or three harmonics. Therefore, it is possible to filter it out leaving the signal due to the reflex orbital motion induced by a planetary companion, unless its orbital period falls close to such frequencies, and to fit simultaneously activity and planetary signals (see Boisse et al. 2011,for details). This approach has been recently applied, among others, to detect a second non-transiting planet around the active star CoRoT-7 and to refine the measurement of the mass of its transiting telluric planet CoRoT-7b (see Queloz et al. 2009; Boisse et al. 2011).

In addition to those methods, another approach to reduce the impact of stellar activity can be based on the models of the active region distribution derived from high-precision photometry such as that provided by the space-borne telescopes of MOST, CoRoT, or Kepler for their target stars. From their light curves, it is possible to derive the longitudes of the active regions and their area, tracing their evolution along successive rotations. This information can be used to model the RV perturbation induced by stellar activity. Lanza et al. (2010) have applied this approach to the light curve of CoRoT-7 to predict the power spectrum of the RV perturbation induced by stellar activity and derive the significance of the signals attributed to the telluric planets CoRoT-7b and CoRoT-7c. The main limitation of that study was the lack of simultaneity between the CoRoT light curve and the HARPS RV time series analysed to extract the reflex motions induced by the two planets.

Here we further explore this approach by means of simultaneous photometric and RV measurements of HD 189733 making use of the MOST light curve and the SOPHIE RV time series published by Boisse et al. (2009). They span about one month in the summer of 2007, corresponding to two full rotations of the star with a period of P_rot ~ 12 days (Henry & Winn 2008). HD 189733 is accompanied by a Jupiter-mass transiting planet with an orbital period of 2.218 days. With an effective temperature of  ~5000 Kand an average chromospheric index S_HK = 0.525, this K1 main-sequence star is remarkably active (Bouchy et al. 2005; Moutou et al. 2007; Boisse et al. 2009). Its radial velocity (hereafter RV) measurements are then perturbed with variations up to  ~20–30 m s^-1, as evident in the residuals of the RV best fit of the planetary orbit (Boisse et al. 2009).

Our spot modelling technique allows us to retrieve the location and the evolution of the photospheric active regions of HD 189733 finding also evidence of a remarkable stellar differential rotation, as recently found by Fares et al. (2010) by means of Doppler Imaging techniques. The information on the starspot longitude and evolution is used to synthesize the RV variations that can be compared with the simultaneous SOPHIE measurements to assess the advantages and limitations of the technique. We find that the modulation of the RV perturbation induced by stellar activity is reproduced remarkably well on timescales comparable with the stellar rotation period but it is not possible to model the variations with a timescale of 1–2 days, i.e., much shorter than the stellar rotation period. This is an unavoidable consequence of the fact that the rotational modulation of the optical flux is used to map starspots. On the other hand, the main advantage of this method is the remarkable reduction of the perturbations at frequencies comparable with the stellar rotation frequency and its low-order harmonics, that is an important improvement over previously proposed approaches.

The next section describes the observational dataset. Then, we define our model: 1) how active region locations are derived from lightcurve; 2) how we compute the distortion of the line profiles due to the brightness perturbations owing to dark spots and bright faculae and the reduction of the convective blueshifts in the active regions; and 3) how the distributions of active regions derived from the spot model fitted to the light curve are used to synthesize RV variations. The fourth section gives the parameters chosen for our modelling of HD 189733. The results, given in Sect. 5, are discussed in the final section.

2. Observations

2.1. MOST photometry

In 2007, HD 189733 was monitored by the MOST (Microvariability and Oscillations of STars) microsatellite for 30.4994 days covering  ~2.5 rotation cycles of the star from HJD 2 454 298.55874 to 2 454 329.05818. MOST was launched in 2003 and was designed to observe nearby, relative bright stars through a single broadband (350–700 nm) filter using a 15/17.3-cm Rumak-Matsutov telescope (Walker et al. 2003; Matthews et al. 2004). The photometry was collected in the so-called MOST’s Direct Imaging mode, where a defocussed stellar image was recorded on a subraster of a single CCD. The data were reduced as described in Rowe et al. (2006). The stellar fluxes were extracted by aperture photometry and the raw instrumental light curve was decorrelated against the sky background and also against the location of the centroid of the point-spread function on the CCD. The satellite had an 820 km altitude Sun-synchronous polar orbit of 101.4 min and the data were binned to this period. The transits of the known planet were removed from the MOST lightcurve. This led to a time series consisting of 414 points. As seen in the top panel of Fig. 3, HD 189733 exhibited in 2007 a  ~1% variation in flux, which is much less than the  ~3% variation observed in 2006 by the same telescope (Miller-Ricci et al. 2008; Croll et al. 2007). Moreover, a decreasing long-term trend was seen that was not exhibited by the comparison stars, suggesting that the drift was real.

We assume the maximum observed flux at HJD 2 454 300.31968 as the unit of measure for the flux and normalize the light curve to that value. The mean error of the binned points is 1.15 × 10^-4 in relative flux units. However, since MOST photometry is non-differential and can experience drifts on daily or weekly timescales, such a value should be considered as a lower limit to the photometric precision and does not take into account systematic errors. For example, Croll et al. (2007) pointed out that the stray light background can be modulated at a period of 1 day and its first harmonic (due to the Sun-synchronous orbit). The flux level corresponding to the star without spots is unknown, therefore we assume the maximum observed flux as the reference level for the spot modelling, i.e., we assume that it corresponds to the flux of the unspotted star in the MOST passband. We can also mention that the M-dwarf companion of HD 189733 (Bakos et al. 2006) does not contaminate appreciably the photometric signal.

2.2. SOPHIE radial velocity measurements

SOPHIE is a high-resolution spectrograph mounted on the 1.93-m telescope at the Observatoire de Haute-Provence, France (Perruchot et al. 2008; Bouchy et al. 2009b). The SOPHIE observations of HD 189733 were conducted simultaneously with the MOST photometry. Usually two exposures per night were gathered during more than one month (from 12 July to 23 August 2007) to sample both the short orbital period (~2.2 days) and the longer rotation period of the star (~12 days). The observations were conducted in the high-resolution mode (λ/Δλ ~ 75   000) and with an average signal-to-noise ratio (SNR) per pixel of about 80 at λ = 550 nm. The journal of the observations and the instrument setup were reported in Sect. 2 of Boisse et al. (2009).

For the present work, the original spectra were reduced again with a new pipeline that, among others, corrects for the effects of the charge transfer inefficiency of the CCDs (Bouchy et al. 2009a) and removes part of the spectrum that may contain telluric lines before applying the correlation mask. This gives slightly different weights to the red and blue parts of the spectrum for the determination of the RV through the standard cross-correlation method, resulting in a RV offset of +21.6 m s^-1 with respect to the RV values reported in Boisse et al. (2009). The recently reported seeing effect (Boisse et al. 2010) was not well constrained and observed in these data, therefore it is not corrected.An average error of 3.9 m s^-1 for the RV measurements is then a good estimate of the accuracy of our values accounting for the photon-noise error, the uncertainty in the wavelength calibration, and SOPHIE systematics. The effects of the systematics is taken to be 3 m s^-1 for all the measurements, derived from the standard deviation of the RV of the most stable stars observed with SOPHIE, some of which have been observed with HARPS to have a stability better than 1 m s^-1. These stars have more than 10 measurements and have been observed for more than six months in order to take into account the long-term systematics. We note that new fiber scramblers have just been installed on SOPHIE (June 2011), and preliminary tests showed a significant improvement in the RV accuracy. To compute the RV perturbations due to stellar activity, the orbital solution has been subtracted. The adopted orbital elements are those given by Boisse et al. (2009).

3. Light curve and RV modelling

3.1. Spot modelling of wide-band light curves

The reconstruction of the surface brightness distribution from the rotational modulation of the stellar flux is an ill-posed problem, because the variation of the flux vs. rotational phase contains information only on the distribution of the brightness inhomogeneities vs. longitude. The integration over the stellar disc effectively cancels any latitudinal information, particularly when the inclination of the rotation axis along the line of sight is close to 90°, as in the present case (see Sect. 4 and Lanza et al. 2009a). Therefore, we need to include a priori information in the light curve inversion process to obtain a unique and stable map. This is done by computing a maximum entropy (hereinafter ME) map, which has been proven to best reproduce active region distribution and area variations in the case of the Sun (cf. Lanza et al. 2007). For a different modelling approach, based on discrete starspots, see, e.g., Mosser et al. (2009).

In the present model, the stellar surface is subdivided into elements, i.e., in 200 squares of side 18°, with each element containing unperturbed photosphere, dark spots, and facular areas. The fraction of the kth element covered by dark spots is indicated by its filling factor f_k, the fractional area of the faculae is Qf_k, and the fractional area of the unperturbed photosphere is 1 − (Q + 1)f_k. The contribution to the stellar flux coming from the kth surface element at time t_j, where j = 1,...,N, is an index numbering the N points along the light curve, is given by: (1)where I₀ is the specific intensity in the continuum of the unperturbed photosphere at the isophotal wavelength of the observations, c_s and c_f are the spot and facular contrasts, respectively (cf. Lanza et al. 2004), A_k is the area of the kth surface element, (2)is its visibility, and (3)is the cosine of the angle ψ_kj between the normal to the surface element and the direction of the observer, with i being the inclination of the stellar rotation axis to the line of sight, θ_k the colatitude and ℓ_k the longitude of the kth surface element, Ω the angular velocity of rotation of the star (Ω ≡ 2π/P_rot), and t₀ the initial time. The specific intensity in the continuum varies according to a quadratic limb-darkening law, as adopted by Lanza et al. (2003) for the case of the Sun, viz. I₀ ∝ a_p + b_pμ + c_pμ². The stellar flux computed at time t_j is then: F(t_j) =  ∑ _kΔF_kj. To warrant a relative precision of the order of 10^-5 in the computation of the flux F, each surface element is further subdivided into elements of 1° × 1° and their contributions, calculated according to Eq. (1), are summed up at each given time to compute the contribution of the 18° × 18° surface element to which they belong.

We fit the light curve by changing the value of the spot filling factor f over the surface of the star, while Q is held constant. Even fixing the rotation period, the inclination, and the spot and facular contrasts (see Lanza et al. 2007,for details), the model has 200 free parameters and suffers from non-uniqueness and instability. To find a unique and stable spot map, we apply ME regularization, as described in Lanza et al. (2007), by minimizing a functional Z, which is a linear combination of the χ² and the entropy functional S; i.e., (4)where f is the vector of the filling factors of the surface elements, η > 0 a Lagrangian multiplier determining the trade-off between light curve fitting and regularization, andthe expression of S is: (5)where w_k is the relative area of the kth surface element (total surface area of the star  = 1) and m the default spot covering factor that fixes the limiting values for f_k as: m < f_k < (1 − m). In our modelling we adopt m = 10^-6 (cf. Lanza et al. 1998).The entropy functional S is constructed in such a way that it attains its maximum value when the star is virtually immaculate (f_k = m for every k). Therefore, by increasing the Lagrangian multiplier η, we increase the weight of S in the model and the area of the spots is progressively reduced. This gives rise to systematically negative residuals between the observations and the best-fit model when η > 0.

The optimal value of η depends on the information content of the light curve, which in turn depends on the ratio of the amplitude of its rotational modulation to the average standard deviation of its points. In the case of HD 189733, the amplitude of the rotational modulation is  ~0.015, while the nominal standard deviation of the points is  ~1.15 × 10^-4 in relative flux units (see Sect. 2.1), giving a signal-to-noise ratio of  ~130. However, since short-term instrumental fluctuations cannot be excluded in the case of MOST photometry (cf. Sect. 2.1), we prefer to consider the standard deviation of the residuals obtained with the unregularized best fit (i.e., with η = 0; see Sect. 5.1). In such a way, we find that the relative accuracy of MOST photometry is overestimated by at least a factor of  ~3, i.e. the mean standard deviation of the photometric points is  ~4 × 10^-4, leading to a signal-to-noise ratio of  ~40–50.

To fix the optimal value of the Lagrangian multiplier η, we compare the modulus of the mean of the residuals of the regularized best fit |μ_reg| with the standard error of the residuals themselves, i.e., , where σ₀ is the standard deviation of the residuals of the unregularized best fit and N is the number of points in each fitted subset of the light curve having a duration Δt_f (see below). We iterate until |μ_reg| ≃ βϵ₀, where β is a numerical parameter that will be fixed a posteriori according to the requisites of an acceptable quality for the fit and a regular evolution of the spot pattern from one spot map to the next (cf. the case of CoRoT-4 in Lanza et al. 2009b).

In the case of the Sun, by assuming a fixed distribution of the filling factor, it is possible to obtain a good fit of the irradiance changes only for a limited time interval Δt_f, not exceeding 14 days, which is the lifetime of the largest sunspot groups dominating the irradiance variation (Lanza et al. 2003). In the case of other active stars, the value of Δt_f must be determined from the observations themselves, looking for the maximum data extension that allows us a good fit with the applied model (see Sect. 4).

The optimal values of the spot and facular contrasts and of the facular-to-spotted area ratio Q in stellar active regions are unknown a priori. In our model the facular contrast c_f and the parameter Q enter as the product c_fQ, so we can fix c_f and vary Q, estimating its best value by minimizing the χ² of the model, as shown in Sect. 4. Since the number of free parameters of the ME model is large, we make use of the model of Lanza et al. (2003)to fix the value of Q. It fitsthe light curve by assuming only three active regions to model the rotational modulation of the flux plus a uniformly distributed background to account for the variations of the mean light level. This procedure is the same as adopted for CoRoT-2 and CoRoT-4 to fix the value of Q (cf. Lanza et al. 2009a,b).

We shall assume an inclination of the rotation axis of HD 189733 of i ≃ 85° (see Sect. 4). Since the information on spot latitudes that can be extracted from the rotational modulation of the flux for such a high inclination is negligible, the ME regularization virtually puts all the spots at the sub-observer latitude (i.e., 90° − i ≈ 5°) to minimize their area and maximize the entropy. Therefore, we are limited to mapping only the distribution of the active regions vs. longitude, which can be done with a resolution of at least  ~50° (cf. Lanza et al. 2007, 2009a). Our ignorance of the true facular contribution to the light modulation may lead to systematic errors in the active region longitudes derived by our model, as discussed by Lanza et al. (2007) in the case of the Sun.

3.2. Modelling the line distortion induced by activity

Surface magnetic fields in late-type stars produce brightness and convection inhomogeneities that distort their spectral line profiles leading to apparent RV variations (cf., e.g., Saar & Donahue 1997; Saar et al. 1998). To compute the apparent RV variations induced by stellar active regions, we adopt a simple model for each line profile considering the residual profile R(λ) at a wavelength λ along the line, i.e.: R(λ) ≡  [1 − I(λ)/I_c], where I(λ) is the specific intensity at wavelength λ and I_c the intensity in the continuum adjacent to the line. The local residual profile is assumed to be a Gaussian with thermal and microturbulent width Δλ_D (cf., e.g., Gray 1988,Ch. 1), i.e.: (6)where λ₀ is the local central wavelength. The different effective temperatures in dark spots or bright faculae modify the depth of the local line profile. To take into account this effect, we multiply the unperturbed profile given by Eq. (6) by (1 − C_df_k), where C_d is a coefficient that depends on the specific line considered and the physical properties of the stellar atmosphere, and f_k is the spot filling factor in the specified surface element. This factor accounts for the temperature reduction in dark spots, but does not consider the increased effective temperature in faculae because the facular contrast is assumed to be solar-like (cf. Sect. 4), i.e., corresponding to a modest temperature increase of  ≈ 100 K in the layers where the optical spectrum is produced (cf., e.g., Unruh et al. 1999). Since the RV measurement is based on a cross-correlation function, C_d must be averaged over a great number of spectral lines having in general quite different temperature responses. For this reason, we consider C_d as a free parameter in our model and adjust it to fit the observed RV variations induced by stellar activity (see below for the effect of the variation of C_d on the RV perturbation induced by a starspot). In addition to the variation of the local line profile, the local specific intensity along a spectral line I(λ) depends on the variation of the continuum intensity I_c owing to limb-darkening and the effects of dark spots and bright faculae, as specified by Eq. (1).

To include the effects of surface magnetic fields on convective motions, we consider the decrease of macroturbulence velocity and the reduction of convective blueshifts in active regions. Specifically, we assume a local radial-tangential macroturbulence, as introduced by Gray (1988), with a distribution function Θ for the radial velocity v of the form: (7)where ζ_RT is the macroturbulence and μ is given by Eq. (3). The profile emerging from each surface element is obtained by convolving the local profile in Eq. (6) with the Doppler shift distribution as generated by the macroturbulence function in Eq. (7). The parameter ζ_RT is assumed to be reduced in spotted and facular areas, according to their total filling factor, i.e., ζ_RT = ζ₀ [1 − (Q + 1)f], where ζ₀ is the unperturbed macroturbulence. This is a convenient parameterization of the quenching effect of surface magnetic fields on convective turbulence, at least in the framework of our simplified model, and is supported by the observations of the reduction of turbulent velocity fields in the plage regions of the Sun (cf., e.g., Title et al. 1992).

Convective blueshifts arise because in stellar photospheres most of the flux comes from the extended updrafts at the centre of convective granules. At the centre of the stellar disc, the vertical component of the convective velocity produces a maximum blueshift, while the effect vanishes at the limb where the projected velocity is zero. The cores of the deepest spectral lines form in the upper layers of the photosphere where the vertical convective velocity is low, while the cores of shallow lines form in deeper layers with higher vertical velocity. Therefore, the cores of shallow lines are blueshifted with respect to the cores of the deepest lines. Gray (2009) shows this effect by plotting the endpoints of line bisectors of shallow and deep lines on the same velocity scale. He shows that the amplitude of the relative blueshifts scales with the spectral type and the luminosity class of the star. For the G8V star τ Ceti, which has a spectral type not too much different from the K1V of HD 189733, the convective blueshifts are approximately similar to those of the Sun, so we adopt solar values in our simulations. In active regions, vertical convective motions are quenched, so we observe an apparent redshift of the spectral lines in spotted and facular areas in comparison with the unperturbed photosphere. Meunier et al. (2010b) quantify this effect in the case of the Sun, and we adopt their results, considering an apparent redshift ΔV_f = 200 m s^-1 in faculae and ΔV_s = 300 m s^-1 in cool spots.

In principle, the integrated effect of convective redshifts can be measured in a star by comparing RV measurements obtained with two different line masks, one including the shallow lines and the other the deep lines (cf. Meunier et al. 2010b). For HD 189733, which lacks such measurements, we apply the results of Gray (2009) and adopt solar-like values as the best approximation.

Considering solar convection as a template, intense downdrafts are localized in the dark lanes at the boundaries of the upwelling granules, but they contribute a significantly smaller flux because of their lower brightness and smaller area. While a consideration of those downdrafts is needed to simulate the shapes of line bisectors, it is beyond the scope of our simplified approach that assumes that the whole profile of our template line forms at the same depth within the photosphere. Therefore, we restrict our model to the mean apparent RV variations by neglecting the associated variations of the bisector shape and do not include the effect of convective downdrafts, as well as those of other surface flows typical of solar active regions, such as the Evershed effect in sunspots (cf., e.g., Meunier et al. 2010b).

The local central wavelength λ₀ of the kth surface element at time t_j is given by λ_0kj = λ_R(1 + v_kj/c), where v_kj is its radial velocity, (8)with c the speed of light, vsini the stellar projected rotational velocity, λ_R the rest wavelength of the spectral line, and ΔV_cs the apparent convective redshift arising from the reduction of convective blueshifts in spots and faculae, which is parameterized as (9)where f_k is the spot filling factor of the element, Q the facular-to-spotted area ratio, and μ_kj is given by Eq. (3).

We integrate the line specific intensity at a given wavelength over the disc of the star using a subdivision into 1° × 1° surface elements to obtain the flux along the line profile F(λ,t_j) at a given time t_j. To find the apparent stellar RV, we can fit a Gaussian to the line profile F(λ,t_j), or we can determine the centroid of the profile as (10)where ℛ ≡  [1 − F(λ,t_j)/F_c(t_j)] is the residual profile of the spectral line computed from the ratio of the flux in the line F(λ,t_j) to the flux F_c(t_j) in the adjacent continuum at any given time t_j. For HD 189733, the two methods lead to almost equal values of the RV perturbation with differences never exceeding  ~ 3 percent.

Fig. 1 Upper panel: the RV variation produced by a single equatorial spot on HD 189733 vs. the rotation phase for different values of the parameter Cd according to the different linestyles: Cd = 0.0 (no temperature effect on the line profile): solid line; Cd = 0.2: dotted line; and Cd = 0.4: dashed line. Lower panel: the corresponding modulation of the optical flux vs. the rotation phase. Note that the three plots overlap because the modulation is independent of Cd (see the text).

A single line profile computed with the above model can be regarded as a cross-correlation function (hereafter CCF) obtained by cross-correlating the whole stellar spectrum with a line mask consisting of Dirac delta functions giving the rest wavelength and depth of each individual line (e.g., Queloz et al. 2001). Therefore, to derive the RV from a single synthetic line profile is equivalent to measuring the RV from a CCF, either by fitting it with a Gaussian or by computing its centroid wavelength. A better approach would be to simulate the whole stellar optical spectrum or an extended section of it to account for the wavelength dependence of the spot and facular contrasts, as well as of the convective inhomogeneities (cf., e.g., Desort et al. 2007; Lagrange et al. 2010; Meunier et al. 2010b). Again, in view of our limited knowledge of the distribution of active regions on the stellar surface, we believe that our simplified approach is adequate for estimating the RV perturbations induced by magnetic fields.

An illustration of the effects of dark spots, bright faculae, and macroturbulence perturbations on the stellar RV measurements is presented in Lanza et al. (2010) together with a comparison of the results of the present model with those of previous approaches. Here we only show in Fig. 1 (upper panel) the effect of the variation of the parameter C_d on the RV curve of a single spot that was not considered by Lanza et al. (2010). The adopted spot is located on the equator of the star and has a square shape with a side of 18°, a covering factor f = 0.99, and a contrast c_s = 0.482 without any facular contribution. The stellar and line profile parameters are assumed to be the same as adopted for the modelling of the RV variations of HD 189733 (see Sect. 4). The increase of the amplitude of the RV variation with C_d is a consequence of the higher bump associated with the spot on the line profile when the depth of the local profile is reduced by the factor (1 − C_df). On the other hand, the modulation of the optical flux associated with the transit of the spot across the stellar disc is independent of the value of C_d because it depends on the flux in the continuum outside the spectral line (cf. Fig. 1, bottom panel).

The amplitude of the RV variation produced by an active region depends on several parameters that are poorly known, i.e., the latitude of the active region, the spot and facular contrasts, and the macroturbulence parameter which is difficult to disentangle from rotational broadening in a slowly rotating star such as HD 189733 (Léger et al. 2009). Moreover, the spot and facular contrasts depend on the wavelength (Lanza et al. 2004), leading to a difference of  ≈ 10 percent in the RV variations as derived from different orders of an echelle spectrum (cf., e.g., Desort et al. 2007). The simultaneous presence of several active regions gives rise to a complex line profile distortion in the case of a slowly rotating star because the perturbation of each active region is not separated from the other by the rotational broadening. This implies an additional 10–15 percent uncertainty in the determination of the central wavelength of the profile by the Gaussian fit or the centroid method, even for an infinite signal-to-noise ratio and perfectly regular sampling, as in the case of simulated line profiles. In consideration of all these limitations, the absolute values of the RV variations computed with our model are uncertain by 10–20 percent in the case of complex distributions of active regions, such as those derived by our spot modelling technique as applied to MOST and CoRoT light curves.

3.3. RV variations from spot modelling

We can use the distribution of active regions as derived from our spot modelling to synthesize the corresponding RV variations according to the approach outlined in Sect. 3.2. Since our spot model assumes that active regions are stable for the time interval of each fitted time series Δt_f, the distribution of surface inhomogeneities can be used to synthesize the RV variations having a timescale of Δt_f or longer. To improve the time resolution, we shall consider models obtained from time intervals whose initial time is shifted by Δt_f/2 from the previous interval, and average the RV variations computed at the same time to account for the spot evolution in a simple way. Our ME modelling minimizes the area of the starspots and puts them close to the equator. Since the spot latitudes cannot be determined from photometry (cf. Sect. 3.1), we assume that all the spots are located close to the equator, which maximizes the computed RV variations (cf. Eq. (8)). Note that the frequency content of the simulated RV time series is not significantly modified by this assumption because the duration of the transit of a given spot across the stellar disc is pratically independent of its latitude given that the inclination of the stellar rotation axis is close to 90° and our reference frame is rigidly rotating with the mean stellar rotation period. Moreover, the occurrence times of the maximum and minimum of the RV variations associated with a given spot are mainly determined by its longitude and its projected area that our modelling approach consistently retrieves from the modulation of the optical flux. Therefore, we do not expect any significant modification of the power spectrum of the synthetic RV variations associated with the assumption of equatorial spots.

The active regions with lifetimes shorter than Δt_f/2 produce a photometric modulation that appears in the residuals of the best fit to the lightcurve. As we shall see in Sect. 5.1, most of the short-term variability occurs on timescales of 1–2 days, i.e., significantly shorter than the rotation period of HD 189733, so we can neglect, as a first approximation, the variation in the disc position of active regions due to stellar rotation and estimate their area as if they were located at the centre of the disc. Of course, the actual RV perturbation depends on the position of the active regions on the disc, but this information is lacking. Therefore, the approach we describe below cannot be used to synthesize the actual RV variation, but only to estimate its maximum amplitude.

First, we express the photometric residuals as the relative deviation ΔF/F₀ between the observed flux and its best fit measured in units of the reference unspotted flux F₀. Then we subtract its mean value μ_res ≡  ⟨ ΔF/F₀ ⟩  from ΔF/F₀ because the mean value corresponds to a uniformly distributed population of active regions that does not produce any RV variation. At each given time t, we adopt as a measure of the filling factor of the active regions producing the short-term RV variations. In doing so, we neglect limb-darkening effects and assume that those active regions consist of dark spots with a contrast c_s. Finally, we compute the amplitude of the RV perturbation due to such brightness inhomogeneities by means of Eq. (1) of Desort et al. (2007) obtaining (11)where the RV perturbation is measured in m s^-1 and the vsini of the star in km s^-1. Equation (11) gives an upper limit for the amplitude of the RV perturbation because it assumes that an active region spans the diameter of the stellar disc. On the other hand, the effects of convective redshifts in short-lived active regions is estimated as (12)By adding both short-term perturbations, we compute the amplitude of the total variation induced by the active regions evolving on timescales shorter than the time resolution of our model, and add it to the standard deviation of the RV variations computed with the slowly evolving active regions to compute the total uncertainty of the synthesized RV values.

4. Model parameters

The basic stellar parameters are taken from Bouchy et al. (2005) and Triaud et al. (2009) and are listed in Table 1. The limb-darkening parameters have been derived from Kurucz (2000) model atmospheres for T_eff = 5050 K, log g = 4.53 (cm s^-2) and solar abundances, by adopting the MOST transmission profile in Walker et al. (2003). We compared our limb-darkening profile with that derived by Miller-Ricci et al. (2008), who adopted a different limb-darkening law and used MOST data acquired in 2006, and found a maximum relative difference of 3 percent.

The rotation period adopted for our spot modelling has been derived from the optical light curve (P_rot = 11.953 ± 0.009 days) by Henry & Winn (2008). The polar flattening of the star due to the centrifugal potential is computed in the Roche approximation with a rotation period of 11.95 days. The relative difference between the equatorial and the polar radii is ϵ = 2.58 × 10^-5 which induces a completely negligible relative flux variation of  ≈ 10^-6 for a spot coverage of  ~2 percent, as a consequence of the gravity darkening of the equatorial regions of the star.

The inclination of the stellar rotation axis is constrained by the observations of the Rossiter-McLaughlin effect by Triaud et al. (2009) who found that the angle β_RM between the projected directions on the plane of the sky of the orbital angular momentum and the stellar spin is . Although the possibility that the stellar spin be not aligned with the orbital angular momentum cannot be excluded with certainty because only the sky-projected angle between the two vectors has been measured, we regard this occurrence as highly improbable and assume that the rotation axis is inclined with respect to the line of sight of the same angle as the orbit normal.

The maximum time interval that our model can accurately fit with a fixed distribution of active regions Δt_f has been determined by dividing the total interval, T = 30.4994 days, into N_f equal segments, i.e., Δt_f = T/N_f, and looking for the minimum value of N_f that allows us a good fit of the light curve, as measured by the χ² statistics. We found that for N_f < 4 the quality of the best fit degrades significantly with respect to higher values, owing to a substantial evolution of the pattern of surface brightness inhomogeneities. We adopt N_f = 4, and therefore Δt_f = 7.625 days is the maximum time interval to be fitted with a fixed distribution of surface active regions in order to estimate the best value of the parameter Q (see below).

To compute the spot contrast, we adopt the same mean temperature difference as derived for sunspot groups from their bolometric contrast, i.e. 560 K (Chapman et al. 1994). In other words, we assume a spot effective temperature of 4490 K, yielding a contrast c_s = 0.482 in the MOST passband (cf. Lanza et al. 2007). A different spot contrast changes the absolute spot coverages, but does not significantly affect their longitudes and their time evolution, as discussed in detail by Lanza et al. (2009a). Therefore, adopting a spot temperature deficit of, e.g.,  ~ 1000 K, as discussed in the analysis of HD 189733 HST data by Pont et al. (2007), will decrease the filling factor in proportion to the lower contrast, but does not change our results in any significant way.

In the models computed with a facular contribution, the facular contrast is assumed to be solar-like with c_f = 0.115 (Lanza et al. 2004). The best value of the area ratio Q between the faculae and the spots in each active region has been estimated by means of the three-spot model by Lanza et al. (2003,cf. Sect. 3.1. In Fig. 2, we plot the ratio of the total χ² of the composite best fit of the entire time series to its minimum value , versus Q, and indicate the 95 percent confidence level as derived from the F-statistics (e.g., Lampton et al. 1976). The choice of Δt_f = 7.625 days allows us to fit the rotational modulation of the active regions for the longest time interval during which they remain stable, modelling both the flux increase due to the facular component when an active region is close to the limb as well as the flux decrease due to the dark spots when the same region transits across the central meridian of the disc. In such a way, a measure of the relative facular and spot contributions can be obtained leading to an estimate of Q. Figure 2 shows that the best value of Q is Q = 0, with an acceptable range extending from  ~0 to  ~5. Therefore, we adopt Q = 0 for all our modelling in Sect. 5 unless otherwise indicated. We shall further comment on the value of Q in Sect. 6.

To compute the RV variations induced by surface inhomogeneities, we assume a line rest wavelength of 600 nm and a local thermal plus microturbulence broadening Δλ_D = 2.3 km s^-1. The vsini = 3.316 ± 0.07 km s^-1 is derived from the observations of the Rossiter-McLaughlin effect (Triaud et al. 2009). Although such a measurement can be affected by a systematic error, we prefer to use this value, that is larger than the 2.97 ± 0.22 km s^-1 adopted by Winn et al. (2006), to maximize the synthesized RV variations. Adopting vsini = 2.97 km s^-1 would imply an amplitude of the RV variations lower by  ≈ 10 percent for the same value of the parameter C_d. On the other hand, we could still reproduce the observed variations by changing the average value of C_d (see Sect. 5.3), so that the choice of the vsini value is not very critical for our modelling. The radial-tangential macroturbulence velocity is assumed ζ₀ = 1.2 km s^-1 after Gray (1988) because it is difficult to disentangle from the rotational broadening owing to the slow rotation of the star.

Fig. 2 The ratio of the χ2 of the composite best fit of the entire time series of HD 189733 to its minimum value vs. the parameter Q, i.e., the ratio of the area of the faculae to that of the cool spots in active regions. The horizontal dashed line indicates the 95 percent confidence level for , determining the interval of acceptable Q values.

Fig. 3 Upper panel: the ME-regularized composite best fit to the out-of-transit light curve of HD 189733 obtained for Q = 0.0. The flux is in relative units, i.e., measured with respect to the maximum observed flux along the light curve. Lower panel: the residuals from the composite best fit versus time.

5. Results

5.1. Light curve models

We applied the model of Sect. 3.1 to the out-of-transit light curve of HD 189733, considering time intervals Δt_f = 7.625 days. The initial epoch of each fitted interval is shifted by Δt_f/2 = 3.81 days, to sample the evolution of the starspots in a more continuous way. The best fit without regularization (η = 0) has a mean μ_res = 2.15 × 10^-6 and a standard deviation of the residuals σ₀ = 3.99 × 10^-4 in relative flux units. The Lagrangian multiplier η is iteratively adjusted until the mean of the residuals , where N = 103 is the mean number of points in each fitted light curve interval Δt_f; the standard deviation of the residuals of the regularized best fit is σ = 4.86 × 10^-4. In other words, we fix β ≃ 2 in the regularization procedure which gives the best compromise between the accuracy of the fit and the smoothness of the spot distribution and evolution.

The composite best fit to the entire light curve is shown in the upper panel of Fig. 3 while the residuals are plotted in the lower panel. The computed flux values at the same time from successive spot models are averaged. In spite of that, at the matching points between successive best fits we observe some discontinuities in the first derivative of the flux variation vs. time. This may in principle affect the rate of change of the RV perturbation induced by the spots vs. the rotation phase. However, we do not worry about correcting for this effect because the discontinuities always occur at epochs when no RV observation has been obtained. Therefore, the synthesizes RV values to be compared with the SOPHIE RV time series are not appreciably affected by those discontinuities.

The residuals show oscillations with a typical timescale of  ~1–2 days that can be related to the rise and decay of active regions that cover  ≈ 0.2–0.3 percent of the stellar disc, i.e., comparable with the largest sunspot groups. The projected area of those active regions is estimated from the amplitude of the flux residuals and the adopted spot contrast, while their lifetimes are estimated by the duration of the residual fluctuations. These active regions cannot be modelled by our approach because they do not produce a significant rotational flux modulation during the  ~12 days of the stellar rotation period as they move across the disc by only  ≈ 60° in longitude. Their estimated size and lifetime are not significantly degenerate with each other because our model has 200 degrees of freedom to fit any minute flux variation with a timescale of several days by adjusting the covering factors f_k of its 18° × 18° surface elements. In other words, by subtracting our model we filter out all the variations with timescales longer than 4–5 days leaving only the effects of the active regions evolving on shorter timescales.

By decreasing the degree of regularization, i.e., the value of η, we can marginally improve the best fit, but at the cost of introducing several small active regions that wax and wane from one Δt_f time interval to the next and are badly constrained by the rotational modulation. Nevertheless, the oscillations of the residuals do not disappear completely even for η = 0, indicating that HD 189733 has a population of short-lived active regions with typical lifetimes of 1–2 days. Spots can be mapped along the strip occulted by the planet, as shown by Pont et al. (2007), but this is beyond the scope of the present investigation.

5.2. Longitude distribution of active regions and stellar differential rotation

Fig. 4 The distributions of the spotted area vs. longitude at the labelled times (HJD−2  454  000.0) for Q = 0.0. The plots have been vertically shifted to show the migration of individual spots (relative maxima of the distributions) versus time. The vertical short-dashed lines mark longitudes 0° and 360°, beyond which the distributions have been repeated to help following the spot migration. The long-dashed lines labelled with capital letters from A to E trace the migration of the most conspicuous spots detected in the plots, respectively (see the text).

Fig. 5 The total spotted area as derived from the regularized ME models vs. time for Q = 0.0.

The distributions of the spotted area vs. longitude are plotted in Fig. 4 for the seven mean epochs of our individual subsets adopting a rotation period of 11.953 days. The longitude zero corresponds to the point intercepted on the stellar photosphere by the line of sight to the centre of the star at HJD 2 454 298.55874, i.e., the sub-observer point at the initial epoch. The longitude increases in the same direction as the stellar rotation. This is consistent with the reference frames adopted in our previous studies (Lanza et al. 2009a, 2009b), but does not allow a direct comparison of the mapped active regions with the dips in the light curve.

Our maps show that the distribution of the spotted area versus longitude evolves rapidly, which makes difficult to trace the evolution and migration of the active regions in an unambiguous way. Nevertheless, five main active regions can be identified in Fig. 4 and their migration is traced with dashed lines labelled with letters from A to E. Their rotation rates with respect to the adopted reference frame have been estimated by a linear best fit, i.e., by assuming a constant migration rate, and range from  − 5.9 ± 1.5 deg/day for that labelled with D to 1.1 ± 1.5 deg/day for that labelled with A. Note that active longitudes D and E are not detectable in some of the ME distributions plotted in Fig. 4, but we decided to trace their evolution in a continuous way attributing the change of their visibility to a combination of their fast intrinsic evolution and the differential rotation.

Fig. 6 Upper panel: the observed RV variations due to stellar activity (filled dots) and the corresponding synthesized variations (open diamonds) versus time for the spot models with Q = 0. The vertical dotted lines mark the time interval with the steepest variations that cannot be fitted by our model (see the text). Lower panel: the difference between observed and synthesized variations vs. time.

The relative amplitude of the surface differential rotation, as estimated from the difference between the maximum and the minimum migration rates reported above, is ΔΩ/Ω = 0.23    ±    0.10. Unfortunately, no information is available on the spot latitudes, thus our ΔΩ/Ω is only a lower limit. Since HD 189733 is more active than the Sun, in principle its active regions may cover a latitude range larger than in the Sun where sunspot groups are confined to  ± 40° from the equator (see, e.g., Strassmeier 2009). Fares et al. (2010) using Zeeman Doppler Imaging and adopting a solar-like differential rotation profile, estimated equatorial and polar rotation periods of 11.94 ± 0.16 and 16.53 ± 2.43 days, giving a relative amplitude ΔΩ/Ω = 0.38 ± 0.18. Adopting their differential rotation profile with the above uncertainties, we estimate a latitude ranging from 40° to 80° for the highest latitude spots, i.e., those in the active longitudes B and D. On the other hand, spots in the active longitudes A and C are close to the equator and those in E are at intermediate latitudes. Since the RV variations synthesized with all the spots located close to the equator reproduce the observations well (cf. Sect. 5.3), the lower value for the latitudes of the spots in B, D, and E, i.e.,  ≈ 30°−40°, is favoured. Note that assuming, in the next sub-section, those spots close to the equator leads to an overestimate of their effects in the induced RV variations of  ~20–25 percent which is within our typical errorbars of  ± 5 m s^-1.

The spot longitudes in our model show in general a dependence on the facular-to-spotted area ratio Q (cf., e.g., Lanza et al. 2007), therefore we computed regularized models with the maximum allowed facular contribution corresponding to Q = 5. The relative amplitude of the differential rotation is not significantly modified with respect to the case with Q = 0, i.e., it is still compatible with the result of Fares et al. (2010).

The active longitudes traced in Fig. 4 show significant changes on a timescale as short as 3–4 days although their overall evolution may exceed the duration of the observations, i.e.,  ~30 days. The variation of the total spotted area vs. time is plotted in Fig. 5 and shows a monotonous increase along the MOST observing run corresponding to the overall decrease of the mean flux during the same interval. Note that the area per longitude bin and the total area depend on the spot and facular contrasts adopted for the modelling. Specifically, darker spots lead to a smaller total area while the effect of the facular contrast is more subtle and influences somewhat the longitudinal distribution of the active regions (see, Lanza et al. 2007, 2009a,for more details).

Finally, we note that the time resolution of our spot models is not adequate to look for a possible modulation of the stellar activity with the orbital period of its close-in massive planet, as suggested by Shkolnik et al. (2008). Given the short orbital period of the planet, a different approach should be used to search for signatures of a possible star-planet interaction, as in the case of CoRoT-2 (cf., e.g., Pagano et al. 2009).

5.3. Activity-induced RV variations

To simulate the apparent RV changes induced by the distribution of active regions derived from our light curve modelling, we consider a spectral line with a rest wavelength of 600 nm that is close to the isophotal wavelength of MOST observations for which our contrast coefficients are given. The spectral resolution of the line profile is λ/Δλ = 75   000, i.e., comparable with that of SOPHIE.

We fit the observed RV residuals by adjusting an offset value and the parameter  ⟨ C_d ⟩  that measures the average reduction of the line depth inside starspots. We exclude from the fit the four measurements between HJD 2 454 308.0 and 2 454 310.0 because those residuals show a very steep variation with a decrease of the RV of  ~25 m s^-1 in only two days, that is impossible to reproduce with our model. This fast change and other similar rapid variations in the observed RV time series exceed the accuracy of the SOPHIE velocimetry (cf. Sect. 2.2). They could be due to photospheric velocity fields of several km s^-1 localized in active regions and displaying short-term variations, such as those observed in the Sunduring flares or the emergence of new magnetic flux, or in syphon flows (e.g., Rüedi et al. 1992), or in the time-varying features of the Evershed flows (e.g., Cabrera Solana et al. 2007). We looked at the variability of the residuals in the other published RV datasets (Bouchy et al. 2005; Winn et al. 2006). They also show RV changes on short timescales with night-to-night variations of several tens of m s^-1, in agreement with the present SOPHIE observations (see, e.g., the bottom panel of Fig. 1 in Bouchy et al. 2005). The chromospheric activity, as monitored by the core emission of the Ca II H&K lines, shows remarkable variations on timescales as short as one hour, suggesting frequent flaring events that can contribute to the short-term RV variations (Moutou et al. 2007; Fares et al. 2010).

The best fit for the spot model with Q = 0.0 is found for a RV offset of 21.6 ± 1.0 m s^-1 and  ⟨ C_d ⟩  = 0.25 ± 0.05. We plot the observed and simulated RV residuals vs. time in the upper panel of Fig. 6 together with their differences in the lower panel. The error bars of the synthesized residuals take into account the effects of the photometric accuracy, the differential rotation, and the photometric residuals. The nominal accuracy of MOST photometry is  ~1.1 × 10^-4 in relative flux units that implies an error in the estimated spot area of the order of  ~2.5 × 10^-4 leading to a RV error of  ~0.75 m s^-1 when Eq. (1) in Desort et al. (2007) is applied to compute the maximum deviation. The differential rotation changes the longitudes of the starspots during the time intervals considered for the photometric modelling. Since our modelling assumes a fixed pattern of spots during each interval of 7.625 days, our synthesized RV variations have systematic errors produced by the migration of the active regions during those time intervals. We estimate a maximum deviation of 1.5 − 2.0 m s^-1 depending on the spot migration rate. Finally, the small active regions that evolve on timescales of 1–2 days as revealed by the residuals of the photometric best fit (see Fig. 3 lower panel) also contribute to the uncertainty of our synthesized RV values. Since their estimated areas are of the order of 0.2–0.3 percent of the stellar disc, the estimated fluctuations are of  ~1.8 m s^-1 using the model of Sect. 3.3. In conclusion, the overall standard deviation of our synthesized RV values is  ~2.8–3.0 m s^-1.

The agreement between our synthesized RV variations for Q = 0 and the observations is generally good, although our model cannot reproduce the oscillations shown by the data on timescales as short as 1–2 days because it assumes a fixed configuration of active regions that is changed with a timestep of 3.8 days.Assuming the standard deviation of the RV measurements given by the SOPHIE reduction pipeline, we find a reduced χ² = 1.57 for the best fit of the 23 simultaneous data points. On the other hand, considering the observed RV fluctuations on timescales shorter than 1.0 day, we derive a standard deviation of 5.55 m s^-1 giving a reduced χ² = 1.18.

Including a facular contribution, i.e., considering models with Q > 0, leads to an increase of the reduced χ² of the best fit. This happens because faculae produce an increase of the amplitude of the RV variations in our model leading to values that are systematically below the observations before HJD 2 454 308.0 and systematically above for epochs later than 2 454 310.0. The deviations are of about 10–15 m s^-1 for the spot model with Q = 5, and are already significant for Q = 1 − 2. We conclude that the observed RV fluctuations induced by stellar activity favour a spot model with a negligible facular contribution, supporting the result found in Sect. 4 for the facular-to-spotted area ratio considering only the photometric observations.

Boisse et al. (2011) noticed that most of the power of the RV variations induced by stellar activity is concentrated at the rotation frequency and its first two or three harmonics. Although the number of data points is limited to 23 in our case, we applied their technique to test whether the subtraction of the synthesized RV variations is capable of reducing the impact of stellar activity on the RV measurements. We plot in Fig. 7 the Lomb-Scargle periodogram of the time series of the RV measurements simultaneous with MOST photometry (solid line) and the periodograms of the synthetic RV time series (dash-dotted line, upper panel) and of the residuals obtained after subtracting the synthetic series from the observations (dotted line, middle panel), excluding in all the cases the four points between HJD 2 454 308.0 and 2 454 310.0. All the periodograms have the same normalization to allow an immediate comparison of their power levels. The observed RV variations show the maximum power close to the rotation frequency with the successive harmonics having a slightly increasing power up to the third. This effect is likely due to the presence of multiple, partially unresolved peaks among which the power is split, as clearly evident in the case of the first harmonic. Therefore, the peak of the third harmonic is higher than those of the first and the second because it is narrower. The spectral window of the time series is plotted in the lower panel of Fig. 7 and shows the presence of significant sidelobes, notably at 0.07 d^-1 from the main peak. This, in combination with the remarkable differential rotation of HD 189733, can explain the multi-peaked features in the power spectrum, particularly evident in the signals of the first and second harmonics.

When we subtract the synthesized RV variations, the power is significantly reduced at the rotation frequency, its first, and second harmonics, although some residual power is still present between the fundamental frequency and the first harmonic and between the first and the second harmonics, probably as a consequence of the multiple frequency nature of those peaks in combination with the spectral window. The reduction at the third harmonic after the subtraction of the synthetic RV time series is less remarkable. This is expected because our approach cannot model the variations on timescales shorter than 4 − 5 days. Unfortunately, the small number of simultaneous RV and photometric observations limits our capability to test the performance of our approach. Nevertheless, the present application demonstrates that it can be highly advantageous even for a small sample of data with limited precision, i.e., a photometric precision of a few percent of the flux modulation and a RV precision of 3 − 5 m s^-1.

Fig. 7 Upper panel: Lomb-Scargle periodogram of the observed RV variations produced by stellar activity excluding the four data points between HJD 2 454 308.0 and 2 454 310.0 (solid line) and the periodogram of the corresponding synthetic time series (dash-dotted line); Middle panel: Lomb-Scargle periodogram of the observed RV variations produced by stellar activity excluding the four data points between HJD 2 454 308.0 and 2 454 310.0 (solid line) and the periodogram obtained after subtracting the synthesized RV variations from the observed RV time series (dotted line). The vertical dashed lines mark the rotation frequency and its harmonics for a rotation period of 11.953 days. Lower panel: the spectral window of the time series of the RV observations considered in the upper panel.

6. Discussion and conclusions

We have applied the ME spot modelling method introduced by Lanza et al. (2009a) and Lanza et al. (2009b) to the MOST photometry of HD 189733. It is among the brightest and most active stars accompanied by a transiting hot Jupiter and displays activity-induced RV variations with an amplitude reaching several tens of m s^-1. The simultaneous RV measurements collected with SOPHIE offer us a unique opportunity to compare the models proposed to account for the RV variations induced by activity with the observations. Both the MOST photometry and the SOPHIE RV measurements have errors comparable with those expected with larger telescopes on fainter stars. In the case of CoRoT or Kepler, it is possible to reach a photometric relative accuracy of 1.5 × 10^-4 and 2.0 × 10^-5 on a G2V star of magnitude V = 12 in one hour integration time, respectively. Using HARPS, we can improve the radial velocity accuracy down to 2 − 5 m s^-1 on stars of that apparent magnitude, making our test particularly relevant for those targets, especially for the most active ones, given that HD 189733 is among the most active planetary hosts found so far.

Our spot models show several active regions with lifetimes comparable or longer than the duration of MOST observations, i.e.,  ~ 30 days, that evolve remarkably and migrate at different rates with respect to a reference frame rotating with the mean stellar rotation period. This is interpreted as evidence of surface differential rotation with a relative amplitude of ΔΩ/Ω = 0.23 ± 0.10, in agreement with the result obtained from Doppler Imaging (Moutou et al. 2007; Fares et al. 2010). It is interesting to note that, although the MOST observations cover only  ~ 2.5 stellar rotations, the differential rotation signal is fairly evident thanks to its large amplitude. Moreover, a population of small active regions with a surface of 0.1 − 0.3 percent of the stellar disc and typical lifetimes of 1 − 2 days is observed, as indicated by the oscillations in the residuals of the photometric model (cf. Fig. 4).

We introduce a method to synthesize the RV variations induced by stellar activity from the longitude distribution of the stellar active regions derived by the ME spot modelling. The model assumes that the stellar active regions have the same contrasts and turbulent convection properties of the solar active regions. Starspots are assumed to be close to the stellar equator because we cannot derive their latitudes from the high-precision photometry. Indeed, in the case of HD 189733 the spot latitudes could be estimated from the migration of the active longitudes, given the differential rotation law of Fares et al. (2010), but the results are uncertain and we prefer to compute the RV variations adopting nearly equatorial spots to test our approach in the general case when the latitudinal dependence of the rotation period is unknown.

Our model has only two free parameters, i.e., an offset applied to match the zero points of the synthesized and observed RV scales, and a coefficient  ⟨ C_d ⟩  measuring the average reduction of the depth of the spectral lines in dark spots. They are derived by a best fit to the RV observations. The synthesized RV variations obtained from the starspot distributions show a remarkably good agreement with the simultaneously observed RV variations, although the latter display remarkable changes on timescales as short as 1 − 2 days that cannot be accounted for by our model whose spot distributions can be updated only with a time interval of  ~4 days. This is a limitation of every model based on the rotational modulation of spot visibility because it requires a minimum time interval, not too short in comparison with the stellar rotation period, to map the brightness inhomogeneities on the stellar surface.

The RV variations set also a constraint on the facular effects suggesting that dark spots are dominating both the photometric and the activity-induced RV variations, at least during the time interval covered by the present observations. Of course, this does not exclude the possible presence of bright faculae in the photosphere of HD 189733, but their contribution cannot be detected in the available data. We show that the best value of Q, the facular-to-spotted area ratio, is Q = 0, with an acceptable range extending from  ~0 to  ~5 when only the modelling of the optical flux variations is considered. Considering the modelling of the RV variations, an upper limit of Q ~ 2−3 is derived. In the case of the Sun, the best value of Q is 9 (Lanza et al. 2007). Thus our results indicate a lower relative contribution of the faculae to the light variation of HD 189733. The amplitude of the rotational modulation of our star was  ~0.015 mag during the present MOST observations and  ~0.03 mag during the 2006 observations described in Miller-Ricci et al. (2008), i.e.,  ~5−10 times that of the Sun at the maximum of the eleven-year cycle. This indicates that HD 189733 is remarkably more active than the Sun, which may account for the undetectable facular contribution to its light variations, as suggested by Radick et al. (1998) and Lockwood et al. (2007).

We conclude that the approach presented in this paper can reduce the power of the activity-induced RV variations at the rotation frequency and its low-order harmonics by a factor ranging from 2 to 10, which represents an important improvement for the detection of exoplanets corotating with the star or with an orbital frequency in the 2:1 ratio with the rotation frequency, as suggested by the statistical studies of Lanza (2010) and Lanza et al. (2011). Thanks to the high-precision space-borne photometry already available with MOST, CoRoT, and Kepler, it is possible to map the longitudinal distribution of active regions in late-type stars and apply the method presented in this paper to reduce remarkably the impact of stellar activity on their RV jitter allowing us to confirm the detection of telluric planets or refine the measurement of their mass (see, e.g., the application to CoRoT-7 in Lanza et al. 2010).

Acknowledgments

The authors are very grateful to the MOST team for making available the observations analysed in this work. They wish to thank the anonymous referee for valuable comments that helped to improve this work. Active star research and exoplanetary studies at INAF-Osservatorio Astrofisico di Catania and Dipartimento di Fisica e Astronomia dell’Università degli Studi di Catania are funded by MIUR (Ministero dell’Istruzione, dell’Università e della Ricerca) and by Regione Siciliana, whose financial support is gratefully acknowledged. This research has made use of the ADS-CDS databases, operated at the CDS, Strasbourg, France. I.B. would like to thank the support by the Fundação para a Ciência e a Tecnologia (FCT), Portugal, through a Ciência 2007 contract funded by FCT/MCTES (Portugal) and POPH/FSE (EC), and in the form of grants reference PTDC/CTE-AST/098528/2008 and PTDC/CTE-AST/098604/2008.

References

All Tables

All Figures

Fig. 1 Upper panel: the RV variation produced by a single equatorial spot on HD 189733 vs. the rotation phase for different values of the parameter Cd according to the different linestyles: Cd = 0.0 (no temperature effect on the line profile): solid line; Cd = 0.2: dotted line; and Cd = 0.4: dashed line. Lower panel: the corresponding modulation of the optical flux vs. the rotation phase. Note that the three plots overlap because the modulation is independent of Cd (see the text).

In the text

	Fig. 2 The ratio of the χ2 of the composite best fit of the entire time series of HD 189733 to its minimum value vs. the parameter Q, i.e., the ratio of the area of the faculae to that of the cool spots in active regions. The horizontal dashed line indicates the 95 percent confidence level for , determining the interval of acceptable Q values.
In the text

	Fig. 3 Upper panel: the ME-regularized composite best fit to the out-of-transit light curve of HD 189733 obtained for Q = 0.0. The flux is in relative units, i.e., measured with respect to the maximum observed flux along the light curve. Lower panel: the residuals from the composite best fit versus time.
In the text

Fig. 4 The distributions of the spotted area vs. longitude at the labelled times (HJD−2  454  000.0) for Q = 0.0. The plots have been vertically shifted to show the migration of individual spots (relative maxima of the distributions) versus time. The vertical short-dashed lines mark longitudes 0° and 360°, beyond which the distributions have been repeated to help following the spot migration. The long-dashed lines labelled with capital letters from A to E trace the migration of the most conspicuous spots detected in the plots, respectively (see the text).

In the text

	Fig. 5 The total spotted area as derived from the regularized ME models vs. time for Q = 0.0.
In the text

Fig. 6 Upper panel: the observed RV variations due to stellar activity (filled dots) and the corresponding synthesized variations (open diamonds) versus time for the spot models with Q = 0. The vertical dotted lines mark the time interval with the steepest variations that cannot be fitted by our model (see the text). Lower panel: the difference between observed and synthesized variations vs. time.

In the text

Fig. 7 Upper panel: Lomb-Scargle periodogram of the observed RV variations produced by stellar activity excluding the four data points between HJD 2 454 308.0 and 2 454 310.0 (solid line) and the periodogram of the corresponding synthetic time series (dash-dotted line); Middle panel: Lomb-Scargle periodogram of the observed RV variations produced by stellar activity excluding the four data points between HJD 2 454 308.0 and 2 454 310.0 (solid line) and the periodogram obtained after subtracting the synthesized RV variations from the observed RV time series (dotted line). The vertical dashed lines mark the rotation frequency and its harmonics for a rotation period of 11.953 days. Lower panel: the spectral window of the time series of the RV observations considered in the upper panel.

In the text