© ESO, 2013

1. Introduction

The solar dynamo generates the magnetic field of the Sun and causes its 11-year activity cycle. One of its major ingredients is believed to be differential rotation (DR) of the surface. The interaction of rotation and convection produces turbulence in the convection zone, leading to non-uniform rotation (Kitchatinov 2005). On the Sun the angular velocity Ω decreases from the equator to the poles. Assuming an initial poloidal magnetic field with frozen field-lines, DR winds up the lines and transforms an initial poloidal field into a toroidal field (Ω-effect). The opposite effect, which transforms the toroidal field back into a poloidal one by twisting of the field-lines, is called α effect. In contrast to the Ω-effect, the α-effect is able to produce a toroidal field from a poloidal one, and vice versa. Furthermore, the strength of DR varies with spectral type. Barnes et al. (2005) found that ΔΩ strongly increases with effective temperature, supplying the power law $\mathrm{\Delta \Omega }\mathrm{~}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{8.92}}$. For temperatures above 6000 K this trend was confirmed by Reiners (2006). Collier Cameron (2007) combined results from Doppler imaging (DI) and the Fourier-transform method (see below), yielding the equation ΔΩ = 0.053(T_eff/5130)^8.6. This could be a hint that there are different dynamo mechanisms, but the final role of DR is yet to be understood.

The relation between rotation period and DR has been studied by several authors. Hall (1991) found that the relative horizontal shear α increases for longer rotation periods. Donahue et al. (1996) confirmed this trend, finding ΔP ~ ⟨P⟩^{1.3    ±    0.1}, independent of the stellar mass. Using the Fourier-transform method, Reiners & Schmitt (2003) also found that α increases with rotation period for F-G stars. This result has been confirmed by Ammler-von Eiff & Reiners (2012), who compiled previous results and new measurements for A-F stars. Barnes et al. (2005) found the relation ΔΩ ~ Ω^0.15, which shows the weak dependence of DR on the rotation-rate. We believe that the potential of high-precision photometry (e.g. provided by Kepler) has not been fully exhausted for measuring DR, especially for cooler stars.

The solar DR has been theoretically studied for a long time. Kitchatinov & Rüdiger (1999) computed DR models for late-type (G2 and K5) stars. They found that the relative shear α increases with rotation period. They also showed that α increases for cooler stars. Küker & Rüdiger (2005a) computed models for an F8 star and found a weak dependence of the total horizontal shear on rotation period, which was confirmed by subsequent studies for F, G, and K stars (Küker & Rüdiger 2005b), showing that the dependence of the horizontal shear on temperature is much stronger. This result holds for different main-sequence-star models, according to Küker & Rüdiger (2007). These authors found that above a temperature of 5500 K the strong temperature dependence of the horizontal shear (~${\mathit{T}}_{\mathrm{eff}}^{\mathrm{8.92}}$) from Barnes et al. (2005) fits the model data reasonably well, whereas below 5500 K the data lie far off the fit. Recent studies (Küker & Rüdiger 2011) have shown that the temperature dependence of ΔΩ cannot be represented by one single power law for the whole temperature range from 3800–6700 K. The weak dependence of ΔΩ on the rotation period was confirmed for different solar-mass models. Hotta & Yokoyama (2011) modeled DR of rapidly rotating solar-type stars and found that DR approaches the Taylor-Proudman state, meaning that ΔΩ/Ω decreases with angular velocity as long as the rotation-rate is faster than the solar value. Browning (2011) provided an explanation for the increase of DR for cooler stars. M dwarfs exhibit low luminosities, and therefore low convective velocities. Thus, they strongly depend on the rotation rate even at solar rotation-rates. Magnetic fields cannot quench DR effectively, so they exhibit strong solar-like DR.

DR can be measured in different ways. Ammler-von Eiff & Reiners (2012) summarized the results using the Fourier-transform method which analyzes the shape of Doppler-broadened spectral lines (Reiners & Schmitt 2002). DI tracks the migration of individual active regions over time to draw conclusions about the stellar rotation law. This method has successfully been used for instance by Donati & Collier Cameron (1997) and Collier Cameron et al. (2002). DR can also be detected by analyzing the star’s light curve, as done, for instance, by Hall (1991). An analytical spot model (Budding 1977; Dorren 1987) was fit to the light curve of ϵ Eridani in Croll et al. (2006), accounting for different spot periods. The same light curve has been analyzed by Fröhlich (2007) using the Markov chain Monte Carlo method to estimate the parameters in a Bayesian way. This approach has also been applied to the CoRoT-2 light curve in Fröhlich et al. (2009). Kipping (2012) presented a faster-operating version of the spot model from Budding (1977) and Dorren (1987), which accounts for DR and spot evolution in the parameter space.

Following another approach, Lanza et al. (1993) created light curves of spotted stars and detected different periods by taking the Fourier-transform. Walkowicz et al. (2013) fit an analytical spot model to synthetic light curves of spotted stars to see whether the model can break degeneracies in the light curve, especially accounting for the ability of determining the correct rotation periods, both with and without DR. In this paper we follow the same approach: we ran an extensive Monte Carlo simulation, producing 100 000 light curves in an attempt to cover a significant portion of the parameter space. All light curves were then analyzed using a Lomb-Scargle periodogram that incorporates a global sine fit. This is a fast and easy tool for deriving periods from photometric data. Thus, the fundamental idea is to determine how accurately this method can detect DR. This is necessary because with the advent of space telescopes such as CoRoT and Kepler it became possible to study stellar flux variations of thousands of stars at the level of milli-magnitude precision. Different effects, such as co-rotating spots on the surface, pulsations, spot evolution, or instrumental effects can introduce periodic variations in the light curves as well. Therefore, the accuracy of this method needs to be tested for spot-induced signals alone.

The paper is organized as follows: in Sect. 2 the routine MODSTAR is introduced, which synthesizes the light curves for the Monte Carlo simulation. Section 3 describes the period detection based on the Lomb-Scargle periodogram, with special focus on prewhitening, the selection process of the periods, and the sample properties. The results are found in Sect. 4, starting with the most significant period P1_out of the star in Sect. 4.1. In Sect. 4.2 the DR results are presented. In Sect. 4.3 we compare our method with previously published differentially rotating Kepler stars. The results are summarized in Sect. 5, briefly discussing the model parameters and analysis method.

2. Model description

Stellar rotation in the presence of active regions leads to photometric variability in the stars’ light curve. MODSTAR is our basic routine, which creates a model star to simulate the photometric signal of a rotating spotted star. The star is modeled as a sphere with a fixed resolution of the surface pixels and inclination of the rotation axis. The intensity I and projected area of each surface element depends on the value of μ, which is the cosine of the angle between the pixel’s surface normal and the line of sight. A quadratic limb-darkening law is used, $\mathit{I}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{a}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\mu }\mathrm{)}\mathrm{+}\mathit{b}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\mu }{\mathrm{)}}^{\mathrm{2}}\mathrm{)}\mathit{,}$(1)with I₀ being the intensity at the star’s center. In view of the Kepler mission we used the values a = 0.5287 and b = 0.2175 from Claret (2000), which relate to solar-like stars in the V-band.

Active regions can be placed on the surface. We used circular spots with the desired longitudes, latitudes, radii, and a fixed intensity-contrast. For the intensity-contrast between the spots and the quiet photosphere we used a value of 0.67, which is approximately the solar penumbra-to-photosphere contrast. For simplicity, the spot pixels all had the same contrast value, meaning that we neglected the umbra and penumbra structure. The stellar flux was integrated over the whole surface by summing the pixel intensities weighted by projected area. Since the star can be rotated, the flux was calculated at each rotation step, which produces a light curve. To describe the implemented rotation law, we quantify the amount of shear by $\mathit{\alpha }\mathrm{=}\frac{{\mathit{P}}_{\mathrm{pole}}\mathrm{-}{\mathit{P}}_{\mathrm{eq}}}{{\mathit{P}}_{\mathrm{pole}}}\mathit{,}$(2)with α = 0 supplying rigid body rotation. P_pole and P_eq are the rotation periods at the pole and the equator, respectively. α > 0 means that the equator rotates faster than the poles (solar-like DR), whereas α < 0 describes the opposite effect (anti-solar-like DR). For our simulations, we only considered solar-like DR since one cannot distingish between the two effects from the light curve alone. The rotation period of a spot centered at a certain latitude θ is given by a common solar-like DR law: ${\mathit{P}}_{\mathrm{spot}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{eq}}}{\mathrm{1}\mathrm{-}\mathit{\alpha }{\sin }^{\mathrm{2}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\mathrm{·}$(3)According to Eq. (3), a spot would be torn apart after some rotation cycles. To avoid this, the spots were fixed on the surface for all phases. In this way, we achieved long-lived spots that produce a stable beating pattern in the light curve. Spot evolution is not included in our model so far. The spots are allowed to overlap with no further contrast reduction.

2.1. Monte Carlo simulation

The Kepler mission provides light curves of all kinds of stellar activity, especially rotation-induced variability. In some stars DR has been detected (Frasca et al. 2011; Fröhlich et al. 2012), and many other light curves exhibit similar patterns. Inspired by this potpourri of active stars we wondered how accurately DR can be measured solely from photometry if we allow for different kinds of stellar properties and spot configurations.

We ran a Monte Carlo simulation, producing 100 000 light curves of spotted stars to account for a large fraction of possible realizations. The most important stellar parameters are the inclination i, the number of spots on the surface, and the amount of DR α. All parameters were uniformly distributed with sin(i) ∈  [0,1], the number of spots between 1 and 10, and α ∈  [0,1/3]. The inclination covered the whole parameter space from pole-on (i = 0°) to edge-on (i = 90°) view. The spot positions were chosen at random, and the spot radii were between 2° and 21° (Tables 1 and 2). The number of spots was limited to ten because the spot radii can be very large. These two limits prevent the star from being completely covered with spots, which would result in a darker and therefore cooler star. The smaller the spot radii, the more pixels are needed in our model to resolve individual surface elements. We fixed the minimum spot radius to 2° to limit the computational effort. α covers the parameter space from rigid body rotation (α = 0) up to α = 1/3, including the solar value of α_⊙ = 0.20. According to Eq. (3), the period is a function of latitude with P = 1 cycle at the equator to P = 1.5 cycles at the poles. A higher α value, e.g., α = 0.5 would result in rotation periods of P = 1 cycle at the equator to P = 2 cycles at the poles. This would lead to problems in distingishing between harmonics and true rotation-rates. In observations values of α > 0.5 have been found (Ammler-von Eiff & Reiners 2012). DI usually yields lower values of about α ≲ 0.01, see Table1 in Barnes et al. (2005), and references therein. Our method is sufficient to produce light curves of differentially rotating stars for a wide range of α, but becomes problematic for either very high or very low α values. This problem is discussed in Sect. 5.

The spot-to-photosphere contrast was fixed because there is a degeneracy between spot size and contrast. Choosing the contrast as free parameter in the simulation would just change the depth of the spot signature in the light curve, but does not affect the period. The limb-darkening coefficients are fixed for all light curves. Each light curve consists of 300 data points that cover ten rotation periods P_eq to see how the light curves evolve in time.

In the following we consider different noise levels in the light curves: the noise-free case, 100 ppm, 1000 ppm, and 10 000 ppm Poisson noise, which is added to the light curves. A minimum noise level of 100 ppm was chosen because it is lower than the depth produced by the smallest spot (2°), which we found to be approximately 250 ppm. The same argument applies to the highest noise level, which is lower than the depth of the largest spot (21°) and is approximately 48 000 ppm. An example of a 1000 ppm model light curve is shown in Fig. 1. This light curve looks similar to what we see in Kepler data, therefore we are optimistic that the parameter selection is sufficient for our main purposes – the production of a variety of light curves with periodic variability and the detection of DR. We aimed to keep the model simple and were interested to see whether this can reproduce real light curves. A larger parameter space can be tackled if we find that our model cannot reproduce the data. A similar argument applies to spot evolution, which so far we have excluded from our approach.

3. Period determination

We detect DR by measuring different periods in a light curve. Because we are facing a large sample, a fast and reliable frequency-analysis tool is needed. We chose the Lomb-Scargle periodogram, which is widely used in time-series analysis. Calculating the periodogram of a single light curve takes only one second. Although it is a purely mathematical tool, the program is sufficient to find different periods in the data. Fitting an analytical spot model also supplies rotation periods and several other stellar parameters, but severely slows down the analysis process.

Fig. 1 Upper panel: model light curve with 1000 ppm Poisson noise added and the best global sine fit overplotted (green). The spot periods of the model and those returned by our method are given in Table 3. Lower panel: residuals of the best fit subtracted from the data; no periodicity visible anymore.

3.1. Lomb-Scargle periodogram

The generalized Lomb-Scargle periodogram (Zechmeister & Kürster 2009)¹ is a powerful spectral-analysis tool for unevenly sampled data. It fits the data using a series of sines and cosines. The frequency grid used for the fit is sampled equidistantly. Its range has an upper limit because of the Nyquist frequency. The lower limit is given by the inverse product of the time span and a desired oversampling factor to achieve a proper frequency resolution. We used a factor of 10, i.e., a minimum frequency of 0.01/cycle. Depending on the goodness of the fit, one obtains peaks with different powers – the better the fit, the higher the peak in the periodogram. The periodogram is normalized to unity. The period P (or frequency f = 1/P) associated to the highest peak is the most dominant one in the data. In some cases an alias of P/2, P/3, etc. may also produce a peak with high power. One reason for an alias period can be the presence of two active longitudes separated by approximately 180° from each other. Another one can be the wrong shape of a sine wave to fit the spot signature. A single spot does not produce a sinusoidal shape except for a pole-on view. In the frequency domain aliases are equidistant (2f, 3f, etc.) and can easily be detected by eye because the peak height usually decreases towards higher harmonics. We excised most alias periods as described in Sect. 3.3.

3.2. Prewhitening

In most cases our model stars are covered by several spots, adding their signals to one light curve. Thus, we are facing the challenge to fit a mixture of period signatures. These periods are extracted from the light curve in a successive way called prewhitening. First, we adopted the period associated to the highest peak in the periodogram and fit a sine wave to the light curve. The initial sine function was subtracted from the data and another periodogram was taken from the residuals. Again, we fit a sine function and subtracted it from the data. This prewhitening process can be repeated as often as desired, i.e., until there is no periodicity present anymore. On the one hand, a high number of prewhitening steps is crucial for the detection of several periods, but on the other hand, prewhitening is computationally intensive and one has to be careful to select the correct periods afterwards (Sect. 3.3), which becomes more difficult with a larger set of periods. Since the stars are covered by ten spots or fewer, we repeat this procedure ten times for each light curve. Finally, all ten periods detected during prewhitening were used as input parameters for a global sine fit, which is the sum of ten sine functions with different periods, amplitudes, phases, and one total offset. The result of this last step is a best-fit set of parameters found through χ² minimization. A model light curve with 1000 ppm noise and the fit obtained through this procedure are displayed in Fig. 1. In the first column of Table 3 the actually contained spot periods are shown; this particular example has only three spots. The second column contains all detected periods, which are arranged as follows: the first period has the highest power found in the prewhitening process, the second one belongs to the second highest power, and so on. We see that in this case the third and fifth period are harmonics of the first and second one, respectively, the fourth, and the last five periods just fit the remaining noise. The fit in Fig. 1 agrees well with the light curve, and the residuals carry no more periodicity in the domain of one or more cycles.

3.3. Period selection

The fitting procedure described above returns ten periods for each light curve. In this section we assign a physical meaning to the most significant period (compare Table 3), and show how to detect additional periods adjacent to it as evidence for DR. Even though we know that in our model the spot periods range from 1−1.5 cycles (Eq. (3)), we did not constrain our algorithm to this range since in reality we do not know the correct period either.

The first sine period is the most significant one in the data. In several cases this period is equal to the first harmonic (P/2) of the true spot rotation-rate due to certain spot configurations, e.g., two spots located on opposite sides of the star. In these cases the second sine period is most likely the correct period P. To minimize the number of alias periods we compared the double of the first period to the second sine period. If these two differ by less than 5%, the second sine period was chosen. The finally selected period is our primary period P1_out. In Sect. 4.1 we show that this period yields the best estimate of the stellar rotation-rate. In most cases the light curves result from several spots on the surface, therefore we need to compare P1_out with the spot periods according to Eq. (3). This set is called the input periods since they belong to spots inherent in the light curves. Of these, we took the one closest to P1_out and call it P1_in.

Based on the period P1_out, we searched for a second period, which we call P2 for the moment. To attribute this period to a second spot, one has to achieve three tasks: 1) find a second period that is no harmonic of P1_out; 2) try to exclude as few spot periods as possible; and 3) try to dismiss all period artifacts that come from fitting a sine wave to the light curve and not a spot model. Therefore, P2 should hold the relations $\mathrm{0.01}\mathrm{\le }\frac{\mathrm{|}\mathit{P}{\mathrm{1}}_{\mathrm{out}}\mathrm{-}\mathit{P}\mathrm{2}\mathrm{|}}{\mathit{P}{\mathrm{1}}_{\mathrm{out}}}\mathrm{\le }{\mathit{\alpha }}_{\max }\mathit{,}$(4)with α_max = 1/3. The value α_max = 1/3 corresponds to the maximum α in our model. As mentioned above, a higher value near α_max ≈ 0.5 would yield ambiguous results for the DR. Imagine the two cases: 1) P1_out = 1,P2 = 0.5; and 2) P1_out = 1,P2 = 1.5. In both cases a second period would be selected, but would come from completely different origins. In case 1) the first harmonic of P = 1 would be misinterpreted as DR, whereas case 2) results from real spot configurations. We chose the value α_max = 1/3 because it excludes harmonics and covers a wide range for a second spot period. For example, in the extreme case of P1_out = 1, the harmonic at P2 = 1/2 is excluded, but we are unable to find a spot period longer than P = 1.33 although there might be spots with longer periods. The lower limit in relation (4) accounts for the fixed frequency resolution in the Lomb-Scargle periodograms (Sect. 3.1). If two spot periods differ by less than 1%, they cannot be resolved. Again, we might miss some spot periods lying closer together than 1% with our method.

When one or more periods were found that fit both criteria (compare Table 3), we chose that P2 associated to the period with the lowest row-index in the table of remaining² sine periods. This period is called P2_out. In 88 330 cases of our models we found a P2_out that fulfills these criteria. Again, P2_out was compared with all input periods, and the period closest to P2_out was called P2_in. When the closest input period picked was again P1_in, P2_out was discarded. Finally, we were left with 64 172 stars with two periods that belong to two different spots.

For our example light curve in Fig. 1, the periods P1_out = 1.016,P2_out = 1.184,P1_in = 1.018, and P2_in = 1.186 were selected from Table 3, resulting in α_in = 0.14 and α_out = 0.14. Although our method is able to recover the correct α_in value, the total equator-to-pole shear equals α = 0.30 in this case (compare Eq. (3)), and thus is underestimated by more than 50%. For this specific spot configuration it is impossible to obtain the correct α value since the highest spot latitude equals θ = 46.3°, generating the longest period. This is a general problem of DR measurements from photometric observations due to the initial spot configuration on the surface. The measured shear will always yield a lower value than the total one.

3.4. Sample properties

This section is meant to be a consistency check of our model and the selection algorithm. The stellar parameters (Table 1) of two mutually exclusive samples, S2 and S1, are compared. The S2 sample consists of all light curves with two detected periods (64.2%) that come from two distinct spots, whereas S1 contains all cases where only one spot period could be associated (32.0%). Because of a combination of a low number of spots, certain spot latitudes, and a low inclination in 3.8% of all cases no spot was visible.

The goal is to point out the stellar models where DR can likely be found, compared with the cases where detecting DR is challenging or even impossible with our method. For example, one would expect that it is easier to detect DR in a highly spotted star than in a star covered by only two close-in spots because it will probably be hard to resolve individual periods in the latter case.

In Fig. 2 we compare the inclination, number of spots, spot radii, and differential rotation α of the two samples. The S2 sample is shown in the left and the S1 sample in the right column. The colors correspond to the different noise levels: noise-free (green), 100 ppm (yellow), 1000 ppm (orange), and 10 000 ppm (red) Poisson noise.

Fig. 2 Comparison of basic stellar parameters for the two samples S2 and S1 for different noise cases: noise-free (green), 100 ppm (yellow), 1000 ppm (orange), and 10 000 ppm (red) Poisson noise. Left panel: stars with two spot periods (S2). Right panel: stars with only one spot period (S1).

Due to the small difference of detections for the noise-free (green) and the 100 ppm noise case (yellow), there is basically no difference visible in the histograms. All trends have similar shape and become more distinct at higher noise levels.

In the first row the inclination for both samples is plotted. We recall that the sin(i) of the whole sample has a flat distribution. Starting from the edge-on view (i = 90°), we see a continuous decrease (increase) of detections in the S2 (S1) sample. Around inclinations lower than i = 10°, the number of detections in S2 decreases significantly. The opposite effect applies to S1.

In the second row we show distributions of the actual number of spots of the models. Only in a very few cases the models in S2 can be attributed to only two spots on the surface. In the majority of all cases the light curve is composed of the signature from more than five spots! Clearly, the actual number of spots decreases in the S1 sample.

The third row shows the distribution of the spot radii. In the left panel one clearly sees that the radii of the two spots increases for higher values because it is more likely to detect more than one period if the star has large spots. In the right panel we plot the radii of all visible spots except for the one associated to the one period found. These residual or unresolved spots show a shallow decrease in radii.

Finally, the last row shows histograms of the α value. Both distributions are basically flat except for the 10 000 ppm noise case and for low values of α. The distribution of S2 decreases while that of S1 increases for lower α values.

All histograms show consistent results for each sample, supporting the underlying model. With a focus on the S2 sample the selection process convincingly picks mostly models for which the detectability of DR is expected. The derived periods and accuracy of our tool are discussed in the following section.

4. Results

In this section we compare the outcome of our analysis with the periods from our model. First, we present the basic results for the most dominant period P1_out. In Sect. 4.2 the detection of DR is discussed considering different noise levels. Finally, we compare our method with three Kepler stars with confirmed DR (Sect. 4.3).

Fig. 3 Main plot: comparison of the exactly known input periods P1in and output periods P1out for the noise-free case. The most significant period P1out can be recovered fairly well. P1out is on average 2.4% shorter than the actually contained period P1in. Small plots: the distribution of P1out is shown with increasing noise from top to bottom. For higher noise levels the fraction of periods lower than 0.5 cycles increases because the algorithm interprets the noise as short periods (compare Table 4).

4.1. Rotation periods P1_out

96.2% of all light curves have one detected periodicity P1_out that is the most significant one in the data. In Fig. 3 we compare the periods P1_out with the input periods P1_in for different noise levels to see how well our selection process works. The distribution of input periods P1_in is shown in gray, the shaded distribution shows the output periods P1_out. Because we only considered solar-like DR, the P1_in cannot be shorter than one cycle. The shaded gray area shows that the two histograms match quite well. For the noise-free case the number of wrong detections is negligible, with 3.3% of all periods shorter than 0.9 cycles and 2.9% longer than 1.6 cycles. This is no longer true for higher noise levels. The region from 0−0.5 cycles becomes populated, as shown in the three small plots in Fig. 3. We compare the weighted means ⟨P1_in⟩ and ⟨P1_out⟩ in the above range in Table 4. In the noise-free case P1_out is on average 2.4% shorter than the actually contained period P1_in. For higher noise levels ⟨P1_out⟩ decreases because the algorithm interprets the noise as short periods. In the noise-free case the wrong detections are due to a low stellar inclination close to the pole-on view, or in some cases due to higher harmonics. P1_out is shorter than one cycle in some cases because a sine function has the wrong shape to fit spot signatures in a light curve. For example, a detected period of P1_out = 0.98 cycles will be considered as a valid rotation period although it is not possible for the spots to rotate this fast in our model. This will not be noticed in real data because we do not have information on the real rotation of a star. Only a small fraction of harmonics remained (less than 0.5%) around a period of 0.5 cycles after identification and correction (Sect. 3.3).

4.2. Differential rotation

In this section we show in which situations DR can successfully be detected and which cases lead to wrong interpretations. In our model the detection of a second period P2_out adjacent to the primary one and associated to a second spot period P2_in is considered as evidence for DR. To investigate where this selection process is acceptable we considered the two pairs (P1_out,P2_out) and their associated periods (P1_in,P2_in). We estimated the amount of DR by sorting these pairs and computing their α_in,out value, (Table 3): ${\mathit{\alpha }}_{\mathrm{in}\mathit{,}\mathrm{out}}\mathrm{=}{\frac{\mathit{P}\mathrm{1}\mathrm{-}\mathit{P}\mathrm{2}}{\mathit{P}\mathrm{1}}|}_{\mathrm{in}\mathit{,}\mathrm{out}}\mathit{,P}\mathrm{1}\mathit{>}\mathit{P}\mathrm{2.}$(5)The α_in value is always lower than the inherent α value of each light curve since we can only measure the rotational shear at two defined latitudes. If α_in is calculated from the spots with the largest separation in latitude on a certain hemisphere, then α_in is the maximum shear that can be detected by our method. The distribution of α_out − α_in gives a statistical measure of the robustness of our period-selection process. In Fig. 4 we show the distribution of the differences between α_out and α_in in our set of models for different noise levels. All distributions exhibit an asymmetric shape for too high α_out values. A reason why α_out can be too high is given at the end of this section.

Fig. 4 Histograms of αout − αin for all noise levels. From the upper left to the lower right panel the noise increases. For the noise-free case (upper left panel) in 64 172 stars (at least) two periods are detected. The distribution is centered at ⟨αout − αin⟩ = 0.032 with a width of σ(αout − αin) = 0.058. With higher noise levels the total number of findings decreases and the error increases (Table 5).

The total number of light curves with two (or more) detected periods, the weighted mean, and the error of the distributions are given in Table 5. In general, the total number of stars with two detected periods decreases with increasing noise. All differences between the noise-free and the 100-ppm case are marginal. The 1000-ppm case has slightly fewer detections, and the most significant decrease occurs for the 10 000-ppm case (compare Fig. 2). For all cases the weighted mean and the error increase with noise. For the noise-free distribution (upper left) we found a weighted mean ⟨α_out − α_in⟩ = 0.032 and a width of σ(α_out − α_in) = 0.058. The error is given in absolute units of α_out − α_in regardless of the true shear of α. We also considered the relative errors | α_out − α_in | /α for α > 0.05. Each distribution is proportional to constant/α, i.e., σ(α_out − α_in) does not scale with α. For each noise level the number of stars decreases with increasing relative error. We only considered α > 0.05 to avoid large errors. The statistical error resulting from the sample of the light curves used is negligible. The above calculations were made for two additional sets of light curves (each 100 000 in total). We found that ⟨α_out − α_in⟩ and σ(α_out − α_in) are almost equal for each set and that the largest statistical uncertainty in the number of stars with two detected periods is about 0.3%.

Fig. 5 Distribution of αout − α for the noise-free case. It is centered at ⟨αout − α⟩ =  −0.056 with a width of σ(α − αout) = 0.098. Comparing this plot with Fig. 4 shows that the total amount of DR is underestimated by 8.8% because αin ≤ α by definition. The detections where α < 0.06 are overplotted in red. Clearly, in these cases the detection of DR is difficult. In most cases, αout > α which demonstrates the limits of our method.

We have shown that DR can be measured with high accuracy for a broad noise range. In Fig. 5 we compare α_out with the total equator-to-pole shear α of the star. In the noise-free case the distribution has a weighted mean ⟨α_out − α⟩ =  −0.056 and an error σ(α_out − α) = 0.098. Considering the mean values yields $\begin{array}{ccc}\mathrm{0.032}& \mathrm{=}& ⟨{\mathit{\alpha }}_{\mathrm{out}}\mathrm{-}{\mathit{\alpha }}_{\mathrm{in}}⟩\mathrm{\ge }⟨{\mathit{\alpha }}_{\mathrm{out}}\mathrm{-}\mathit{\alpha }⟩\mathrm{=}\mathrm{-}\mathrm{0.056}\\ & \mathrm{\Rightarrow }& ⟨{\mathit{\alpha }}_{\mathrm{out}}\mathrm{-}\mathit{\alpha }⟩\mathrm{\le }\mathrm{0.088}\mathit{,}\end{array}$since α_in ≤ α by definition. This means that the total amount of DR is underestimated by 8.8% by our method. Models with low α values exhibit spot periods very close to each other. Therefore, we tested whether these are prone to misidentification of DR. In Fig. 5 the models with an equator-to-pole shear α < 0.06 are overplotted in red. We found that α_out yields wrong values higher than α itself. This results from the difficulty of resolving two distinct peaks in the Lomb-Scargle periodogram. Since the peak width is proportional to the inverse time span of the light curve, we are unable to resolve two periods within ten cycles. One broadened peak appears, which results in a mean period. In these cases the initial sine wave does not have a proper shape to fit the mixture of spot signatures, therefore there remain artifacts that are corrected for by fitting more sine waves. If one of the residual periods fulfills the selection criteria, the algorithm selects it as P2_out, yielding too high α_out values. This behavior partly applies to the distributions in Fig. 4, especially for the highest-noise case. A second spot period cannot be resolved properly and is lost in the noise, yielding an asymmetric distribution.

4.3. Comparison with Kepler data

We have seen that our method performs well for simulated data. To test its reliability, we applied it to three previously studied Kepler stars (Frasca et al. 2011; Fröhlich et al. 2012) for which DR is the favored explanation for the light curve variability. In both papers synthetic light curves from an analytical spot model (Dorren 1987) were fit to the data, and the parameters were estimated in a Bayesian way using Markov chain Monte Carlo methods. The stars KIC 8429280 and KIC 7985370 were fit with seven spots, and the star KIC 7765135 was fit with nine spots. In the model of Frasca et al. (2011), Fröhlich et al. (2012) each spot has a certain rotation period. Additional model parameters are the equator-to-pole differential rotation dΩ, and Fröhlich et al. (2012) also used the equatorial period P_eq. We used our method to determine P1_out and P2_out to compare them with the reported results. We call the spot periods closest to our findings P1_closest and P2_closest. From these periods we calculated α_out and α_closest. The last row in Table 6 contains the calculated equator-to-pole shear dΩ/Ω_eq of the models that equal α in our simulations. Since KIC 8429280 has no parameter for the equatorial period, we used the period from the spot closest to the equator (θ =  −0.5°). For all three stars we determined consistent results for the horizontal shear (Table 6).

Some spots fitted in the model have very close or even identical periods. According to relation (4) we are not able to detect periods with a smaller difference than 1%. Even without this restriction one would need a very high frequency resolution in the Lomb-Scargle periodograms to see these shallow differences. Therefore, the most significant periods in the data have probably been successfully detected by our method because the returned α values are very similar. Recently, Roettenbacher et al. (2013) determined a surface rotation period of 3.47 days for the Kepler target KIC 5110407 and found evidence for DR and spot evolution using light curve inversion. We found P1_out = 3.61 and P2_out = 3.42 days, yielding α_out = 0.052. Roettenbacher and collaborators DR coefficient k (equal to α in our case) strongly depends on the assumed inclination angle and varies between k = 0.024 − 0.118. Using i = 45° these authors found k = 0.053, which is very close to our result. Hence, our method successfully returned a consistent DR coefficient using a faster and much easier method than light curve inversion.

5. Summary

We have run a comprehensive Monte Carlo simulation of differentially rotating, spotted stars covering a significant fraction of the parameter space. The resulting light curves were analyzed with the Lomb-Scargle periodogram in a prewhitening approach. The returned periods were used for a global sine fit to the data. The main goal was to determine under which stellar conditions and how accurately DR can be detected.

The most significant periodicity P1_out was detected in 96.2% of all light curves. A second period close to P1_out was attributed to a second spot found in 64.2% of all cases. The latitudinal shear of the two spots associated to the above periods has been calculated. We found that the shear of the two spots α_in was on average 3.2% lower than α_out. Furthermore, comparing α_out with the total equator-to-pole shear, we found that α was underestimated by 8.8%. In particular, we found that the detection of DR is challenging for stars with an equator-to-pole shear lower than 6%, which usually yields large errors.

In our model each light curve is composed by a fixed number of spots rotating at defined latitudes. On the Sun the situation is different. Spots vanish or are created while rotation takes place. Their preferred latitudes of occurrence are around θ = 30°, and during the 11-year activity cycle they migrate toward the equator. So far, our program does not account for spot lifetimes or meridional drifts.

The least known model parameter is the number of spots on the surface and their associated size. Alternative to our approach, one could also use more spots with smaller radii, keeping the spot filling-factor constant, i.e., the fraction of the surface covered with spots. But this also increases the number of available spot periods. Comparing them with the results of our analysis bears the risk of choosing a spot that is not responsible for the detected signal.

On average, the stars in our simulation exhibit five spots with radii of 12°, which is due to the chosen parameter range and its uniform distribution. It would be more realistic to link the spot radii with the activity level of the stars. It is well known that younger stars are more active than the Sun, and exhibit on average larger active regions for instance, whereas on the Sun small spots at preferred latitudes are observed. Our Monte Carlo simulation does not assign activity levels to the stars, which means that DR studies of inactive stars such as the Sun are not fully covered by our model so far.

Our method is limited to a relative DR of α < 0.5. Allowing for α ≥ 0.5 in the model would lead to confusion between real spot periods and aliases of faster-rotating spots with P1 ≤ 1/2P2. This is a general problem of DR determination from photometric data. Our lower limit in relation (4) accounts for the frequency resolution in the periodograms. This limit is only relevant for a reasonable computation time and could in principle be discarded, unlike the upper limit. Observations of DR cover a wide range of α values. The Doppler-imaging technique is particularly sensitive to low equator-to-pole shear (α ≲ 0.01) limited by the spot lifetimes, although there are measurements (Donati et al. 2003) that yield α ≈ 0.05 for LQ Hya, and new measurements (Marsden et al. 2011) that found values between 0.005 ≲ α ≲ 0.14. The Fourier-transform method (e.g. Reiners & Schmitt 2003) is sensitive to α > 0.1 and has been used to determine surface shears as strong as 50% for some A-F stars.

We conclude that our results confirm the possibility to reliably detect DR from photometric data using a fast tool with relatively simple mathematical assumptions. In real data it will be difficult to correctly interpret the results because of additional effects such as spot evolution, pulsations, instrumental effects, and combinations of them. This does not diminish the power of our analytic tool, but its applicability to real data needs to be tested on a complex sample of high-quality light curves. We will apply our method to the exquisite sample of the Kepler satellite in a forthcoming publication.

Acknowledgments

T.R. acknowledges support from the DFG Graduiertenkolleg 1351 Extrasolar Planets and their Host Stars. A.R. acknowledges financial support from the DFG as a Heisenberg Professor under RE 1664/9-1.

References

All Tables

All Figures

	Fig. 1 Upper panel: model light curve with 1000 ppm Poisson noise added and the best global sine fit overplotted (green). The spot periods of the model and those returned by our method are given in Table 3. Lower panel: residuals of the best fit subtracted from the data; no periodicity visible anymore.
In the text

	Fig. 2 Comparison of basic stellar parameters for the two samples S2 and S1 for different noise cases: noise-free (green), 100 ppm (yellow), 1000 ppm (orange), and 10 000 ppm (red) Poisson noise. Left panel: stars with two spot periods (S2). Right panel: stars with only one spot period (S1).
In the text

Fig. 3 Main plot: comparison of the exactly known input periods P1in and output periods P1out for the noise-free case. The most significant period P1out can be recovered fairly well. P1out is on average 2.4% shorter than the actually contained period P1in. Small plots: the distribution of P1out is shown with increasing noise from top to bottom. For higher noise levels the fraction of periods lower than 0.5 cycles increases because the algorithm interprets the noise as short periods (compare Table 4).

In the text

Fig. 4 Histograms of αout − αin for all noise levels. From the upper left to the lower right panel the noise increases. For the noise-free case (upper left panel) in 64 172 stars (at least) two periods are detected. The distribution is centered at ⟨αout − αin⟩ = 0.032 with a width of σ(αout − αin) = 0.058. With higher noise levels the total number of findings decreases and the error increases (Table 5).

In the text

Fig. 5 Distribution of αout − α for the noise-free case. It is centered at ⟨αout − α⟩ =  −0.056 with a width of σ(α − αout) = 0.098. Comparing this plot with Fig. 4 shows that the total amount of DR is underestimated by 8.8% because αin ≤ α by definition. The detections where α < 0.06 are overplotted in red. Clearly, in these cases the detection of DR is difficult. In most cases, αout > α which demonstrates the limits of our method.

In the text