Free Access
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/0004-6361/201321161
Published online 13 August 2013

© 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 . 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(Teff/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 ~ ⟨P1.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 (~) 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 P1out 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, (1)with I0 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 (2)with α = 0 supplying rigid body rotation. Ppole and Peq 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: (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.

Table 1

Stellar simulation parameters.

Table 2

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 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 Peq 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.

thumbnail 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.

Open with DEXTER

3.1. Lomb-Scargle periodogram

The generalized Lomb-Scargle periodogram (Zechmeister & Kürster 2009)1 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 χ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.

Table 3

Left: periods of the three spots from the light curve in Fig. 1 (left column) and output periods returned by the prewhitening analysis for the fit in Fig. 1 (right column).

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 P1out. 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 P1out 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 P1out and call it P1in.

Based on the period P1out, 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 P1out; 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) P1out = 1,P2 = 0.5; and 2) P1out = 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 P1out = 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 remaining2 sine periods. This period is called P2out. In 88 330 cases of our models we found a P2out that fulfills these criteria. Again, P2out was compared with all input periods, and the period closest to P2out was called P2in. When the closest input period picked was again P1in, P2out 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 P1out = 1.016,P2out = 1.184,P1in = 1.018, and P2in = 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.

thumbnail 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).

Open with DEXTER

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 P1out. 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).

thumbnail 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).

Open with DEXTER

4.1. Rotation periods P1out

96.2% of all light curves have one detected periodicity P1out that is the most significant one in the data. In Fig. 3 we compare the periods P1out with the input periods P1in for different noise levels to see how well our selection process works. The distribution of input periods P1in is shown in gray, the shaded distribution shows the output periods P1out. Because we only considered solar-like DR, the P1in 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 ⟨P1in⟩ and ⟨P1out⟩ in the above range in Table 4. In the noise-free case P1out is on average 2.4% shorter than the actually contained period P1in. For higher noise levels ⟨P1out⟩ 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. P1out 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 P1out = 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).

Table 4

Weighted means ⟨P1in⟩ and ⟨P1out⟩ and their associated errors σ(P1in) and σ(P1out) 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 P2out adjacent to the primary one and associated to a second spot period P2in is considered as evidence for DR. To investigate where this selection process is acceptable we considered the two pairs (P1out,P2out) and their associated periods (P1in,P2in). 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 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.

thumbnail 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).

Open with DEXTER

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 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%.

thumbnail 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.

Open with DEXTER

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 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 P2out, 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 Peq. We used our method to determine P1out and P2out to compare them with the reported results. We call the spot periods closest to our findings P1closest and P2closest. 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 P1out = 3.61 and P2out = 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.

Table 6

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 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 P1out was detected in 96.2% of all light curves. A second period close to P1out 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.


2

When P1out was a harmonic, the first two sine periods were excluded.

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

Table 1

Stellar simulation parameters.

Table 2

Spot simulation parameters.

Table 3

Left: periods of the three spots from the light curve in Fig. 1 (left column) and output periods returned by the prewhitening analysis for the fit in Fig. 1 (right column).

Table 4

Weighted means ⟨P1in⟩ and ⟨P1out⟩ and their associated errors σ(P1in) and σ(P1out) for both input and output periods for each noise level.

Table 5

Weighted mean ⟨αout − αin⟩ and error σ(αout − αin) of αout − αin for all stars with two detected periods for different noise levels.

Table 6

Comparison of periods returned by our method with spot periods from previously studied Kepler stars.

All Figures

thumbnail 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.

Open with DEXTER
In the text
thumbnail 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).

Open with DEXTER
In the text
thumbnail 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).

Open with DEXTER
In the text
thumbnail 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).

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.