1 - INAF, Osservatorio Astrofisico di Catania, Via Santa Sofia, 78, Città Universitaria, 95123, Catania, Italy
2 - Dipartimento di Fisica e Astronomia, Università degli Studi, Via Santa Sofia, 78, Città Universitaria, 95123, Catania, Italy
3 - LAM, Laboratoire d'Astrophysique de Marseille, BP 8, 13376 Marseille Cedex 12, France

Abstract
We analyse the time variability of the total solar irradiance (TSI) as measured by VIRGO/SoHO in order to model the variability of the Sun as a star. Apart from the phases near the minimum at the beginning of activity cycle 23, the period of the rotational modulation is significantly different from the solar synodic period as a consequence of the growth and decay of active regions on time scales shorter than a solar rotation. In order to model the variability of the TSI, we have considered the contributions of discrete active regions and a uniformly distributed background emission. To reproduce the rotational modulation of the TSI, we used three active regions, the areas and coordinates of which were changed every seven days to account for their evolution. The simultaneous presence of dark spots and bright faculae was considered by means of appropriate contrast functions which took into account the observed center-to-limb dependence of their contrast with respect to the unperturbed photosphere. The method proved to be capable of modelling the variability of the TSI on time scales going from 7-10 days up to the solar cycle. The relative amplitude of the residuals was of the order of $(1{-}2)\times 10^{-4}$ with the larger values observed during the phases of maximum solar activity of cycle 23. The application of a similar technique to solar-like stars, such as those that will be observed by the next generation of space-borne photometers, should allow us to minimize the effects of stellar magnetic activity on the detection of planetary transits. Moreover, the availability of long-term highly accurate light curves will allow us to measure stellar rotation period, detect stellar activity cycles, and derive information on the inclination of the stellar rotation axis. The location of the active regions and their irradiance properties can also be retrieved with moderate accuracy from single-band light curves. However, a combination of multi-band photometry and spectroscopy will allow us to constrain some of the free parameters of the model and improve the mapping of stellar surfaces.

Key words: Sun: activity - Sun: rotation - stars: activity - stars: rotation - stars: planetary systems

The surface brightness of the Sun and solar-like stars is a complex function of the position on the disk and the time because of the fluctuations associated with convection, p-mode oscillations and magnetic fields. Specifically, solar active regions, arising from photospheric concentration of magnetic flux, contain both bright faculae and dark sunspots. The time evolution of the magnetic field in each active region, the solar rotation and the long-term changes of magnetic activity level associated with the 11-year cycle make the changes of the flux integrated over the solar disk a very complex function of the time. Radiometers on board of several satellites have monitored the changes of the bolometric flux of the Sun since the late seventies (e.g., Chapman 1987; Fröhlich & Lean 2002). The bolometric flux measurements are usually converted into total solar irradiance (TSI) values, which are independent of distance and other systematic effects, thus revealing the intrinsic variability of the Sun (cf. Fröhlich et al. 1997a). Most of the TSI time variability on time scales of the order of the day or tens of days was explained by models that accounted for the bolometric flux changes produced by the active regions observed on the solar disk by means of empirically derived proxy irradiance functions for sunspots and faculae (e.g., Sofia et al. 1982; Chapman 1987). The modelling of the variability on time scales comparable with the 11-yr activity cycle required the introduction of a uniformly distributed component likely to be associated with the network outside active regions (Foukal & Lean 1988; Foukal et al. 1991).

The recent TSI data from the VIRGO experiment on board the satellite SoHO are of significantly higher quality than those of previous instruments. Moreover, the sunphotometers of the VIRGO experiment have provided us with simultaneous data in three optical passbands centered at 402, 500 and 803 nm with a band width of 5 nm, allowing to study the wavelength dependence of flux changes (Fröhlich et al. 1997a). The VIRGO data nicely confirm the validity of the above mentioned modelling approach, i.e., that the active regions are responsible for most of the variability of the TSI, as well as of the fluxes simultaneously monitored by the sunphotometers, on time scales ranging from a few days up to a few hundred days (see, e.g., Fligge et al. 1998; Unruh et al. 1999; Fligge et al. 2000).

Stars with an outer convective envelope and a sufficient rotation rate to excite a solar-like dynamo are expected to show the same kind of variability of the Sun (e.g., Godoli 1968; Parker 1979). Ground-based optical measurements are limited to an accuracy of a few 10^-3 mag which is about one order of magnitude larger than the typical amplitude of the solar variability. However, several stars with a color index $0.45 \la B-V \la 1.2$ are observed to be variable both on a year-to-year basis and on a seasonal time scale. The correlation with the chromospheric proxy indicators suggests that a solar-like picture is indeed appropriate and that the optical variability can be explained as the result of active region growth and decay and rotational modulation (Radick et al. 1998).

In the forthcoming years the availability of space-borne photometers dedicated to stellar observations will provide us with thousands of optical time series going from one month up to a few years with a sampling rate from a few seconds to a few tens of minutes and an accuracy of 10^-6-10^-5 mag (Baglin et al. 2001; Favata & Aigrain 2002). The main purpose of such instruments is to detect stellar oscillations and to find out Earth-like planets orbiting stars through the periodic light dimming produced by their transits across the disk of the parent star (e.g., Jenkins et al. 1996; Defay et al. 2001; Jenkins 2002).

The planned missions dedicated to stellar photometry include COROT, KEPLER and Eddington. They should be launched between 2005 and 2008. For instance, COROT will perform CCD photometry of stellar fields providing also time series of wide-band color index variations. The solar data are bolometric data, however, because the total solar variability closely mimics the changes in the optical band (e.g., Fröhlich et al. 1997a; Fligge et al. 1998), they can be used as a good proxy for future stellar time series. In such a way, we can test possible techniques for modelling the variability of solar-like stars to reduce the impact of magnetic activity on the detection of planetary transits. Moreover, the analysis of time series should provide information on the rotation period and the inclination of the rotation axis of solar-like stars, as well as on their surface differential rotation and activity cycles, when a long-term coverage is available.

The magnetically active stars studied so far are single objects or components of close binaries with flux variation amplitudes of the order of 0.1 mag in the optical passbands. Their light variations are modelled by assuming huge, dark regions on their photospheres (Vogt 1983; Rodonò et al. 1986; Henry et al. 1995). However, when modelling the variability of the Sun, the activity level of which is much lower, it is necessary to account for the presence of faculae and for a diffuse network. Moreover, the area of the solar active regions is always smaller than 1-2% of the solar surface, so that a simpler approach can be adopted to treat the variability. Its application allows us to model the effects of rotational modulation and active region intrinsic evolution on time scales ranging from about one week up to the solar activity cycle. The method is effective in reducing the amplitude of the solar flux variability, up to an order of magnitude, while leaving unaffected the variability on short scale, due to short-lived active regions with characteristic time scales of about 1-5 d.

In the next section, we present the main characteristics of the VIRGO/SoHO TSI data set and analyse it with an approach we adopted for studying the time series of active stars. In Sect. 3, we introduce our method to model the rotational modulation of the Sun as a star, and in Sect. 4 we present the results of its application. In Sect. 5 we discuss the basic assumption of our method and consider its applicability to other stars. We address its advantages and drawbacks and present a preliminary application to the problem of planetary transit detection. The main conclusions of our study are summarized in Sect. 6.

2 The total solar irradiance data set and its variability

The time series that we analyse in the present study is the level 2.0 series of TSI made available by the VIRGO team on the official web site of the mission. The VIRGO experiment is based on two radiometers DIARAD and PMO6-V, each of which has two channels, one in nearly continuous operation and the other used occasionally to evaluate the degradation of the operating channel due to the exposure to the solar radiation (see Fröhlich et al. 1995 and Fröhlich et al. 1997a, for details). Several corrections must be applied to the measurements to derive the intrinsic solar irradiance and, unfortunately, not all of them are known to the desirable level of accuracy.

The irradiance of the Sun as measured by the radiometers was corrected for the changes in the infrared thermal radiation of the environment within each instrument and then normalized to a reference distance of 1 AU and zero orbital velocity in order to reveal the intrinsic variability of the Sun. The effect of the long-term degradation of each radiometer, which was measured using the occasionally exposed channel, was corrected. This correction was performed for both types of radiometers separately, yielding two independent time series which were then combined into a more reliable estimate of the TSI with a reduced long-term uncertainty. Finally, a comparison with the data of other radiometers and the indication coming from empirical models of the TSI changes were applied to reduce further residual systematic long-term trends. This work led to the so-called level 2.0 time series of the TSI, which is the best approximation currently available for the solar intrinsic bolometric irradiance variations (Fröhlich et al. 1997a; Fröhlich 2002).

The time series analysed in the present study consists of one TSI measurement per hour, each measurement having a relative accuracy of $2.0\times 10^{-5}$ . No obvious outliers are present in the level 2.0 data, in contrast to the DIARAD data set analysed by Jenkins (2002). However, several gaps are present, going from occasional gaps not exceeding a few hours, up to the four month gap due to the loss of the satellite tracking in 1998.

A first analysis of the characteristic time scales of variability contained in the TSI time series can be performed by computing the so-called pooled variance of the data (Dobson et al. 1990; Donahue & Baliunas 1992; Donahue 1993). For a given time interval, the data are binned into segments having the length of that interval and the average variance of the set of these segments, that is the pooled variance, is computed. The pooled variance is plotted as a function of the time interval length in Fig. 1, where four ranges of variation can be identified.

$\begin{figure} \par\includegraphics[width=7.3cm,clip]{3322.f1}\end{figure}$

Figure 1: The pooled variance of the relative variations of the TSI as a function of the time interval $\tau$ chosen to bin the data. The relative variations are computed with respect to a reference value S₀= 1367.3 W m^-2. The main ranges of variation of the pooled variance are indicated with A, B, C and D (see the text).

The pooled variance analysis is independent of the shape of the flux time variation. Since we are interested in modelling the rotational modulation, it is useful to apply a different tool to look for a periodicity of the data in the characteristic time range of solar active region rotation. The latitude variation of the activity belts and the different mean lifetime of the surface features, as a function of the activity cycle phase, produce a variation of the period of the rotational modulation along the 11-yr cycle (e.g., Donahue & Keil 1995). This effect is worth to be taken into account since our time series goes from January 1996 to September 2001, i.e., from the activity minimum, at the beginning of cycle 23, up to its maximum and beyond. A useful tool to study the time variability of the rotational modulation is provided by the Morlet wavelet, as described in detail by Hempelmann & Donahue (1997) and Hempelmann (2002). We computed the Morlet wavelet amplitude for the entire TSI time series by adopting the definition of Hempelmann & Donahue (1997). The relative amplitude of the wavelet transform is plotted in Fig. 2, together with the data sequence, so that a readily straightforward comparison between them is possible.

$\begin{figure} \par\includegraphics[angle=90,width=17.2cm,clip]{3322.f2}\end{figure}$

Figure 2: Upper panel: the amplitude of the Morlet wavelet of the TSI time series normalized to its maximum value. Lower panel: the TSI time series of level 2.0 as derived by the measurements of VIRGO/SoHO described in Sect. 2.

During the phase of solar cycle minimum (1996-1998), the peak amplitude of the wavelet is $\sim$ 40% of its maximum, but the period is well defined at about 27.7 d, which corresponds to the synodic rotation period of the Sun in the latitude range of the sunspot groups. The large data gaps in 1998 are likely to be responsible for the period shifts and for the splitting observed from 1998.0 to 1999.0, so that the wavelet analysis in this time interval can not be trusted. Starting from 1999.1, the data continuity is restored. The results, however, can be regarded with confidence only from mid 1999, because the wavelet feels the effects of the gap up to 150-200 d after it ended. In the time interval from 1999.5 to the beginning of 2001, the wavelet transform amplitude is about two times larger than during the minimum phase, because the amplitude of the rotational modulation due to active regions is now larger. However, the simple structure of the wavelet transform is now lost. There is a splitting of the amplitude into two ridges, one centered around a period of 26.2 d and another around 31 d, with the latter having the largest amplitude. Such a splitting can be interpreted as a consequence of the phase shifts produced by the growth and decay of activity complexes at different longitudes, and over times scales shorter than the time resolution of the Morlet wavelet, as discussed by Hempelmann (2002). Therefore, when the Sun is near the maximum of the 11-yr cycle, the rotational modulation signal can be masked by the active region evolution. A similar conclusion was reached by LaBonte (1982, 1984) by the analysis of the absolute value of the magnetic flux and the 2.8 GHz solar flux time variations, the latter being a good tracer of sunspot groups. The rotation period as derived by means of standard Fourier analysis techniques applied to 1 yr data sets, turned out to be highly variable with values differing from the true solar synodic period up to about 6 d. Donahue & Keil (1995) pointed out that a crucial parameter for the proper determination of the rotational modulation is the length of the data set over which the Fourier spectrum is computed and suggested that it should not exceed 150-200 d for the Sun in order to detect systematic changes related to the evolution of the activity cycle. In view of such considerations, we have re-analysed the TSI time series after splitting it into thirteen subsets of length 158 d. In order to attain a better estimate of the rotational modulation period and of the confidence level of the detected periodicity, we applied the Lomb-Scargle periodogram method to each subset (Scargle 1982; Horne & Baliunas 1986). The results are listed in Table 1, where the MJD of the center of each subset (the first column), the year (the second column) and the period $P_{\rm mod}$ corresponding to the main peak of the periodogram of the data in the subset, are given. The confidence level of all the peaks is larger than 99.9%. The effective resolution of the period determination is of the order of 1.0 d, except for the sixth and the seventh subsets for which the large data gaps degrade the resolution by a factor of about two.

Table 1: The period corresponding to the highest periodogram peak in each of the subsets of 158 d derived by dividing the entire TSI series into 13 time intervals.
MJD Year $P_{\rm mod}$

(JD - $2.4\times 10^{5}$ ) (d)

50180.77 1996.27 26.2

50338.77 1996.70 26.7

50496.77 1997.13 27.2

50654.77 1997.56 31.4

50812.77 1998.00 30.3

50970.77 1998.43 29.3

51128.77 1998.86 25.0

51286.77 1999.30 24.5

51444.77 1999.73 28.5

51602.77 2000.16 29.9

51760.77 2000.59 31.9

51918.77 2001.03 22.0

52076.77 2001.46 26.2

**Table 1:** The period corresponding to the highest periodogram peak in each of the subsets of 158 d derived by dividing the entire TSI series into 13 time intervals.
MJD	Year	$P_{\rm mod}$
(JD - $2.4\times 10^{5}$ )		(d)
50180.77	1996.27	26.2
50338.77	1996.70	26.7
50496.77	1997.13	27.2
50654.77	1997.56	31.4
50812.77	1998.00	30.3
50970.77	1998.43	29.3
51128.77	1998.86	25.0
51286.77	1999.30	24.5
51444.77	1999.73	28.5
51602.77	2000.16	29.9
51760.77	2000.59	31.9
51918.77	2001.03	22.0
52076.77	2001.46	26.2

3 The model of the solar rotational modulation

In order to model the variation of the solar irradiance due to active regions, we adopted the approach described by Oster et al. (1982) and Sofia et al. (1982). The difference between the solar irradiance in the presence of active regions S and a reference value S₀ is given by:

Solar observations show that the contrast of faculae and, to a lesser extent, of the sunspots varies significantly among active regions (Chapman 1987; Chapman et al. 1994). Moreover, systematic changes with the phase of the solar cycle are also possible (Albregtsen et al. 1984). In particular, the average value of Q is a function of the phase of the solar cycle, with Q decreasing for larger active regions appearing around the maximum of activity (Foukal 1998). We did not include these effects in our model because we are interested in analysing the TSI by adopting a simple stellar-like approach. In the case of the other stars, we have no a priori information on the phase of the activity cycle or the faculae-to-cool-spots area ratio.

The value of the reference bolometric irradiance S₀, which corresponds to the solar flux without active region and uniform network contributions, is also unknown in the stellar case. S₀ is the value with respect to which the model TSI is normalized before fitting the observed variations. In modelling stellar wide-band light curves it is customary to adopt as a reference value the flux corresponding to the brightest magnitude of the star in the considered time series, i.e., at the minimal coverage of the stellar disk by cool spots (cf., e.g., Messina et al. 1999). In the present analysis, we adopt a stellar-like approach and assume S₀ = 1367.3 W m^-2, i.e., the maximum TSI value in our time series. It is important to note that the values of the active region areas A_i's and of the uniformly distributed network contribution $S_{\rm r}$ are referred to such a value of the solar TSI. However, their variations are largely independent of S₀ and represent the truly meaningful quantities to characterize the stellar activity behaviour, as will be shown by our modelling results (cf. Sect. 4).

The position of an active region on the Sun is defined by its longitude $\lambda_{i}$ and latitude $\theta_{i}$ in the Carrington reference frame. The value of $\mu_{i}$ is given by:

It is important to note that the angular velocity of rotation of the reference frame of our model $\Omega$ does not coincide with that of the Carrington reference frame. Actually, it is an additional free parameter of the model because the rotational modulation of the solar flux by the active regions does not show a constant periodicity during the time interval we analysed, as shown in Sect. 2.

An important point regards the choice of the number of active regions to be adopted to model the TSI rotational modulation. A simple semi-sinusoidal modulation can be modelled by only one active region, but the observed modulation is usually of considerably more complex shape. This implies that at least two or, more frequently, three active regions are needed to obtain a good fit of the TSI modulation.

The analysis of the most accurate stellar light curves obtained so far supports such a choice also in the case of magnetically active rotators (Budding 1995). Moreover, the validity of the assumption of three active regions is confirmed a posteriori by the goodness of fit to the entire TSI time series (cf. Sect. 4). Therefore, our model requires the specification of eleven free parameters, that is the area and the coordinates of each active region plus the background irradiance $S_{\rm r}$ and the period of rotation of the pattern of the active regions $P\equiv 2\pi/\Omega$ . Their best values are evaluated by fitting the model to the observed TSI, $S_{{\rm obs}}$ , by minimizing the $\chi ^{2}$ statistics:

In order to fit the observations with our model, it is important to consider a data subset spanning a time interval shorter than the time scale for a significant evolution of the active regions on the Sun. After a series of tests, we found that a time interval of 14.0 days is the best compromise between the requirement of avoiding significant active region changes and that of sampling adequately the rotational modulation of the active regions themselves. As a matter of fact, the evolution of the active regions is so fast on the Sun that it is usually impossible to model adequately an entire rotation by using a fixed active region pattern, whereas with our choice of the length of the subset interval the fit turned out to be always good.

It is important to note that a numerical instability may arise in the fitting of a given data set because a linear combination of active region longitude and angular velocity appears in Eq. (2). Therefore, we have always fitted each data set by fixing a test value for the rotation period and then finding the best values of the remaining parameters that minimized the $\chi ^{2}$ by applying the constrained minimization algorithm by Byrd et al. (1995), as implemented by Zhu et al. (1994). Then the rotation period P was iteratively adjusted in the interval 23.0-33.5 d by performing successive fitting to locate the absolute minimum of the $\chi ^{2}$ by means of a golden section algorithm (Press et al. 1992).

We used the best fit parameters of a given data set as the starting point for the analysis of the successive set, by adopting a spacing between consecutive sets of 7.0 d, unless a larger gap was imposed by the lack of data. In this way we could follow the evolution of the pattern of active regions on a time interval shorter than that adopted to sample their rotational modulation. This procedure led to a further improvement of the quality of the fit.

It is worth noticing that the parameters we held fixed specify the limb-darkening of the unperturbed stellar disk and the contrast of the surface inhomogeneities. They are ill-constrained by bolometric flux variations alone. In fact, only multi-band data may provide information for their determination (cf. Sect. 5). On the other hand, the geometrical parameters specifying the active region areas and coordinates and the uniform network contributions are better-constrained by the TSI modulation and are assumed as variable parameters in our fitting procedure.

From a mathematical point of view, the solution obtained by minimizing the $\chi ^{2}$ is usually unique when the TSI time variation is dominated by the rotational modulation of a few and stable active regions, and the contrast and limb-darkening parameters are known a priori as well as the inclination of the rotation axis and the rotation period (e.g., Eker 1996). Conversely, when the TSI variation is produced by a complex pattern of many small active regions or is greatly affected by the fluctuations due to active region growth and decay, multiple solutions may be found and the convergence and stability of the fitting algorithm may not be warranted.

Our procedure of fitting the geometrical parameters of the active regions only for fixed values of P and then minimizing the $\chi ^{2}$ with respect to P, eliminates a primary source of instability in the convergence procedure. However, in order to obtain a locally stable solution, we found that also the values of the active region areas A_i's and the contribution of the uniformly distributed network $S_{\rm r}$ must be constrained from the outset. Otherwise, there can be a trade-off between the active region total area and the variation of $S_{\rm r}$ , which makes the resulting solution non-unique and unstable. This is a well known inconvenient in the analysis of wide-band light curve sequences and we adopted the same method used in stellar light curve analysis to overcome the problem. Specifically, the depth of the largest dip of the TSI time series was used to fix limits for the spot area: $0\leq A_{i} \leq 8.2\times 10^{-4}$ ; while the range spanned by the upper envelope of the TSI series in the 1996-2001 interval was used to constrain the variation of $S_{\rm r}$ : $-1.5\times 10^{-4}\leq S_{\rm r}/S_{0} \leq 1.5\times 10^{-4}$ .

The above mentioned constraints proved to be sufficient to warrant convergence and local stability for the determination of the geometrical free parameters in the case of three active regions. However, from a global point of view, sometimes multiple solutions with comparable values of the $\chi ^{2}$ may be present and the fitting algorithm may become unstable when their minima are so close that their basins of attraction intersect one with the other. Actually, we found that our best fitting procedure converged in more than 80% of the cases by adopting as an initial guess, for the analysis of a given TSI subset, the best fit parameters of the previous subset. In the remaining cases, we needed to change the initial values of the parameters to avoid the solution to get stuck around a local minimum with an unacceptably high value of the $\chi ^{2}$ . In these cases, by applying repeatedly a bootstrap routine for the choice of new initial parameters (i.e., starting values of the active region areas and longitudes), we searched for other minima in the $\chi ^{2}$ landscape until a solution with an acceptable $\chi ^{2}$ value was found. Such a multi-stage procedure yielded a good fit in all the cases, i.e., for all the analysed TSI subsets ranging from mid January 1996 to the end of August 2001.

$\begin{figure} \par\includegraphics[angle=90,height=14.3cm,width=13.2cm,clip]{3322.f3}\end{figure}$

Figure 3: a) The TSI data sequence for the time interval from mid-January 1996 to June 1997 (open diamond) and the best fit computed with our model (solid line); b) the deviation of the best fit from the TSI data in the same time interval; c) the TSI data sequence from June 1997 to October 1998 (open diamond) and the best fit (solid line); d) the deviation of the best fit from the TSI data in the same time interval; e) the TSI data sequence from October 1998 to March 2000 (open diamond) and the best fit (solid line); f) the deviation of the best fit from the TSI data in the same time interval; g) the TSI data sequence from March 2000 to August 2001 (open diamond) and the best fit (solid line); h) the deviation of the best fit from the TSI data in the same time interval. The model is always embedded into the data sequence except when gaps are present (cf., e.g., the interval 1999.12-1999.13).

$\begin{figure} \par\includegraphics[height=14.3cm,angle=90,width=13.9cm,clip]{3322.f4}\end{figure}$

Figure 4: A zoom-in to eight short time intervals along the whole data sequence plotted in Fig. 3. In each couple of panels, the upper one shows the TSI data (open diamonds) and the best fit computed with our model (solid line) vs. time, whereas the lower one shows the deviation between data and fit vs. time. The sequence of panel couple runs from activity cycle minimum (upper left couple) to the rising phase of cycle 23 (lower left couple) and then from the rising phase (upper right couple) to activity maximum of cycle 23 (lower right couple). The model is always embedded into the data sequence, except when data gaps or a significant deviation is present (cf., e.g., the 1996.89-1996.93 time interval in the second couple of panels on the left column).

4 Results

The observed relative variation of the TSI together with our model best fit is plotted in Fig. 3 with a zoom-in to selected time intervals shown in Fig. 4. The value of the fit for a given datum was computed with the parameters determined by the best fit of the 14-d data subset the midpoint of which was the closest preceding the datum. Thus the time interval between a given datum and the nominal epoch of the best fit parameters was always smaller than 7.0 d. Of course, successive fits had slightly different endpoints, but their relative differences did not exceed $5\times 10^{-5}$ and were randomly distributed. We chose to retain such small differences and the small changes in the slopes of the fitting curves near the endpoints to avoid the introduction of any forced matching in the data space, which should have been difficult to justify on a geometrical or physical ground. In any case, the cadence of the endpoint mismatchs can be easily changed by modifying the time interval between successively fitted data sets. The fit computed according to these prescriptions was always good with relative deviations from the data <10^-4, during the activity minimum phases, and of $\sim$ $2\times 10^{-4}$ during the phases with the largest active regions, close to cycle 23 maximum. Such deviations were due to the growth and decay of small active regions on time scales of the order of a day that could not be modelled by our approach because the active region configuration we adopted was changed every seven days.

The $\chi ^{2}$ of the best fits to the individual data subsets were plotted in Fig. 5. Most of the values were around the unity. However, significantly larger values of the $\chi ^{2}$ were obtained for the data sets with the largest active complexes around the activity maximum. This is likely to be due to the evolution of several short-lived active regions that could not be modelled by our approach.

The best values of the rotation period, as derived by means of our fitting procedure, are plotted in Fig. 6. They show a large dispersion around a mean value of 30.25 d and no systematic trend with the solar cycle is apparent. A time baseline of fourteen days is certainly not adequate to measure the rotation period with high accuracy. Moreover, the large variations between the fits to successive data subsets suggest that our code tried to deal with the growth and decay of the active regions by changing the period of rotation of the model star in order to vary their projected area during the transit over the disk and their epochs of rising and setting at the limb. The fact that the average period is about 10% longer than the synodic period of rotation of the Sun is due to the accumulation of the period values toward the upper limit for $P_{{\rm rot}}$ that was fixed at 33.5 d. This indicates that, on the average, the photometric effect of a model active region tends to persist for a somewhat longer time than the actual rotation period of the Sun, an effect that may be related to our assumption that spots and faculae were co-spatial and that the extension of the active regions was negligible.

It is interesting to notice that the clustering of the points in Fig. 6 around some discrete values is an artifact produced by the golden section method used to search for the best period. The accumulation toward the upper limit is likely related to the persistence of the photometric effect of the active regions mentioned above. It seems more pronounced in the interval 1999-2001, when the Morlet wavelet transform had a significant amplitude for periods well in excess of the solar synodic period. We notice that surface differential rotation alone can not account for the range of the rotation period because its amplitude is of only $\sim$ 5% in the latitude range of solar active regions.

$\begin{figure} \par\includegraphics[width=7.3cm,clip]{3322.f6}\end{figure}$

Figure 6: The solar rotation period from our fitting procedure described in Sect. 3 for each of the fourteen-day TSI data subset vs. time.

The total sunspot area as determined from our model is plotted in Fig. 7 together with the observed monthly mean sunspot group area compiled by Hathaway (2002). The model values follow rather well the overall variation trend of the observed values. It is important to note that the observed areas of the sunspot groups reported in different compilations differ systematically from each other, even if approximately constant conversion factors between different compilations can be assumed (e.g., Hathaway 2002). The lower panel of Fig. 7 shows that the model and observed area values are reasonably well correlated. The ratio of about 2 between the model areas and the observed areas during the rising phase of cycle 23 and close to its maximum could be a consequence of the different weight of the spot penumbra and different methods in measuring the group areas. However, a larger systematic deviation is apparent close to the minimum, when the model spot areas are remarkably larger than the observed ones. Such a behaviour is a consequence of the fixed ratio between the facular and the sunspot areas in the model. As a matter of fact, the rotational modulation of the TSI during the phases of minimum activity was dominated by faculae with only a very small contribution from sunspots. In our model, a fixed proportion is assumed between facular and spotted areas, which leads to an overestimate of the sunspot area during those intervals in which the actual ratio of facular to sunspot area is larger than assumed. Since this ratio varies inversely with the level of activity (Foukal 1998), the model best reproduces the observations during the epochs of maximum activity, when Q=10 should have been a more adequate assumption (Chapman et al. 1997; Fligge et al. 1998; Chapman et al. 2001).

$\begin{figure} \par\includegraphics[height=11.3cm,width=7.4cm,clip]{3322.f7}\end{figure}$

Figure 7: Upper panel: the monthly average of the observed total sunspot group area vs. time from Hathaway (2002) (solid line) and from our best fit model (dotted line). The areas are in unit of the solar surface. Lower panel: the total best fit sunspot area vs. observed sunspot group area. The dashed line indicates a constant ratio of 2.0 between best fit model and observed areas.

$\begin{figure} \par\includegraphics[width=7.8cm,clip]{3322.f8}\end{figure}$

Figure 8: The value of the model parameter $S_{\rm r}$ , as derived from the best fit of the fourteen-day data subsets vs. time.

The contribution to the TSI of structures uniformly distributed in longitude was modelled by the term $S_{\rm r}$ in Eq. (1). As shown in Fig. 8, it was always negative. Apart from some isolated values, its variation mimics that of the upper envelope of the rotational modulation of the TSI, as expected on the basis of our model assumptions. As a matter of fact, the variation of the average TSI level along the five years covered by our data set could be modelled by an increase of a uniform facular contribution toward the activity maximum or, alternatively, by an increase of a uniform sunspot contribution around the activity minimum. It was not possible to distinguish between such two scenarios on the basis of bolometric data alone. We fixed the reference value S₀ equal to the TSI maximum, which led to negative values for $S_{\rm r}$ in order to account for the lower average irradiance level close to the minimum and during the rise to the cycle maximum.

$\begin{figure} \par\includegraphics[angle=90,width=14.8cm,clip]{3322.f9}\end{figure}$

Figure 9: The latitude of the observed sunspot groups (filled circles) and of the three active regions from our best fit model to the TSI variations (open triangles) vs. time. The dimensions of the symbols are proportional to the area of the sunspot group or active region, respectively.

The model and observed active region latitudes of sunspot groups are compared in Fig. 9. The typical structure of the butterfly diagram is clearly traced by the observed groups, with a remarkable increase of the latitude at the beginning of 1997, when cycle 23 started after the 1996 minimum. The latitudes of the model active regions are systematically higher than those of the observed sunspot groups. The rotational modulation of the TSI did not constrain adequately the latitude of the surface inhomogeneities because the inclination of the solar rotation axis was always close to $90^{\circ}$ implying that the duration of the transit of a surface feature across the disk was almost independent of its latitude. An active region at high latitude produces a photometric contribution dominated by the faculae because they are more contrasted toward the limb than spots. Hence, large active regions at high latitudes could be used to effectively increase the contribution of faculae with respect to sunspots in our model, as needed during the activity minimum. Conversely, during the maximum of activity, the modulation of the TSI was dominated by cool spots and most of the model active regions fell at low and intermediate latitudes in better agreement with the latitudes of the observed groups. Such a result may be a consequence of the fact that Q=10 was more appropriate for phases of maximum activity. We note that a small active region near the visible rotation pole could be used to finely tune the fit, particularly during phases of minimum activity, when a smaller projected sunspot area was needed.

The situation is different for the longitude of the active regions because the rotational modulation of the TSI can effectively constrain the variation of the surface brightness of the Sun with longitude. In Fig. 10 the longitudes of the observed sunspot groups, together with the longitude of the model active regions, are plotted for a time interval of about one year during the rise toward the maximum of cycle 23. We excluded active regions at latitudes > $60^{\circ }$ because their longitudes were ill-defined since the inclination of the rotation axis of the Sun is close to $90^{\circ}$ . Taking into account that the time extension of the fitted subsets was 14.0 days (i.e., 0.038 yr), the general correlation between the longitudes of the observed and the modelled active regions was quite satisfactory, with differences usually not exceeding $\pm$ $15^{\circ}$ . However, the number of sunspot groups actually present on the Sun at a given time could be larger than three and significant modifications of the sunspot groups took place on time scales shorter than one week, the time interval between successively modelled data sets. The combination of these effects produced in several cases a large deviation between model and observed active region longitudes, because our code preferred to place active regions between large sunspot complexes in an attempt to fit the flux modulation arising from a more complicated distribution of surface features.

It is interesting to compare the power spectrum of the observed TSI variations with that of the residuals, after our best fit was subtracted from the whole data set. The results are presented in Fig. 11 for the five time intervals that are labelled, considering a frequency range from a 60-d periodicity up to the Nyquist frequency of our data sampling (one value per hour). The amplitude of the TSI fluctuations for frequencies corresponding to periods between 10 and 60 d in the phase of maximum activity of cycle 23 (1999-2001) was larger by about a factor of two with respect to minimum activity (1996-1997), because the amplitude of the rotational modulation during the former interval was correspondingly larger. When the rotational modulation was modelled following our approach, the amplitude of the residual fluctuations was remarkably reduced by a factor of $\approx$ 20-30 in the same frequency range.

5 Discussion

Our model for the rotational modulation of the Sun as a star proved to be successful in fitting most of the solar variability in the range from 5-7 d up to the 11-yr solar cycle time scale.

$\begin{figure} \par\includegraphics[angle=90,width=7.5cm,clip]{3322.f10}\end{figure}$

Figure 10: The longitude of the observed sunspot groups (filled circles) and of the three active regions found by the best fit model to the TSI variations (open triangles) vs. time; the dimensions of the symbols are proportional to the area of the sunspot groups or of the active regions, respectively. Only active regions with latitudes $\pm$ $60^{\circ }$ from the equator are plotted.

$\begin{figure} \par\includegraphics[angle=90,width=15cm,clip]{3322.f11}\end{figure}$

Figure 11: The power spectrum of the relative variation of the TSI (solid line) and of the residuals with respect to our best fit models (dotted line), in the labelled time intervals, vs. frequency.

A comparison with the observed areas of the sunspot groups showed that our approach is capable of determining the long-term variation of the total spotted area, although the values found by the model differ systematically from the observed ones because of the model assumptions on spot irradiance, the fixed facula-to-spot area ratio, the neglected penumbra effects and, least but not last, the different methods adopted to estimate the observed group area. The longitudes of the active regions, as derived from our model, compared well with those of the observed sunspot groups, although the simplification introduced by the assumption of only three active regions produces significant differences in the case of a pattern consisting of several groups. The fixed ratio between facular and spotted area may be an additional source of systematic errors. However, even in the ideal case of a noiseless light curve, numerical simulations show that it is not possible to go beyond a resolution of $\sim$ $70^{\circ}{-}90^{\circ}$ in longitude, when the surface pattern consists of a continuous distribution of brightness inhomogeneities (cf., e.g., Lanza et al. 1998). Therefore, a number of model active regions greater than three would generally not improve the effective map resolution, but increases the number of free parameters giving rise to non-uniqueness and instability of the solution. On the other hand, solutions with only two active regions sometimes provided us with good fits. However, the coordinates of the model active regions were in better agreement with those of the observed active regions only when the surface pattern was dominated by one or two large regions, as happened during a few rotations close to the minimum of activity (cf., e.g., Oláh et al. 1999). In general, the fit provided by models with three active regions turned out to be significantly better during most of the analysed time sequence.

The length of the fitted subsets was fixed to 14.0 d, as the best compromise between the opposite requirements for sampling the rotational modulation and avoiding sizeable active region evolution. Longer sequences could be fitted by increasing the number of active regions or by allowing their area to change versus time. Again, this would imply an increase of the number of free parameters which leads to the above mentioned problems.

It is interesting to compare the results of our analysis with those of previous modellings of the solar rotational modulation. In an attempt to fit the chromospheric emission of the Sun as a star, Catalano et al. (1998) found a good agreement between the longitude of the model active regions and the H $\alpha$ faculae by using a fitting procedure with Maximum Entropy regularization based on the method introduced by Lanza et al. (1998). The main differences with the present approach consist in: a) the use of a continuous brightness distribution instead of a set of discrete brightness inhomogeneities; b) the presence of only one kind of surface inhomogeneities, i.e., bright chromospheric faculae. The Maximum Entropy regularization yielded active region areas systematically larger than those of the H $\alpha$ faculae, although such a result could depend also on the contrast adopted for the faculae. The presence of only one kind of brightness inhomogeneities led to a remarkable simplification, avoiding the systematic errors associated with the simultaneous presence of faculae and spots in our model. Therefore, the distribution of the chromospheric faculae in longitude could be reliably reconstructed by the modelling, although with a limited resolution. Similar conclusions were reached by Oláh et al. (1999) who modelled the rotational modulation of the Sun as a star in the microwaves and in soft X-rays during a phase of minimum activity characterized by a few observed active regions. Their results gave further support to a modelling approach based on discrete active regions and clarified some systematic effects due to the intrinsic evolution of the active regions.

A basic source of uncertainty comes from the fact that the bolometric flux modulation does not contain enough information to fully constrain the bright and dark components of our model. In principle, it is possible to fit the bolometric variation by using only one kind of inhomogeneities, for instance dark spots or, alternatively, bright faculae. Only multiband data can provide us with the additional information required to model both spots and faculae without ambiguity. As a matter of fact, VIRGO/SoHO provided this kind of information through its sunphotometers that monitored the solar irradiance changes in three bands of width 5 nm centered at 402, 500 and 802 nm. Fligge et al. (2000) demonstrated how such additional information could be used to constrain effectively the center-to-limb variation of the facular contrasts which, in turn, arised from the thermal stratification in the facular atmosphere.

In the case of the stars to be observed by space-borne telescopes, wide-band colors should be available to effectively contrain the spotted and facular components of the model. In principle, bolometric data may be used to find the position of the active regions, whereas simultaneous color data may allow us to evaluate the contrast of the inhomogeneities, which depends on their thermal structure. Colors in the optical and UV bands can provide information on the temperature of the hot faculae, as shown by, e.g., Fligge et al. (2000), while infrared colors can give the temperature of cool spots (cf., e.g., Poe & Eaton 1985; Strassmeier & Oláh 1992). The sections of the light curve that are dominated by the rotational modulation of a single active region should be preferentially used for such temperature determination in order to reduce the ambiguity that may arise from the simultaneous presence of several active regions.

The limb-darkening parameters can be derived from state-of-the-art model atmospheres (e.g., Diaz-Cordoves et al. 1995; Kurucz 2000) if information on the effective temperature and surface gravity of the star is available. In addition to the methods based on wide-band colors, the contrast of cool spots may be evaluated by the same kind of models if the spot temperature is known. It can be estimated in very active stars by measuring the intensity of the TiO bands (Neff et al. 1995) or, in moderately active rotators, by the ratio of the depth of temperature sensitive lines (e.g., Catalano et al. 2002). Such a method may be hopefully applied to stars with low activity levels, more similar to the solar case.

The average atmospheric model of the faculae is another important ingredient because its thermal stratification gives rise to center-to-limb variations of facular contrast. Moreover, the facular-to-spot area ratio must be determined for distant stars. It is important to note that such parameters can be reliably estimated from spectroscopic observations, only if they cover a suitably long time interval. The rotational modulation of the flux in the core of the Ca II H&K lines can be used to estimate the coverage factor of the faculae during a given rotation, but the evaluation of the Q parameter variation requires a data set covering at least a complete activity cycle. This is difficult or impossible to obtain for samples as large as those that will be observed by the next generation of space-borne stellar photometers. However, the important result that comes from our study is that the rotational modulation can be fitted even if the true value of Q is different from the adopted one.

In view of such considerations, it may be possible to use a simpler approach to estimate the parameter Q for stars other than the Sun, by extrapolating the solar behaviour to higher levels of activity, as suggested by Solanki (1999). If such an approach is correct, more active stars should have a smaller value of Q. This implies a larger photometric weight of cool spots, as suggested by the observations of Radick et al. (1998) and it could help to simplify the model and obtain more reliable parameters by fitting only single-band modulations.

In the modelling of solar-like stars additional unknown parameters are the period of rotation, the inclination of the stellar rotation axis and the unperturbed light level.

The rotation period can be accurately determined only if the lifetimes of the surface features is longer than the rotation period itself, a circumstance that favour stars that rotate more rapidly than the Sun and have larger long-lasting active regions. We showed that the pooled variance of the time series may give information on the time scales of the different sources of flux variations, allowing us to constrain the range of the rotation period values to the purpose of our modelling. The method can be successfully applied not only to stars with a level of activity comparable to the Sun, but also to highly active rotators (cf. Messina et al. 1999; Rodonò et al. 2000).

By fitting some synthetic time series with two or three active regions, we found that the inclination i of the rotation axis could be retrieved with an accuracy of $\Delta i \simeq \pm 20^{\circ}{-}25^{\circ}$ , if the rotation period is known with an accuracy better than 1-2%, the inclination lies in the range $70^{\circ} < i < 90^{\circ}$ , the active regions do not change, and the noise in the data set does not exceed 10^-6 -10^-5 of the measured flux. However, in the case of real time series, the uncertainty on the period of rotation, the intrinsic active region evolution and the complexity of their pattern will make this $\Delta i$ estimate certainly too optimistic. A sizeable error in the inclination would produce a systematic error in the latitude of the active regions, but should not affect too much their longitude, neither the quality of the fit. Only in the case of a small inclination angle ( $i \la 20^{\circ}$ ), the fitting parameters may become unstable (cf. Rodonò et al. 1986; Kovári & Bartus 1997). Moreover, as pointed out by Schatten (1993), the higher contrast of faculae toward the limb may produce an enhanced variability of solar-like stars observed at low inclination angles.

The unperturbed light level is usually unknown unless a very long photometric time series is available encompassing at least a few activity cycles. However, one can assume the brightest magnitude in the available data set, a choice which would certainly introduce some systematic effect on the value of the fitting parameters, especially on the A_i's and $S_{\rm r}$ , but should not significantly limit the information that can be derived from the variation of the other parameters.

In the light of the results of Sect. 4, our modelling approach may significantly improve planetary transit detection in the case of stellar variability dominated by rotational modulation. As an illustrative example, we show in Fig. 12 the TSI variation with superposed the transit of an Earth-like planet of radius $R=2.3~ R_{\oplus}$ , for which the central amplitude is $4.0\times 10^{-4}$ of the average solar flux. After the subtraction of the best fit to the TSI rotational modulation, the transit feature is clearly evident, even during epochs of maximum solar activity.

We notice that the power of the sharply-shaped periodic transit dips is scattered in the Fourier space over a large number of harmonics. Therefore, the same signal can not be directly detected in the power spectrum of the original TSI time series because the amplitude of its fluctuations was about 20-30 times larger (cf. Fig. 13).

We also notice that the fluctuations with time scales shorter than $\sim$ 5 d can not be reproduced by our approach. Therefore, their amplitudes are not significantly modified in the residual time series. This implies that the shape of any transit light curve is not greatly affected by the subtraction of the best fit model since transit duration around a solar-like star ranges from a few hours to a few tens of hours.

$\begin{figure} \par\includegraphics[width=7.5cm,clip]{3322.f12}\end{figure}$

Figure 12: a) The TSI variation with superposed the transit of an Earth-like planet of radius $R=2.3R_{\oplus }$ and orbital period of 30.0 d; b) The same time series after subtracting the best fit model light curve computed with the present model.

$\begin{figure} \par\includegraphics[width=7.6cm,clip]{3322.f13}\end{figure}$

Figure 13: a) The power spectrum of the relative variation of the TSI with the transit of a planet with $R=2.3R_{\oplus }$ and a period of 30.0 d superposed (cf. Fig. 12a); b) The power spectrum of the same data sequence after the subtraction of the best fit computed with the present model (cf. Fig. 12b). The peaks corresponding to the first four harmonics of the transit signal are now clearly detectable.

It is interesting to analyse in more detail the power spectrum of the TSI variation to investigate the variability on short time scales, i.e., $\leq$ 5 d, that are not accounted for by the rotational modulation. In Figs. 14 and 15 we plotted the amplitude of the power spectrum of the TSI and of the residuals of the fit during the activity minimum phase at the beginning of cycle 23 and at its maximum phase, respectively. The subtraction of the best fit allowed us to unveil the power spectrum of other sources of fluctuations. We notice that its amplitude increased significantly with the level of activity. A broad concentration of power was present between 0.15 d^-1 and 0.6 d^-1and it could be tentatively associated with the growth and decay of individual bipolar active regions on the solar disk, as suggested by the similarity of the time scales (cf. Harvey 1993). A curious and unexplained feature was the sharp peak at 0.94 d^-1 which was present unambiguously only close to the maximum of activity.

In the case of solar-like stars a similar type of variability is likely to be present. Its modelling is complicated because the physical processes responsible for the active region appearance and evolution are only partly known. However, an extrapolation based on the solar behaviour and magnetic flux evolution models may be used to predict stellar variability (Schrijver 2001). Our approach will provide an additional possibility to test these models by looking at the power spectrum of the stellar variability, after subtraction of the rotational modulation fit.

$\begin{figure} \par\includegraphics[width=7.6cm,clip]{3322.f14}\end{figure}$

Figure 14: a) The power spectrum of the relative variation of the TSI during the time interval 1996.0-1997.5, at a phase of minimum activity close to the beginning of solar cycle 23; b) The power spectrum of the same data sequence after subtracting our model best fit.

$\begin{figure} \par\includegraphics[width=7.6cm,clip]{3322.f15}\end{figure}$

Figure 15: a) The power spectrum of the relative variation of the TSI during the time interval 1999.5-2001.0, at a phase of maximum activity during the course of solar cycle 23; b) The power spectrum of the same data sequence after subtracting our model best fit.

6 Conclusions

The modulation of the solar bolometric flux by its active regions is affected both by their evolution and the solar rotation. The period of the modulation is comparable with the solar synodic period only during the phases of minimum activity, when long-lived faculae dominate the flux variations. During phases of more intense activity, the evolution of the large active regions on time scales shorter than one solar rotation produces a more complex variation and the real rotation period can not be retrieved. This warns us that in the case of active stars it may be difficult to measure the real rotation period if the lifetime of the active region pattern that dominates the variability is shorter than the rotation period. However, even in this case, it is possible to fit most of the variability by considering the rotational modulation due to a few discrete active regions, the configuration of which is changed at intervals comparable with their lifetime, and a uniformly distributed, time-dependent background. This simple model is capable to account for the variability of the TSI on time scales ranging from about one week up to the duration of the solar cycle. The amplitude of the residual variations obtained by subtracting the fit is at least one order of magnitude smaller than the observed TSI variations. This will allow us to increase significantly the efficiency of the method proposed by Defay et al. (2001) for the detection of transits of Earth-sized planets across solar-like stars. A comprehensive account of the advantages and drawbacks of our modelling approach for transit detection will be presented in a forthcoming paper (Barge et al., in preparation).

Moreover, an analysis of the residuals of the fit may allow us to study the sources of variability on time scales significantly shorter than the stellar rotation, such as the evolution of individual bipolar regions within large activity complexes. The variation of the total model area of the active regions can be linearly correlated with the variation of the observed area of the sunspot groups, supporting the possibility of detecting long-term changes of the activity level in stars, i.e., activity cycles. However, the derivation of active region characteristics is hampered by the simultaneous presence of bright and dark features on the stellar surface in a proportion which is not known a priori and by the limited information content of wide-band light curves. Only in the case of quite simple active region configurations, our technique is capable of retrieving their longitudes with an accuracy of the order of $15^{\circ}{-}20^{\circ}$ .

An improvement of our model may come from the simultaneous fitting of multi-band light curves that may place additional constraints on the contrast and position of surface inhomogeneities.

Contents

Modelling the rotational modulation of the Sun as a star

1 Introduction

2 The total solar irradiance data set and its variability

3 The model of the solar rotational modulation

4 Results

5 Discussion

6 Conclusions

References