Issue 
A&A
Volume 557, September 2013



Article Number  A11  
Number of page(s)  10  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201321161  
Published online  13 August 2013 
Fast and reliable method for measuring stellar differential rotation from photometric data
Institut für Astrophysik, Universität Göttingen, FriedrichHundPlatz 1, 37077 Göttingen, Germany
email: reinhold@astro.physik.unigoettingen.de
Received: 24 January 2013
Accepted: 10 June 2013
Context. Corotating spots at different latitudes on the stellar surface generate periodic photometric variability and can be useful proxies for detecting differential rotation (DR). This is a major ingredient of the solar dynamo, but observations of stellar DR are very sparse. Because the Kepler space telescope steadily collects more data, we are interested in detecting DR using photometric information of a star.
Aims. The main goal of this paper is to develop a fast method for determining stellar DR from photometric data.
Methods. We ran an extensive Monte Carlo simulation of differentially rotating spotted stars with very different properties to investigate the detectability of DR. For different noise levels the resulting light curves were prewhitened using LombScargle periodograms to derive parameters for a global sine fit to detect periodicities.
Results. We show under which conditions DR can successfully be detected from photometric data, and in which cases the light curve provides insufficient or even misleading information on the stellar rotation law. In our simulations, the most significant period P1_{out} is on average 2.4% shorter than the actual spot rotationrate. This period was detected in 96.2% of all light curves. The signature of DR is a second period close to P1_{out} in our model. For the noisefree case, we found such a period in 64.2% of all stars. Calculating the measured latitudinal shear of two distinct spots α_{out}, and comparing this with the known original spot rotationrates shows that the real value is on average 3.2% lower. Comparing the total equatortopole shear α to α_{out}, we find that α is underestimated by 8.8%, especially the detection of DR for stars with α < 6% is challenging. Finally, we applied our method to four differentially rotating Kepler stars and found close agreement with results from detailed modeling.
Conclusions. The method we developed is capable of measuring stellar rotation periods and detecting DR with relatively high accuracy and is suitable for large data sets. We will apply our analysis to more Kepler data in a forthcoming paper.
Key words: stars: rotation / starspots / methods: analytical / methods: statistical / techniques: photometric
© ESO, 2013
1. Introduction
The solar dynamo generates the magnetic field of the Sun and causes its 11year 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 nonuniform rotation (Kitchatinov 2005). On the Sun the angular velocity Ω decreases from the equator to the poles. Assuming an initial poloidal magnetic field with frozen fieldlines, 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 fieldlines, 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 . For temperatures above 6000 K this trend was confirmed by Reiners (2006). Collier Cameron (2007) combined results from Doppler imaging (DI) and the Fouriertransform 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 Fouriertransform method, Reiners & Schmitt (2003) also found that α increases with rotation period for FG stars. This result has been confirmed by Ammlervon Eiff & Reiners (2012), who compiled previous results and new measurements for AF stars. Barnes et al. (2005) found the relation ΔΩ ~ Ω^{0.15}, which shows the weak dependence of DR on the rotationrate. We believe that the potential of highprecision 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 latetype (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 mainsequencestar 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 (~) 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 solarmass models. Hotta & Yokoyama (2011) modeled DR of rapidly rotating solartype stars and found that DR approaches the TaylorProudman state, meaning that ΔΩ/Ω decreases with angular velocity as long as the rotationrate 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 rotationrates. Magnetic fields cannot quench DR effectively, so they exhibit strong solarlike DR.
DR can be measured in different ways. Ammlervon Eiff & Reiners (2012) summarized the results using the Fouriertransform method which analyzes the shape of Dopplerbroadened 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 CoRoT2 light curve in Fröhlich et al. (2009). Kipping (2012) presented a fasteroperating 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 Fouriertransform. 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 LombScargle 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 millimagnitude precision. Different effects, such as corotating 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 spotinduced 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 LombScargle 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 limbdarkening law is used, (1)with I_{0} 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 solarlike stars in the Vband.
Active regions can be placed on the surface. We used circular spots with the desired longitudes, latitudes, radii, and a fixed intensitycontrast. For the intensitycontrast between the spots and the quiet photosphere we used a value of 0.67, which is approximately the solar penumbratophotosphere 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 (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 (solarlike DR), whereas α < 0 describes the opposite effect (antisolarlike DR). For our simulations, we only considered solarlike 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 solarlike DR law: (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 longlived 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 rotationinduced 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.
Stellar simulation parameters.
Spot simulation parameters.
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 poleon (i = 0°) to edgeon (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 rotationrates. In observations values of α > 0.5 have been found (Ammlervon 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 spottophotosphere 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 limbdarkening 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 noisefree 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 frequencyanalysis tool is needed. We chose the LombScargle periodogram, which is widely used in timeseries 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. LombScargle periodogram
The generalized LombScargle periodogram (Zechmeister & Kürster 2009)^{1} is a powerful spectralanalysis 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 poleon 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 bestfit set of parameters found through χ^{2} 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 rotationrate 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 rotationrate. 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 (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 LombScargle 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 rowindex in the table of remaining^{2} 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 equatortopole 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 closein 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: noisefree (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: noisefree (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 noisefree (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 edgeon 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 P1_{in} and output periods P1_{out} for the noisefree case. The most significant period P1_{out} can be recovered fairly well. P1_{out} is on average 2.4% shorter than the actually contained period P1_{in}. Small plots: the distribution of P1_{out} 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 solarlike DR, the P1_{in} cannot be shorter than one cycle. The shaded gray area shows that the two histograms match quite well. For the noisefree 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 noisefree 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 noisefree case the wrong detections are due to a low stellar inclination close to the poleon 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).
Weighted means ⟨P1_{in}⟩ and ⟨P1_{out}⟩ and their associated errors σ(P1_{in}) and σ(P1_{out}) for both input and output periods for each noise level.
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): (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 periodselection 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 noisefree 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). 
Weighted mean ⟨α_{out} − α_{in}⟩ and error σ(α_{out} − α_{in}) of α_{out} − α_{in} for all stars with two detected periods for different noise levels.
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 noisefree and the 100ppm case are marginal. The 1000ppm case has slightly fewer detections, and the most significant decrease occurs for the 10 000ppm case (compare Fig. 2). For all cases the weighted mean and the error increase with noise. For the noisefree 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 noisefree 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 equatortopole shear α of the star. In the noisefree case the distribution has a weighted mean ⟨α_{out} − α⟩ = −0.056 and an error σ(α_{out} − α) = 0.098. Considering the mean values yields 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 equatortopole 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 LombScargle 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 highestnoise 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 equatortopole 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 equatortopole 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 LombScargle 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.
Comparison of periods returned by our method with spot periods from previously studied Kepler stars.
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 LombScargle 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 equatortopole shear, we found that α was underestimated by 8.8%. In particular, we found that the detection of DR is challenging for stars with an equatortopole 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 11year 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 fillingfactor 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 fasterrotating 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 Dopplerimaging technique is particularly sensitive to low equatortopole 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 Fouriertransform 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 AF 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 highquality light curves. We will apply our method to the exquisite sample of the Kepler satellite in a forthcoming publication.
For different periodogram codes see http://www.astro.physik.unigoettingen.de/~zechmeister/gls.php
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/91.
References
 Ammlervon Eiff, M., & Reiners, A. 2012, A&A, 542, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barnes, J. R., Collier Cameron, A., Donati, J.F., et al. 2005, MNRAS, 357, L1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Browning, M. K. 2011, in IAU Symp. 271, eds. N. H. Brummell, A. S. Brun, M. S. Miesch, & Y. Ponty, 69 [Google Scholar]
 Budding, E. 1977, Ap&SS, 48, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Claret, A. 2000, A&A, 363, 1081 [NASA ADS] [Google Scholar]
 Collier Cameron, A. 2007, Astron. Nachr., 328, 1030 [NASA ADS] [CrossRef] [Google Scholar]
 Collier Cameron, A., Donati, J.F., & Semel, M. 2002, MNRAS, 330, 699 [NASA ADS] [CrossRef] [Google Scholar]
 Croll, B., Walker, G. A. H., Kuschnig, R., et al. 2006, ApJ, 648, 607 [NASA ADS] [CrossRef] [Google Scholar]
 Donahue, R. A., Saar, S. H., & Baliunas, S. L. 1996, ApJ, 466, 384 [NASA ADS] [CrossRef] [Google Scholar]
 Donati, J.F., & Collier Cameron, A. 1997, MNRAS, 291, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Donati, J.F., Collier Cameron, A., & Petit, P. 2003, MNRAS, 345, 1187 [NASA ADS] [CrossRef] [Google Scholar]
 Dorren, J. D. 1987, ApJ, 320, 756 [NASA ADS] [CrossRef] [Google Scholar]
 Frasca, A., Fröhlich, H.E., Bonanno, A., et al. 2011, A&A, 532, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fröhlich, H.E. 2007, Astron. Nachr., 328, 1037 [NASA ADS] [CrossRef] [Google Scholar]
 Fröhlich, H.E., Küker, M., Hatzes, A. P., & Strassmeier, K. G. 2009, A&A, 506, 263 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fröhlich, H.E., Frasca, A., Catanzaro, G., et al. 2012, A&A, 543, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hall, D. S. 1991, in IAU Colloq. 130, The Sun and Cool Stars. Activity, Magnetism, Dynamos, eds. I. Tuominen, D. Moss, & G. Rüdiger (Berlin: Springer Verlag), Lect. Notes Phys., 380, 353 [Google Scholar]
 Hotta, H., & Yokoyama, T. 2011, ApJ, 740, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D. M. 2012, MNRAS, 427, 2487 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L. 2005, Physics Uspekhi, 48, 449 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Rüdiger, G. 1999, A&A, 344, 911 [NASA ADS] [Google Scholar]
 Küker, M., & Rüdiger, G. 2005a, A&A, 433, 1023 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Küker, M., & Rüdiger, G. 2005b, Astron. Nachr., 326, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Küker, M., & Rüdiger, G. 2007, Astron. Nachr., 328, 1050 [NASA ADS] [CrossRef] [Google Scholar]
 Küker, M., & Rüdiger, G. 2011, Astron. Nachr., 332, 933 [NASA ADS] [CrossRef] [Google Scholar]
 Lanza, A. F., Rodono, M., & Zappala, R. A. 1993, A&A, 269, 351 [NASA ADS] [Google Scholar]
 Marsden, S. C., Jardine, M. M., Ramírez Vélez, J. C., et al. 2011, MNRAS, 413, 1939 [NASA ADS] [CrossRef] [Google Scholar]
 Reiners, A. 2006, A&A, 446, 267 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reiners, A., & Schmitt, J. H. M. M. 2002, A&A, 384, 155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reiners, A., & Schmitt, J. H. M. M. 2003, A&A, 398, 647 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Roettenbacher, R. M., Monnier, J. D., Harmon, R. O., Barclay, T., & Still, M. 2013, ApJ, 767, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Walkowicz, L. M., Basri, G. S., & Valenti, J. A. 2013, ApJS, 205, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
Weighted means ⟨P1_{in}⟩ and ⟨P1_{out}⟩ and their associated errors σ(P1_{in}) and σ(P1_{out}) for both input and output periods for each noise level.
Weighted mean ⟨α_{out} − α_{in}⟩ and error σ(α_{out} − α_{in}) of α_{out} − α_{in} for all stars with two detected periods for different noise levels.
Comparison of periods returned by our method with spot periods from previously studied Kepler stars.
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: noisefree (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 P1_{in} and output periods P1_{out} for the noisefree case. The most significant period P1_{out} can be recovered fairly well. P1_{out} is on average 2.4% shorter than the actually contained period P1_{in}. Small plots: the distribution of P1_{out} 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 noisefree 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 noisefree 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.