Multiperiodicity, modulations, and flipflops in variable star light curves
III. Carrier fit analysis of LQ Hydrae photometry for 1982–2014^{⋆,}^{⋆⋆}
^{1}
ReSoLVE Centre of Excellence, Aalto UniversityDepartment of Computer
Science,
PO Box 15400,
00076
Aalto,
Finland
email:
nigul.olspert@aalto.fi
^{2}
Tartu Observatory, 61602
Tõravere,
Estonia
^{3} Department of Physics, PO Box 64, 00014 University of
Helsinki, Finland
^{4}
Finnish Centre for Astronomy with ESO (FINCA), University of
Turku, Väisäläntie
20, 21500
Piikkiö,
Finland
^{5}
Center of Excellence in Information Systems, Tennessee State
University, 3500 John A. Merritt
Blvd., Box 9501, Nashville, TN
37209,
USA
Received: 28 November 2014
Accepted: 25 January 2015
Aims. We study LQ Hya photometry for 1982–2014 with the carrier fit (CF) method and compare our results to earlier photometric analysis and recent Doppler imaging maps.
Methods. As the rotation period of the object is not known a priori, we utilize different types of statistical methods first (leastsquares fit of harmonics, phase dispersion statistics) to estimate various candidates for the carrier period for the CF method. Secondly, a global fit to the whole data set and local fits to shorter segments are computed with the period that is found to be optimal.
Results. The harmonic leastsquares analysis of all the available data reveals a short period, of close to 1.6 days, as a limiting value for a set of significant frequencies. We interpret this as the rotation period of the spots near the equatorial region. In addition, the distribution of the significant periods is found to be bimodal, hinting of a longerterm modulating period, which we set out to study with a twoharmonic CF model. A weak modulation signal is, indeed retrieved, with a period of roughly 6.9 yr. The phase dispersion analysis gives a clear symmetric minimum for coherence times lower than and around 100 days. We interpret this as the mean rotation pattern of the spots. Of these periods, the most significant and physically most plausible period statistically is the mean spot rotation period 1 60514, which is chosen to be used as the carrier period for the CF analysis. With the CF method, we seek any systematic trends in the spot distribution in the global time frame, and locally look for previously reported abrupt phase changes in rapidly rotating objects. During 2003–2009, the global CF reveals a coherent structure rotating with a period of 1 6037, while during most other times the spot distribution appears somewhat random in phase.
Conclusions. The evolution of the spot distribution of the object is found to be very chaotic, with no clear signs of an azimuthal dynamo wave that would persist over longer timescales, although the shortlived coherent structures occasionally observed do not rotate with the same speed as the mean spot distribution. The most likely explanation of the bimodal period distribution is attributed to the high and lowlatitude spot formation regions confirmed from Doppler imaging and Zeeman Doppler imaging.
Key words: stars: activity / starspots
Based on observations made with the Nordic Optical Telescope, operated on the island of La Palma jointly by Denmark, Finland, Iceland, Norway, and Sweden, in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias.
Table of minima is only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/577/A120
© ESO, 2015
1. Introduction
LQ Hya (HD 82558) is a chromospherically active BY Draconistype star of the spectral type K2V (Cutispoto 1991; Covino et al. 2001). It also shows a high level of Ca H&K emission (log , White et al. 2007), manifesting very high level of magnetic activity. With an estimated mass of 0.8 ± 0.1 M_{⊙} and age 51.0 ± 17.5 Myr (Tetzlaff et al. 2011), the star is considered a young solar analog. The star spins very fast, with the estimated rotation period being around 1.6 days (e.g., Jetsu 1993; Berdyugina et al. 2002; Kovári et al. 2004; Lehtinen et al. 2012).
In addition to exhibiting strong magnetic activity indicators, the star shows modulation in its light curve, as first proposed by Eggen (1984) and confirmed by Fekel et al. (1986). This behavior is interpreted as cool spots rotating with the stellar surface. For this reason, photometric light curves have been used to determine the rotation period of the star. This simple picture might not be applicable if the star exhibits latitudinal and/or radial surface differential rotation analogous to the Sun, or latitudinal dynamo waves, which is the solar case, as the sunspots form the wellknown butterfly diagram with cyclic behavior known as Spörer’s law, and/or azimuthal dynamo waves (which can occur in the rapid rotation regime, e.g., Lindborg et al. 2011). Indeed, based on previous studies of LQ Hya photometry, it has become evident that no single period suitable exists for describing the physical system throughout the whole observational time span available. However, for shorter epochs, dominating periods have been found. In Jetsu (1993) a good phase coherence was achieved with a period of , Berdyugina et al. (2002) arrived at a period estimate of , Kovári et al. (2004) report a period of , and the analysis of Lehtinen et al. (2012) gives a period .
The surface differential rotation of the object has been estimated by either using photometric light curves (Jetsu 1993; Berdyugina et al. 2002; You 2007; Lehtinen et al. 2012), or spectroscopic observations analyzed with Doppler imaging (DI) methods (Strassmeier et al. 1993; Saar et al. 1994; Rice & Strassmeier 1998; Kovári et al. 2004) and Zeeman Doppler imaging (ZDI) methods (Saar et al. 1994; Donati 1999; Donati et al. 2003b; McIvor et al. 2004). It is customary to define the differential rotation parameter as (1)describing both the magnitude and the type of the latitudinal rotation law. Large values of k denote strong differential rotation, positive signs corresponding to solarlike profiles with a faster equator and a slower pole, negative to antisolar profiles with faster poles and a slower equator. From photometry, only the magnitude of k can be deduced, while using DI one can also determine its sign. The values obtained from the fluctuations in the photometric period range from k = 0.015 (Jetsu 1993) to 0.025 by (You 2007).
The DI and ZDI results of Saar et al. (1994) initially estimated an upper limit of differential rotation based on polar “smearing” of k ≲ 0.03 and further analysis by Kovári et al. (2004), Donati et al. (2003a) indicate even weaker (k = 0.002...0.006) solarlike differential rotation. Because of the relatively low vsin(i) of the object versus the required value for DI techniques, the results show significant scatter (see, e.g., Barnes et al. 2005). Obtaining a reliable value for differential rotation using DI and ZDI requires a longer baseline of observations than used in most studies (Strassmeier et al. 1993; Rice & Strassmeier 1998).
In general, the amount of differential rotation in this object can be concluded to be very small compared to the solar value of k ≈ 0.2. Hydrodynamical meanfield modeling of the rotation law of this object by Kitchatinov & Olemskoy (2011) agrees with the observations in the sense that the obtained profiles are solarlike, and the magnitude falls within the observational range (k = 0.028).
The DI and ZDI maps (e.g., Strassmeier et al. 1993; Saar et al. 1994; Rice & Strassmeier 1998; Donati 1999; Donati et al. 2003b; Cole et al. 2014a) show both high and lowlatitude spot activity on the object to the extent that Cole et al. (2014a) define the distribution over latitude to be bimodal. The relative strength of the two latitudinal spot regions has been reported to be highly variable over time so that during some epochs the nearequator spots dominate, while during others the nearly polar features are the strongest. The spot distribution from DI and ZDI has been postulated to be concentrated onto active longitudes during some epochs. Therefore, the spot distribution has occasionally been highly nonaxisymmetric (e.g., Saar et al. 1994), but during different epochs no clear signs of active longitudes have been found (e.g., Cole et al. 2014a). In addition, Cao & Gu (2014) found rotational modulation of chromospheric emission during 2006−2012, which indicates the presence and change of chromospheric active regions over the surface of LQ Hya over a considerable timespan. Moreover, a close spatial connection between chromospheric active regions and photospheric regions from previous studies was detected.
Many authors of photometric or DI studies have reported the concentrations of the spot activity on certain longitudes (Jetsu 1993; Berdyugina et al. 2002; Kovári et al. 2004; Lehtinen et al. 2012). The periods producing the largest amount of phase clustering and the phase separation of the “active longitudes” vary significantly depending on the data set span and timing, and also the method used, indicating that these structures are not persistent, but rather may be changing strongly in time. For example, Lehtinen et al. (2012) concluded that during 1988−2011, two different periods describing structures with active longitudes and , the former appearing around 1995 and the latter between the years 2003 and 2009, were needed to describe the phase clustering of the data. In contrast, Berdyugina et al. (2002) used a light curve inversion technique to recover spot phases and postulated the existence of two active longitudes about 180° apart, spanning the entire 20 yr of the data set, and calculated a new rotation period for the spot structure from the drift of these active longitudes, .
The mean brightness of the star exhibits obvious signs of cyclic activity. Various determinations from photometry indicate cycles of around six to seven years (e.g., Jetsu 1993; Strassmeier et al. 1997; Cutispoto 1998; Oláh et al. 2009), while also coexisting shorter (approximately three years, e.g. Messina & Guinan 2003) and longer (approximately 11 yr, e.g., Oláh et al. 2000) cycles have been reported. Berdyugina et al. (2002), who postulated coherent active longitudes during the time span of almost 20 yr, found a 7.7yr cycle in the mean brightness, and additionally a 5.2yr cycle related to the regular change of the activity level of the two active longitudes (flipflop).
From various earlier studies, it is evident that the behavior of the object is extremely complex and a systematic approach to understanding these complexities is yet missing in the literature. This is attempted in our study using the carrier fit (CF) method, the essential properties of which are explained in Sect. 3.1. This method has been previously successfully applied to study spot activity in other types of stars: see Hackman et al. (2013) for the analysis of FK Coma Berenices and Lindborg et al. (2013) for the analysis of II Pegasi. As the rotation period of the star is not known a priori, we search for the optimal carrier frequency first using different kinds of statistical methods described in Sects. 4.1−4.3. We then perform a global CF using the optimal carrier frequency (Sect. 5.1), the aim being to search for any persistent trends and/or phase disruptions. Finally, we perform local CFs to study local segments in an attempt to identify interesting phase behavior. Here we also make a comparison to the continuous period search (CPS) method results by Lehtinen et al. (2012) and the recent DI results of Cole et al. (2014a). In Sects. 6 and 7, we discuss and conclude our findings.
2. Data
We use data consisting of nearly 32 yr of photometry from three different sources, namely the photometry collected and published in Berdyugina et al. (2002) from HJD = 2 445 275 (2 November 1982) to HJD = 2 452 053 (23 May 2001); the published photometry of Lehtinen et al. (2012) from HJD = 2 447 141 (11 December 1987) to HJD = 2 455 684 (2 May 2011) obtained with the T3 0.4 m APT at the Fairborn Observatory, Arizona; and finally unpublished photometry obtained with the same telescope from HJD = 2 455 685 (3 May 2011) to HJD = 2 456 783 (5 May 2014). From these data sources (see Fig. 1) we compiled two input data sets. First, we rescaled T3APT data to fit the data from Berdyugina et al. (2002) and for overlapping observations, computed data points as averages (maximum allowed time difference of 0.1 days was used). As a result of this procedure we got a long data set (hence D1) consisting of 3929 observations covering 11 508 days. We used this data set to perform periodicity analysis and optimal carrier period value estimation; the analysis is presented in Sect. 4. For the CF analysis presented in Sect. 5, we combined all data from T3APT telescope, i.e., the data published in Lehtinen et al. (2012) supplemented with the seasons up to May 2014. As a result we got a homogeneous data set (hence D2), which consisted of 2907 observations covering 9642 days. We denote this data set as homogeneous, as it was observed with the same instrument with the same comparison star (HD 82428). We also note that it is of better quality than the D1 data.
Fig. 1 Combined data set from the three different sources. The data set referred to as D1 consists of all available data, while data set D2 comprises only T3APT data (Lehtinen + unpublished). 

Open with DEXTER 
3. Carrier fit method
3.1. Overview
A detailed description and the discussion of applicability of the CF method can be found in Pelt et al. (2011); here we briefly cover some of key aspects of this method. A CF model can be described as a truncated, slowly modulated harmonic decomposition of the signal: (2)where ν_{0} is the carrier frequency, a_{0}(t) is the timedependent mean level of the signal, K is the total number of harmonics included in the model, describing the overtones of the basic carrier frequency. The lowfrequency signal components, a_{k}(t) and b_{k}(t), can be modeled by either splines or harmonics. Depending on the time series in question, either one of these approaches may be a more suitable choice, e.g., the spline approximation might be more suitable for cases where the signal is known to abruptly change. However, generally the difference between goodness of fits (for definition, see Sect. 3.2) is marginal if the number of free parameters for both methods are approximately the same. In the scope of the current study, we limit ourselves to the option of harmonic modulators.
The value for the carrier frequency used in the above model can be, for instance, the rotation period of the object, or any other clocking frequency describing the system. In the case of LQ Hya, however, the rotation period is not directly known, as it is a single star. Therefore, we discuss various approaches to carrier frequency (or period) selection below.
Following the notation used in Pelt et al. (2011), the trigonometric approximation of the slowamplitude modulation curves can be written as: where L is the total number of harmonics used in the modulator model and ν_{D} = 1 /D = 1 /C·(t_{max} − t_{min}), where D is data period, C is the coverage factor, and [ t_{min},t_{max} ] is the time interval to be fitted with the model. The data period D must be significantly longer than the carrier period P_{0} = 1 /ν_{0}, and preferably it must also be a little longer than the actual time span of the data, i.e., the coverage factor should be .
3.2. Selection of free parameters
One of the crucial things to consider when applying the CF method to real data are the suitable values of the free parameters (ν_{0}, K, L and C). As the selection of the carrier frequency ν_{0} is particularly important, we dedicate Sect. 4 to it. The methods used for determining the values for other parameters will be discussed here. In general we need to run the computations with all possible combinations of parameter values drawn from some meaningful ranges and then estimate the goodness of fit for every run using (5)where n is the number of data points, y_{i} is the value of the ith data point, f_{i} is the value of the fit corresponding to the time moment of ith data point, and is the mean of the values of all data points. Qualitatively speaking, the goodness of fit compares the variance of data points around the model to that of the data itself. In the current study, we use two approaches: we start by analyzing the full set of data as a whole (global fit) followed by the analysis of seasonal data segments (local fit). In both cases, we aim for as high R^{2} values as possible while avoiding the possibility of either overfitting the data or fitting into the gaps.
Before continuing with the methods of parameter selection we note that the suitable value for L is, first of all, dependent on the value of the coverage factor C, which defines the period of the slowest modulator in the model. It is reasonable to adopt a value of the same order or little longer than the length of the whole data set. This way we guarantee that the slowest detectable changes in the data are taken into account by the model. In the current study we have fixed C = 1.2.
The difficulty introduced by the gaps in the data constitutes the socalled cycle count problem. The lowfrequency modulators L introduce variance around the carrier frequency. Here, we need to keep in mind that the difference in cycle counts for these maximum and minimum frequencies during the longest gap in the data should be less than one to avoid phase match indeterminacy. This can be concisely expressed by the following criterion: (ν_{0} + ν_{D})Δ_{gap} − (ν_{0} − ν_{D})Δ_{gap}< 1, where ν_{0} is a highfrequency carrier, ν_{D} is a low frequency modulator, and Δ_{gap} is the length of the longest gap in the given data set. After replacing and simplifying ν_{D} = L/D = L/C·(t_{max} − t_{min}), we obtain an estimate for the upper limit for L: (6)Using this formula, we can estimate the maximum valid L for the data set with the given time span and the longest gap.
In case of local fits, because of the low number of data points, there exists a possibility of an overfit. We use Bayesian information criterion (BIC) to determine the optimal values for K and L in a similar way that was done in Lehtinen et al. (2011), the differences being that we have two parameters in the model and we omit the weights of the data points. More precisely, we search for the minimum of the following criterion: (7)where , and we have used the same notation as in Eq. (5). The first term of this equation describes the quality of the fit and the second one adds a penalty proportional to the total number of parameters in the model. Now, as we need the information for possible primary and secondary minima to be able to detect flipflop type events, we omit K = 1 from the set of trial values. For L we do not impose a lower limit, but we will calculate the upper limit using Eq. (6) to check if the value obtained from BIC is valid. Consequently we will use the minimum of the values obtained from Eqs. (6) and (7) as the optimal value. The practical application of this procedure is detailed in Sect. 5.2.
In case of the global fit, the possibility of overfitting is quite low as the presence of long seasonal gaps in the data significantly lowers the maximum possible values for K and L. It is quite probable that when increasing the number of parameters, the model starts showing big distortions in the regions where data is missing considerably earlier than the high value of R^{2} (e.g. 90 %) is achieved.
Before starting to measure the effect of the gaps on the model, first we need to specify how long a region without data qualifies as a gap. In our case, the data is divided into observational seasons where relatively densely spaced data is alternating with slightly shorter ranges with no data at all. Based on a closer look at the actual spacing of the data, we define the gap as any region without data that is longer than 130 days. This definition leads to 27 segments with data and 26 gaps between them (minimum being 139 and maximum 302 days). The length of the homogeneous data set is 9642 days, thus using Eq. (6) with the maximal gap size of 302 days we obtain a maximum value for L to be 19. As explained in Sect. 4.1, this is too low value to cover the full spectrum of the data. To eliminate this problem, one option would be to leave out all data points preceding the longest gap and work with the truncated data set. By not doing so, on the negative side, we lose the reliability of the model during this single gap, but on the positive side, we still maintain a global description for the maximum available time span. We must note, however, that even when neglecting the cycle count problem, the reliability of the model is still expected to decrease during all seasonal gaps. Based on this consideration, we decided to keep the full data set and expect the error estimates to reflect the reliability of the model sufficiently well (please refer to Sect. 3.4 for more details about error estimation).
The second longest gap in the data is 192 days long, increasing the suitable value for L to 30. This number is also close to the number of observing seasons, so that we could expect good approximation of the seasonal variation by the term a_{0}(t) in Eq. (2). In choosing the suitable values of highfrequency components K, there is not much room: on one hand K = 2 is the lowest value meaningful in our analysis because of the same considerations as pointed out in case of the local fits. On the other hand, tests with K = 3 showed only a small positive effect in the achieved R^{2} value, while significantly increasing the freedom of the model (distortions) in the region of gaps. Based on these arguments, we decided to fix K = 2. The total number of parameters in our model is therefore 305, which means approximately 10 data points per parameter.
After we fixed the optimal parameter values for the model, we can further increase the goodness of fit by removing the 3σ outliers to the initial fit from the data and then refitting again. The outliers are either observationally unreliable points or possible flares, such as one that occurred around April 2000 or HJD 2 451 650. In our case, we detected a total of 22 outliers. The removal of these outliers leads to an approximately 3% increase in R^{2} for global fit. Therefore, in all subsequent CF analysis (for global as well as for local) we used the data set with outliers eliminated.
3.3. Visualization
To visualize the CF model, we use the same technique as introduced in Pelt et al. (2011, p. 4, Sect. 2.6). Firstly, we divide the whole time span into the number of bins with the length of the chosen carrier period. Secondly, for each bin we normalize the signal amplitude into range [−1,1] and then plot it with the corresponding time moment of the bin and the phase relative to the carrier period. This normalization step is useful for making the phase behavior of the signal comparable over the whole time span. Without normalization the features during high amplitudes will dominate the picture. Here we use both approaches for the purpose of obtaining more information about the processes governing the star. At the bottom of the plot we include a socalled “barcode” to give information of the density of the data points around the given time moments. Black indicates densely spaced data, while yellow indicates sparsely spaced or no data at all. Some previous examples using the given technique can be found in Hackman et al. (2013), Lindborg et al. (2013).
3.4. Minima detection, error and significance estimation
Besides using CF method for visualizing the results, we determine primary and secondary minima from the model of global fit and compare them to results obtained from DI analysis and the earlier analysis of a shorter segment of the same data used here with a different method (CPS; Lehtinen et al. 2012). Error estimates for the minima are calculated by generating 1000 bootstrap samples from the original data, by reshuffling the residuals of the data points to the initial CF model allowing recurrences; repeating the CF analysis for each new data set; and finally, obtaining the distributions for the minima. We mark a minimum as being reliable if and only if the following two conditions are satisfied for distributions both in time and magnitude: the KolmogorovSmirnov test with preassigned significance level 0.01 against a normal distribution must pass and the bias of the mean of the distribution from the original estimate should be less than the standard deviation of the distribution. The error and significance estimates for each minimum can be found along with the full global fit data in the material published online.
4. Methods and results for searching the optimal carrier frequency
The selection of the optimal carrier frequency to analyze LQ Hya light curves is far from a trivial task. There is a wide range of different periods obtained from different data sets and using different methods. This is why we need a thorough analysis of the periods to proceed with the CF method. Building on the solar analogy, one might hypothesize that there are at least four types of periodicities that one needs to deal with:

1.
Behind all the observable activity is the stellar surface rotationand its nonuniformities. In the solar case, the sunspots quiteclosely follow the motion of the solar surface seen fromDopplerograms; the only exception are the longitudinal activitynests that have been reported to show motions that differ from thegeneral pattern. In the absence of asteroseismic data, it isimpossible to distinguish between the motion of the plasma andthe spots themselves, and the determination of the rotation periodis out of the scope of this study.

2.
Through photometry, one can hope at least to be able to determine the mean rotation period of the spots on the stellar surface, analogous to the Carrington rotation period of the Sun.

3.
Taking the DI and ZDI results of significant spot activity in the lowlatitude regions and solarlike rotation law of LQ Hya as valid assumptions, one may also be able to determine the spot rotation period near the equatorial region, as this should arguably be the shortest period of rotational origin seen in the periodograms.

4.
Analogous to the periods derived for the longitudinal activity nest on the solar surface for the azimuthal dynamo waves on some more evolved rapid rotators, photometric studies can be used to determine whether any longitudinal clustering occurs and if the period of the active longitudes differs from the mean rotation period of the spot distribution.
4.1. Leastsquares fitting
For irregularly spaced data sets the most common procedure of frequency analysis is a simple leastsquares fit of a harmonic waveform into the data with the range of some trial periods. The particular computational schemes and descriptive statistics vary (see, e.g., Barning 1963; Lomb 1976; Scargle 1982). In our analysis, we chose the simplest statistic, which measures the relevance of the harmonic under discussion, namely its amplitude.
Fig. 2 Amplitude spectrum for the important period range 1.5−1.7 days (gray) and first 100 frequencies obtained by sequential prewhitening procedure (black). 

Open with DEXTER 
Fig. 3 Set of 28 strongest amplitudes (black) and original amplitude spectrum (gray). The horizontal line marks the cutoff level obtained from 100 randomized samples. 

Open with DEXTER 
Before computing the amplitude spectrum, we removed the seasonal means (detrending) of the full LQ Hya Vband photometry data. The complexity of the spectrum, depicted in Fig. 2 in gray, is obvious and it is very hard to single out any prominent peak from the forest of peaks in the spectrum. The analysis can and must be improved by removing the spurious periods rising from the gapped nature of the data. For this, we used the socalled prewhitening method (for a recent application of the method, see, e.g., Reinhold & Reiners 2013). We iteratively removed the most significant harmonics with estimated amplitudes, and proceeded in the next step with the least squares fit residuals. In this way we computed a set of the 100 strongest amplitudes, which are depicted in Fig. 2 in black. These are not, with high probability, aliases due to the most prominent yearly gap structure with frequency offsets of Δν ≈ 1 / 365. Even after this cleaning procedure, we still have a very complex set of different frequencies.
A large part of the estimated frequencies with corresponding amplitudes are well below the noise level and we can disregard them. To estimate a suitable cutoff level, below which the peaks are considered insignificant, we proceeded in the following way. We built artificial data sets from the original by reshuffling them in time. For every new data set, we found its strongest peak. We then chose the amplitude cutoff level as a value of maximum amplitude in 100 different random runs. This criterion is rather conservative and we can be quite confident that the 28 periods, shown in black lines in Fig. 3, whose amplitudes were higher than the cutoff level (0.00385; indicated with a horizontal line in Fig. 3), are not the results of random fluctuations. As seen from the plot, the set of the selected peaks is now much more localized. Among the peaks are practically all the periods that have been proposed so far by different authors for different data sets, see Table 1.
Periods with strongest peaks in the spectrum compared to previous estimates.
The selected set of periods lay in the interval 1.598406...1.622572 days (or in frequency terms 0.6163053...0.6256234 cycles per day). The center of the full range interval is at and this is the first logical candidate for the carrier frequency. The logic behind this choice is obvious: the full frequency range will be covered on equal grounds.
Several authors have already reported on the period variability of LQ Hya. By using local fits, e.g., Messina & Guinan (2003) give the period in the range and You (2007) report range from their analysis. The obtained range of periods is also in good agreement with a set of rotation periods obtained from a simple spot modeling procedure by Alekseev (2005). The time span for our main homogeneous data set (D2) is around 9642 days. Correspondingly, the range of rotation counts for the significant periods is 5942...6032 with a difference of 90 rotations. For the carrier frequency around the center of the estimated range, this allows up to 45 full phase cycle runs in both directions. In Sect. 3.2, however, we concluded that the optimal harmonic count for lowfrequency modulation curve is around 30 cycles. It follows that, in principle, the very sharp (momentous) frequency jumps from one side of the range to the other can become smoothed out to some extent.
On the other hand, the second longest seasonal gap in the data is around 192 days and corresponding cycle counts are 199...120. Obviously, for the carrier periods in the middle of the full range, the phase migration during the gaps is certainly less than a full turn, cf. Eq. (6). For the longest gap in the data (around 302 days), however, we can theoretically have a cycle count error and consequently the approximated solution in this region can not be regarded reliable.
Finally, one interesting observation that can be seen in Fig. 3 is that the distribution of the significant frequencies is bimodal, i.e., there are two bunches of them. As seen from Table 1, typical solutions are concentrated in the rightmost bunch of periods, the cutoff being very sharp on the higher frequency (shorter period) side.
The bimodal structure of the period distribution leads us to another hypothesis that this is a case when a certain frequency in between the two bunches is more or less periodically modulated (with period around 2000−3000 days). This hypothesis was already set up in an earlier paper (see Berdyugina et al. 2002). To check it again, we carried out the corresponding analysis for our significantly longer data set D2.
Fig. 4 Folded and stacked light curve model with periods and . This regular structure helps to describe only R^{2} = 23.542% of overall variability. 

Open with DEXTER 
Fig. 5 Carrier fit obtained with the carrier obtained from the twoperiodic model of Sect. 4.2. 

Open with DEXTER 
4.2. Carrier from a multiperiodic model
The simplest conceivable model for the slowly modulated signal is a time series that depends on the carrier period and modulating period P_{mod} = 1 /ν_{mod} in a coupled way. That means that all positive combination frequencies (8)can take part in waveforming. In the simplest case of N = 2, there will essentially be 25 different trigonometric terms (constant included) that must be fitted into the observed data using the leastsquares method (see Berdyugina et al. 2002). The periods producing the best fit will then be used to build a final model.
We performed this kind of analysis for the data set D2. The resulting pair of periods occurred at and (roughly 6.94 yr). The resulting fit can be visualized by the same devices as the CF. We divide the solution into period length fragments and stack them properly. In Fig. 4 our solution is depicted. The R^{2} value for the least squares fit for these two periods is quite low (23.6%). It is remarkable that the obtained modulation period is very similar to values often quoted as the long cycle length of LQ Hya (e.g., Jetsu 1993; Oláh et al. 2000; Berdyugina et al. 2002; Alekseev 2005). However, it is known that given cycle length is nonstationary (e.g., Oláh et al. 2011), partially explaining the low quality of the fit in our analysis.
If we now base our consideration on the idea that the possible doubly periodic component is a relevant aspect of the overall variability then we can use the computed carrier value as a modelbased carrier for the CF procedure. The result of this analysis is shown in Fig. 5. Unfortunately, the regular structure of the simple model is largely lost and a number of other variability elements dominate the picture.
4.3. Carrier from the phase dispersion analysis
From Fig. 3 we can well see that the period computed as a mean from the formal range of significant periods is not very representative because of the bimodal nature of the distribution. This is also true if the peak with maximum amplitude was selected as a carrier. A slightly better method would be the one used in Lehtinen et al. (2012) where the best period was selected using phase distributions of the light curve extrema. However, in this case, the interplay between different frequencies can shift minima or maxima and obscure the general picture.
To our understanding, the best method for computing the mean period of spots comes from phase dispersion analysis (Pelt 1983; Lindborg et al. 2013). It is based on the following simple statistic: (9)where f(t_{i}),i = 1,...,N is the input time series, σ^{2} is its variance, g(t_{i},t_{j},P,Δt) is the selection function, which is significantly greater than zero only when In the latter condition, Δt is the socalled correlation length. For the particular case when Δt is longer than the full data span, the D^{2}(P) statistic is essentially a slight reformulation of the wellknown Stellingwerf statistic (Stellingwerf 1978). As the correlation length is made shorter, we match nearby cycles in a progressively narrower region, and consequently estimate a certain mean period, which need not be coherent for the full time span. This is well illustrated in Fig. 6, where we show the D^{2}(P) statistic for the range of trial frequencies as function of the correlation length, Δt, using color contours.
Fig. 6 Phase dispersion statistic D^{2}(P) for a range of correlation lengths. The dispersion spectrum becomes strongly asymmetric around 230 days and then splits into set of peaks. 

Open with DEXTER 
We see that for a small enough correlation length the D^{2}(P) statistic yields a rather symmetric single minimum. Above that limit the frequency spectrum starts to distort and eventually splits into separate branches. Finally, at large correlation lengths we obtain a forest of minima similar to the results presented in Sect. 4.1. We interpret this behavior in the following way: for short coherence times, the periodogram is dominated by the mean pattern of spot motions, while at longer coherence times the signatures of more persistent spot structures, the rotation of which differs from the mean spot flow, take over and give the strongest signal.
Our aim is to determine the limiting correlation length at which a single minimum is still obtained and use the corresponding period of the phase dispersion minimum as a plausible carrier period. Here we must point out two problems: firstly the minimum of D^{2}(P) statistic is usually very wide, and secondly the shape of the peak somewhat deviates from symmetric form (see for details Fig. 7).
Fig. 7 Curves of the phase dispersion statistic D^{2}(P) for different correlation lengths in days: 100 (black), 200 (gray), and 300 (light gray). Dashed curve corresponds to the bestfitting Gaussian to the curve with correlation length 100 days. 

Open with DEXTER 
This prevents us from determining a sufficiently accurate value of the minimum directly from the curve of the statistic. Instead of that, we fit a Gaussian profile to the curve of D^{2}(P) statistic. Our task is to estimate the free parameters of this Gaussian curve (mean μ and variance σ^{2}) for which the distance to the curve of the D^{2}(P) statistic is minimal. The value of the mean obtained this way represents the optimal carrier frequency and the variance represents the scatter of the periods around it. From our analysis, it turns out that, for a correlation length of 100 days, the curve of the D^{2}(P) statistic is singular and relatively symmetric. Hence, we choose this correlation length as the limiting length, and determine the mean period from this curve. The bestfitting Gaussian has the mean μ = 0.62300 and standard deviation σ = 0.00325 (both in cycles per day). In the time domain the corresponding values are and . For comparison, a similar value, , was obtained by Lehtinen et al. (2012) by calculating the weighted mean of the periods determined for independent subsets of the whole data.
Fig. 8 Phase diagram for global fit with carrier period 1 60514, K = 2 and L = 30. On left with normalized amplitudes, on right with actual magnitudes 

Open with DEXTER 
We estimated the significance of the mean cycle length by testing the null hypothesis that the peak of the minimum is drawn from the distribution of random fluctuations. For that purpose, we generated 1000 samples from the original data set via reshuffling the measured magnitudes and then calculated the value of the D^{2}(P) statistic for each new data set using correlation length of 100 days. In this manner, we obtained a distribution for the minimum of the D^{2}(P) statistic caused by random fluctuations. The results showed that in our case the null hypothesis can be rejected with preassigned significance level of 0.05 as the minimum of the D^{2}(P) for the original data remains well below the 995th value from the sample distribution, which was located around 0.97.
We believe that the carrier value obtained using the D^{2}(P) statistic, which is computed for short correlation lengths, is well grounded. The other values (the strongest peak in the spectrum, central value of the range, mean value of the set of periods, periods based on extrema) tend to bring in an amount of contingency, as using any of them would in practice mean choosing a certain locally active period for the CF analysis of a global nature.
5. Carrier fit results
In the following we will use the obtained mean cycle length period of the spots as the carrier period in Eq. (2), substituting , and compute models of global as well as local fits.
5.1. Global fit
We perform a global CF analysis with the selected parameters , K = 2, L = 30, and C = 1.2, which gives a model with R^{2} = 91.7%. Consequently, the CF model describes a rather large part of the overall input data variability. The resulting phase diagram is displayed in Fig. 8, from which it is evident that the phase behavior is characterized by up and downward trends, while epochs during which the minima would occur at constant phases are almost totally absent. On the one hand, the slopes of the trends correspond to the different periods, and on the other hand, the lengths and the locations of the trends give hints of the “duty times” of these cycles (in other words when and for how long these periods dominate). The most noticeable feature is a downward trend during the years 2003−2009, which roughly corresponds to the period 16037 days, belonging to the strongest peak found previously by frequency analysis. Another prominent downward trend occurs during the years 1990−1994. Average period describing this feature is around 16033 days. During other epochs the phase behavior is changing more rapidly in time while the abrupt phase shifts seem to appear over the whole time span, even during the strongest trends.
Summary of the local CF analysis results.
Fig. 9 Phase diagrams of local fits with refined carrier period 1 60514. 

Open with DEXTER 
5.2. Local fits
In the current study, local analysis of the data is carried out mainly for the purpose of visualizing the phase behavior of the signal in more detail than it is possible with the global analysis. Because there are insufficient number of modulators in the global fit, the whole spectrum of the data is not covered. This leads to the smoothing of the real signal on small timescales. To get better results, we perform a local CF analysis similar to that of Hackman et al. (2013) and Lindborg et al. (2013). Segments contain a relatively low number of data points and in some cases considerable gaps are present, thus the maximum allowed number of modulators is significantly smaller than that for the global fit. Therefore one could assume that no higher precision in the results is achievable. However, the reality is exactly the opposite, we obtain better coverage of the spectrum for two reasons: firstly, the spectrum of each separate segment is narrower than that of the whole data set, and secondly, the frequencies of the modulators in the model are much higher than those used in case of the global fit.
Splitting the data into segments suitable for the local analysis is done using the following rules: we start with the earliest data point and calculate the difference in time between pairs of data points next to each other d_{i} = t_{i} − t_{i − 1}. If d_{i}< 130d the data point at t_{i} is added to the segment, otherwise a new segment is started and the process is repeated. Using the above algorithm we obtained 27 segments, the details of which are given in Table 2. The CFs for each segment were applied with the period of the mean cycle length . We determined the number of harmonics K and modulators L optimal for each segment by finding the minimum for the BIC. For all segments, K = 2 turned out to be suitable except for the segment 5, which had only 27 points so that realistic modeling was impossible. Because of the cycle count problem, L was further decreased for the segments 1 and 26. The exact number of parameters L used in the model and the R^{2} value achieved for each segment are summarized in Table 2. We can see that for most of the segments the goodness of fit is even higher than 90%, but for segments 3, 10, 11, and 26 it is quite low. In case of the segment 26, this can be explained by a big gap in the data and the low number of points. The other above mentioned three segments are, however, quite densely populated. This might be an indication that the signal is more complex in these segments than what is possible to model with the given number of data points. The resulting phase diagrams for all 27 segments are shown in Fig. 9, where the primary minima can be found by following black or dark blue features, while secondary minima appear either as red features between yellow features, or violet features between red features.
In Table 2 occurrences of flipflops are marked with the “+” symbol, totaling four events in segments 2, 11, 13, and 24. Other types of disrupted phase behavior events are marked with “?”, constituting the segments 1, 6, 8, 9, 12, 15, 22, 24, and 27. These are either single phase jumps of primary minima or swaps between primary and secondary minima less than 0.5 in phase. Clear upward trends can be seen in segments 6, 12 14, 15, 22, and a relatively gentle downward trend in segment 16. These are marked with the symbol “/”’ in the same table.
A similar analysis was carried out by Lehtinen et al. (2012) using the CPS method for the same data set, except for the last three seasons. The results of this study are in agreement with ours. Most of the interesting features can be seen on the phase plots from both studies. Some differences occur for the segments 3 and 6 where some of the minima detected by CF are absent in the case of CPS. For segment 19 there is no secondary minima from CPS, but it can be seen during the first 20 days in the case of CF. No comparison between the results is available for segments 4 and 5 because of the low number of data points.
We have also calculated epochs of possible flipflop events from the global fit using the definition from Hackman et al. (2013):

the region of main activity shifts about 180 deg from the old active longitude and then stays on the new active longitude; or

the primary and secondary minima are first separated by about 180 deg, after which the secondary minimum evolves into a longlived primary minimum, and vice versa.
Two additional restrictions were added to the above scheme: firstly, we counted only those events for which the phase shift lies between 0.45 and 0.5; secondly, the primary and secondary minima at the moment of flipflop must be reliable according to the error estimates from bootstrap runs. The results show that four of the total six flipflop events detected from global CF reside within the data, while two of these events occur in gaps. Moreover, three of these events are located in the segments 2, 13, and 24 for which we have also detected flipflops from local fits. One of the flipflops, namely within segment 11, was detectable only from the local fit. The epochs of all seven detected events are depicted in Fig. 10 with thick green vertical lines. Example of the global CF model around the flipflop event in segment 2 can be seen in Fig. 11. Following the light curve from left to right, we notice that the magnitude of the primary minima decreases while that one of the secondary minima increases. After around HJD 2 447 490 both minima “swap” their magnitudes, the change in phase corresponding to 0.5.
In Berdyugina et al. (2002) a 5.2 yr flipflop cycle was reported. In the light of our current study, this periodicity cannot be confirmed: on the one hand, this is because of the small number of events detected, and on the other hand, some of the detected events are only separated by two or three years.
Fig. 10 Phases of the minima for period . Red dots: primary minima from CF; red pluses: primary minima from CPS; blue points: secondary minima from CF; blue crosses: secondary minima from CPS; orange circles: minima from DI; and bold green vertical lines denote the epochs of possible flipflop events. The black bars with the numbers correspond to the photometric observing seasons, the numbers on top of the figure mark the observing seasons of DI. 

Open with DEXTER 
5.3. Comparison with Doppler imaging
Using our CF model for a global fit with the period of mean cycle length of 1 60514, we determined the epochs of primary and secondary minima. In Fig. 10 these epochs are plotted against the phase of the same period. Larger red and smaller blue dots represent the locations of primary and secondary minima, respectively. As noted above, 1000 bootstrap samples were calculated for the global CF to obtain error estimates for the minima. However, for a clearer visualization these are omitted from the figure, but are included in the online material, where the global CF solution is given.
In Cole et al. (2014a) the DI technique was applied to spectrometry of seven observing seasons. The resulting surface temperature maps were used to determine the epochs of the temperature minima, interpreted as starspots, for every season. Especially high activity level, i.e., a large number of cool spots, was observed to occur during October 1999 to November 2000. In Fig. 10 we plotted the retrieved spot epochs with orange circles, where the size of the circle reflects the temperature of the given spot (a larger circle corresponds to a lower temperature). On the same plot we also included the minima obtained from the CPS method published in Lehtinen et al. (2012). Red pluses and blue crosses represent the primary and secondary minima, respectively.
As expected the minima obtained from the global fit serve as the averaged values for the minima obtained from the CPS. Agreement with the results from DI is also quite satisfactory. There is a quite good match between DI and other models for some of the minima found in seasons 3 to 7. Even though neither active longitudes nor flipflop type events were seen in DI, both global and local CF analysis reveal a possible flipflop. However, it is not reasonable to search for the full agreement between photometry and DI. For instance, the photometry is affected by the limb darkening and surface area projection effects of the active regions, so that onetoone correspondence between the strongest minima in photometry and the lowest temperature regions in DI cannot be expected. We also note that the observing seasons for spectrometry and photometry mostly do not overlap. In addition a low signaltonoise ratio as well as less than ideal phase coverage for several observing seasons was reported, which increases the uncertainties even further (Cole et al. 2014a).
Fig. 11 Zoomin to the global CF near the flipflop event detected in the segment 2. The data points are drawn as small rectangles. 

Open with DEXTER 
6. Discussion
If we make an assumption that the scatter around the mean period is the result of differential rotation, we can use the halfwidth of the minimum as an approximation for the standard deviation. Based on this, we estimate the differential rotation coefficient using the formula by Jetsu (1993)(12)where in our case is taken as the σ of the closest Gaussian curve to dispersion statistic with correlation length of 100 days. Substituting the value into the equation leads to k ≈ 0.032, which is in rough agreement with previous estimates from photometry (e.g., Lehtinen et al. 2012; You 2007; Jetsu 1993), k = 0.015...0.025, but significantly larger than the values obtained from DI (e.g., Kovári et al. 2004; Donati et al. 2003b), k = 0.002...0.006. If we assume that the latitudinal dependence of differential rotation is solarlike, i.e., follows the law (13)an estimate of the latitude corresponding to the mean “Carrington” period can be calculated using the indicative “equatorial” period obtained from the leastsquares method () and the differential rotation coefficient derived above. This procedure yields θ_{C} ≈ 21°. The Carrington latitude would be considerably larger (30° up to the pole) if a differential rotation coefficient from earlier studies, which all are smaller than the value derived here, were used. Both the low and high values of θ_{C} match with the latitudinal location of the spots derived by DI during partially overlapping time epochs (Cole et al. 2014a). Although the bimodal distribution of periods in the leastsquares analysis shows more power in the short period (high frequency) bunch of peaks, indicating for a slight dominance of the lowlatitude spot region, the situation remains inconclusive and speculative. From theoretical grounds, spot formation occurring in comparable latitudes to the Sun (solar Carrington latitude 26°) is rather unlikely, as in rapid rotators the strong influence of Coriolis force is generally thought to lead into highlatitude or polar spot formation(see, e.g., Schuessler & Solanki 1992). However, coexisting polar and lowlatitude spot activity has been found in the models of Işık et al. (2011), which combine simple axisymmetric deepseated dynamo models with 3D flux transport models. This is especially apparent in their models of rapidly rotating solarlike and Ktype mainsequence stars with a comparable rotation rate to LQ Hya. These models also predict a dominance by the polar region.
The Coriolis number is an estimate of the strength of the rotational influence over turbulent convection, and can be written as (Saar & Brandenburg 1999) (14)where P_{rot} is the rotation period of the star and τ_{c} is the convective turnover time. Ossendrijver (1997) presented an extrapolation method to the theoretical calculations of Kim & Demarque (1996) to estimate the convective turnover time from the B − V index. Lehtinen et al. (2012) applied this technique and arrived at an estimate of τ_{c} ≈ 33.5 days for LQ Hya. On the other hand, they used their CPS method to compute the time scale of change, denoted with T_{c}, for each individual data set investigated, and as an average of all the analyzed segments, found a value of 50.5 days. As this quantity describes the typical time in which changes in spot configuration occur, it can be postulated to have some relation to the convective turnover time.
The Coriolis number of LQ Hya, based on the above stated values of the turnover time, lies within the range 260−400. These numbers are huge in comparison to the Sun with Co of roughly 6 with the definition Eq. (14). In the study of Saar & Brandenburg (1999) stars were observed to cluster on certain activity branches, when their Coriolis number and rotational vs. magnetic cycle periods, P_{rot}/P_{cyc}, were plotted. LQ Hya was termed an anomalous object, falling in the transition region between the active (A) and superactive (S) branches, based on the values τ_{c} ≈ 20.9 days by Gunn et al. (1998) and P_{cyc} ≈ 7 yr by Strassmeier et al. (1997) used back then. The majority of magnetic cycle determinations still falling into the range of six to seven years, together with the usage of the alternative methods to determine the convective turnover time to compute the Coriolis number, would place the object to an even more anomalous place in the diagnostic diagram, clearly further away from the other stars (although only a few of them were identified) in the transition region. This makes LQ Hya a very fascinating object to follow up and study further. These considerations, on the other hand, might also indicate that the division into active and superactive branches in the diagnostic diagram is not as meaningful as the separation into the inactive and active branches (the socalled VaughanPreston gap).
The difference between the retrieved mean rotation period of the star and the most pronounced coherent phase structure during 2003 and 2009 () is roughly 0.00171 days. During the aforementioned years, the trend in the phasetime diagram is nearly linear, the spot structure going faster (nearly linear downward trend) with respect to the mean rotation of spots. Previously, this behavior has been seen on more evolved rapid rotators (Lindborg et al. 2011; Hackman et al. 2013), and has been interpreted as being due to either latitudinal differential rotation or an azimuthal dynamo wave. The effect is definitely the clearest in the primary component of the binary system II Peg, where a clear linear trend is visible for over ten years, (see, e.g. Lindborg et al. 2013), whereas only disrupted up and downward trends were seen in the single giant FK Com (Hackman et al. 2013). The complexity level and the type of phase behavior seen in FK Com are similar to those we report for LQ Hya; this could be an indication of the binarity of II Peg having an influence on stabilizing the active longitudes in comparison to the single stars, LQ Hya and FK Com, representing this class.
In the case of II Peg and FK Com, however, it was difficult to rule out the differential rotation scenario definitely, as for instance in the case of II Peg, a similar magnitude of drift could have been caused by an antisolar differential rotation profile with k comparable to the values deduced from observations. Also, in the case of LQ Hya, the situation is somewhat similar: the deduced values of the differential rotation parameter k range between 0.002 ... 0.032, the smallest values being obtained from DI, the largest from photometry. The DI and ZDI studies (e.g., Strassmeier et al. 1993; Rice & Strassmeier 1998; Donati 1999; Donati et al. 2003b; Kovári et al. 2004; Cole et al. 2014a) indicate that the majority of the spot activity of LQ Hya occurs at two different latitudinal regions, namely at high and nearly equatorial regions. It cannot be ruled out that in this kind of system, spots could be drifting from the highlatitude location to the lower latitude location, in which case they would gradually attain a faster rotation rate due to the most likely solarlike rotation law of the object. The maximal latitude range of the drifting structure versus the mean spot latitude would be of the order of π/ 4, and therefore the implied k for the spot structures roughly half of the differential rotation parameter, i.e., ... 0.016. The value that can be directly computed from the period difference of the coherent structure and the mean movement of the spots reads (15)which is close to the lower limit of the values derived from DI. Therefore, again, there is certainly enough differential rotation on the object to be the cause of the observed phasetime drift. One must note that this drift would not cause strictly linear (but curved) trends in the phasetime plots (see Pelt et al. 2011, for a simulated example). One should also expect drifts from the lower latitude spot band to the higher one, with opposite direction of the trend in the phasetime plots. These are indeed seen, but with less pronounced phase coherence.
Alekseev (2005) also proposed, based on photometric spot modeling, that a latitudinal dynamo wave, the spot activity migrating from the equator poleward (the solar butterfly reversed), is present on the object. The analysis of Lehtinen et al. (2012) did not reveal these trends, nor does the CF analysis picture of linear down or upward trends support this picture.
The phase behavior, if not due to differential rotation nor latitudinal dynamo waves, could also be a manifestation of an azimuthal dynamo wave, predicted to be excited in rapid rotators (e.g., Krause & Raedler 1980), verified from meanfield dynamo models (e.g., Moss et al. 1995; Küker & Rüdiger 1999; Mantere et al. 2013), and also now found from direct numerical simulations (Cole et al. 2014b). These dynamo waves most often behave as if detached from the overall rotation of the object, moving with a different speed than the stellar surface. Their rotation is rigid even in a differentially rotating object. Therefore, a systematic linear phase drift could most directly be linked to the presence of azimuthal dynamo waves. Dynamo theory, on the other hand, serves no direct explanation as to why the linear trends are broken and reversed, which is clearly the case for LQ Hya. It is, however, well known that the more supercritical a dynamo is, the more chaotic are the solutions (see, e.g., Brandenburg et al. 1989).
7. Conclusions
In this work we have presented an analysis of LQ Hya photometry for 1982−2014. Several different statistical methods were used first to nail down a suitable carrier frequency for our main analysis tool, the CF method. From this preliminary analysis, we learned several interesting points.
Firstly, there is a certain cutoff in the spectrum at the highfrequency end. This can be interpreted as a limiting value for the spot cycle length at the low latitudes or near the equator of the star. Secondly, an interesting feature appearing as the result of the same analysis is the bimodal shape of the frequency spectrum. The explanation for this can be searched from different causes, e.g., latitudinal distribution of the spots. From DI maps for LQ Hya it has become evident that spot regions tend to lie either on high or low latitudes while there seems to be spotless area on midlatitude range (Cole et al. 2014a; Donati et al. 2003b). Other possible explanations might be either radial differential rotation manifesting itself through different anchoring depths for spots or hemispherical asymmetry.
In previous studies the focus has been on searching for active longitudes on the star (Jetsu 1993; Berdyugina et al. 2002; Kovári et al. 2004; Lehtinen et al. 2012). Here we took a different approach by estimating the mean rotation period of the spot structures on the star using the phase dispersion statistic D^{2}(P). This period is a close analog to the Carrington rotation period of the Sun. In subsequent CF analysis, we used the obtained value as a carrier period and produced the corresponding phase plots. We noticed shorter and longer, nearly linear, trends with different slopes reflecting the “duty times” of certain periods during these time frames. Two epochs (1990–1994; 2003–2009) with a downward trend are especially pronounced. The overall picture is inconsistent with the antisolar butterfly diagram postulated by Alekseev (2005). The possible sources of these trends include disrupted azimuthal dynamo waves and solarlike latitudinal differential rotation.
We also tried a multiperiodic model to describe the photometry of LQ Hya. This, however, led to low R^{2} values and the regular structure of the simple model was lost in the phase diagram of the CF. Therefore, we concluded this model to be barely suitable. However, this analysis led to a phase modulation timescale of roughly 6.94 yr, which is close to values earlier derived from the mean brightness variation of the star (Jetsu 1993; Strassmeier et al. 1997; Cutispoto 1998; Oláh et al. 2009). From the global CF, we calculated the phases of the minima and compared them with the results obtained by Lehtinen et al. (2012) using the CPS method. These two models appeared to be in good agreement. Comparison with the results by Cole et al. (2014a) using DI technique was challenging because of the lack of overlap of the corresponding observing seasons. However, around late 1999 and late 2000 rough agreement between the epochs of the photometric minima and the DI spots can be seen.
We did qualitative analysis of the local CF for 27 observing seasons. We detected four flipflop type events from the phase plots, and three of these matched the epochs of flipflops obtained from the global CF. The timing of the events appears to be random, which excludes the possibility of the 5.2 yr cycle reported by Berdyugina et al. (2002). Comparison of the phase plots with those reported by Lehtinen et al. (2012) showed a good agreement; a majority of features can be detected from both analyses.
Acknowledgments
This work has been supported by the Academy of Finland Centre of Excellence ReSoLVE (N.O., M.J.K., and J.P.). The work of T.H. was financed through the project Active Suns by the University of Helsinki. G.W.H. acknowledges longterm support from Tennessee State University and the State of Tennessee through its Centers of Excellence program.
References
 Alekseev, I. Y. 2005, Astrophys., 48, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, J. R., Collier Cameron, A., Donati, J.F., et al. 2005, MNRAS, 357, L1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Barning, F. J. M. 1963, Bull. Astron. Inst. Netherlands, 17, 22 [NASA ADS] [Google Scholar]
 Berdyugina, S. V., Pelt, J., & Tuominen, I. 2002, A&A, 394, 505 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brandenburg, A., Krause, F., Meinel, R., Moss, D., & Tuominen, I. 1989, A&A, 213, 411 [NASA ADS] [Google Scholar]
 Cao, D.T., & Gu, S.H. 2014, AJ, 147, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Cole, E., Hackman, T., Käpylä, M. J., et al. 2014a, A&A, submitted [arXiv:1504.03673] [Google Scholar]
 Cole, E., Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2014b, ApJ, 780, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Covino, S., Panzera, M. R., Tagliaferri, G., & Pallavicini, R. 2001, A&A, 371, 973 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cutispoto, G. 1991, A&AS, 89, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Cutispoto, G. 1998, A&AS, 131, 321 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Donati, J.F. 1999, MNRAS, 302, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Donati, J.F., Collier Cameron, A., & Petit, P. 2003a, MNRAS, 345, 1187 [NASA ADS] [CrossRef] [Google Scholar]
 Donati, J.F., Collier Cameron, A., Semel, M., et al. 2003b, MNRAS, 345, 1145 [NASA ADS] [CrossRef] [Google Scholar]
 Eggen, O. J. 1984, AJ, 89, 1358 [NASA ADS] [CrossRef] [Google Scholar]
 Fekel, F. C., Bopp, B. W., Africano, J. L., et al. 1986, AJ, 92, 1150 [NASA ADS] [CrossRef] [Google Scholar]
 Gunn, A. G., Mitrou, C. K., & Doyle, J. G. 1998, MNRAS, 296, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Hackman, T., Pelt, J., Mantere, M. J., et al. 2013, A&A, 553, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Işık, E., Schmitt, D., & Schüssler, M. 2011, A&A, 528, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jetsu, L. 1993, A&A, 276, 345 [NASA ADS] [Google Scholar]
 Kim, Y.C., & Demarque, P. 1996, ApJ, 457, 340 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Olemskoy, S. V. 2011, MNRAS, 411, 1059 [NASA ADS] [CrossRef] [Google Scholar]
 Kovári, Z., Strassmeier, K. G., Granzer, T., et al. 2004, A&A, 417, 1047 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krause, F., & Raedler, K. H. 1980, Meanfield magnetohydrodynamics and dynamo theory (Oxford: Pergamon Press) 271 [Google Scholar]
 Küker, M., & Rüdiger, G. 1999, A&A, 346, 922 [NASA ADS] [Google Scholar]
 Lehtinen, J., Jetsu, L., Hackman, T., Kajatkari, P., & Henry, G. W. 2011, A&A, 527, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lehtinen, J., Jetsu, L., Hackman, T., Kajatkari, P., & Henry, G. W. 2012, A&A, 542, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lindborg, M., Korpi, M. J., Hackman, T., et al. 2011, A&A, 526, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lindborg, M., Mantere, M. J., Olspert, N., et al. 2013, A&A, 559, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lomb, N. R. 1976, Ap&SS, 39, 447 [NASA ADS] [CrossRef] [Google Scholar]
 Mantere, M. J., Käpylä, P. J., & Pelt, J. 2013, in IAU Symp. 294, eds. A. G. Kosovichev, E. de Gouveia Dal Pino, & Y. Yan, 175 [Google Scholar]
 McIvor, T., Jardine, M., Collier Cameron, A., Wood, K., & Donati, J.F. 2004, MNRAS, 355, 1066 [NASA ADS] [CrossRef] [Google Scholar]
 Messina, S., & Guinan, E. F. 2003, A&A, 409, 1017 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Moss, D., Barker, D. M., Brandenburg, A., & Tuominen, I. 1995, A&A, 294, 155 [NASA ADS] [Google Scholar]
 Oláh, K., Kolláth, Z., & Strassmeier, K. G. 2000, A&A, 356, 643 [NASA ADS] [Google Scholar]
 Oláh, K., Kolláth, Z., Granzer, T., et al. 2009, A&A, 501, 703 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Oláh, K., van DrielGesztelyi, L., & Strassmeier, K. G. 2011, in Proc. International Astronomical Union, Vol. 7, Comparative Magnetic Minima: Characterizing quiet times in the Sun and Stars, 279 [Google Scholar]
 Ossendrijver, A. J. H. 1997, A&A, 323, 151 [NASA ADS] [Google Scholar]
 Pelt, J. 1983, in Statistical Methods in Astronomy, ed. E. J. Rolfe, ESA SP, 201, 37 [Google Scholar]
 Pelt, J., Olspert, N., Mantere, M. J., & Tuominen, I. 2011, A&A, 535, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reinhold, T., & Reiners, A. 2013, A&A, 557, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rice, J. B., & Strassmeier, K. G. 1998, A&A, 336, 972 [NASA ADS] [Google Scholar]
 Saar, S. H., & Brandenburg, A. 1999, ApJ, 524, 295 [NASA ADS] [CrossRef] [Google Scholar]
 Saar, S. H., Piskunov, N. E., & Tuominen, I. 1994, in Cool Stars, Stellar Systems, and the Sun, ed. J.P. Caillault, ASP Conf. Ser., 64, 661 [Google Scholar]
 Scargle, J. D. 1982, ApJ, 263, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Schuessler, M., & Solanki, S. K. 1992, A&A, 264, L13 [NASA ADS] [Google Scholar]
 Stellingwerf, R. F. 1978, ApJ, 224, 953 [NASA ADS] [CrossRef] [Google Scholar]
 Strassmeier, K. G., Rice, J. B., Wehlau, W. H., Hill, G. M., & Matthews, J. M. 1993, A&A, 268, 671 [NASA ADS] [Google Scholar]
 Strassmeier, K. G., Bartus, J., Cutispoto, G., & Rodono, M. 1997, A&AS, 125, 11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tetzlaff, N., Neuhäuser, R., & Hohle, M. M. 2011, MNRAS, 410, 190 [NASA ADS] [CrossRef] [Google Scholar]
 White, R. J., Gabor, J. M., & Hillenbrand, L. A. 2007, AJ, 133, 2524 [NASA ADS] [CrossRef] [Google Scholar]
 You, J. 2007, A&A, 475, 309 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
Fig. 1 Combined data set from the three different sources. The data set referred to as D1 consists of all available data, while data set D2 comprises only T3APT data (Lehtinen + unpublished). 

Open with DEXTER  
In the text 
Fig. 2 Amplitude spectrum for the important period range 1.5−1.7 days (gray) and first 100 frequencies obtained by sequential prewhitening procedure (black). 

Open with DEXTER  
In the text 
Fig. 3 Set of 28 strongest amplitudes (black) and original amplitude spectrum (gray). The horizontal line marks the cutoff level obtained from 100 randomized samples. 

Open with DEXTER  
In the text 
Fig. 4 Folded and stacked light curve model with periods and . This regular structure helps to describe only R^{2} = 23.542% of overall variability. 

Open with DEXTER  
In the text 
Fig. 5 Carrier fit obtained with the carrier obtained from the twoperiodic model of Sect. 4.2. 

Open with DEXTER  
In the text 
Fig. 6 Phase dispersion statistic D^{2}(P) for a range of correlation lengths. The dispersion spectrum becomes strongly asymmetric around 230 days and then splits into set of peaks. 

Open with DEXTER  
In the text 
Fig. 7 Curves of the phase dispersion statistic D^{2}(P) for different correlation lengths in days: 100 (black), 200 (gray), and 300 (light gray). Dashed curve corresponds to the bestfitting Gaussian to the curve with correlation length 100 days. 

Open with DEXTER  
In the text 
Fig. 8 Phase diagram for global fit with carrier period 1 60514, K = 2 and L = 30. On left with normalized amplitudes, on right with actual magnitudes 

Open with DEXTER  
In the text 
Fig. 9 Phase diagrams of local fits with refined carrier period 1 60514. 

Open with DEXTER  
In the text 
Fig. 10 Phases of the minima for period . Red dots: primary minima from CF; red pluses: primary minima from CPS; blue points: secondary minima from CF; blue crosses: secondary minima from CPS; orange circles: minima from DI; and bold green vertical lines denote the epochs of possible flipflop events. The black bars with the numbers correspond to the photometric observing seasons, the numbers on top of the figure mark the observing seasons of DI. 

Open with DEXTER  
In the text 
Fig. 11 Zoomin to the global CF near the flipflop event detected in the segment 2. The data points are drawn as small rectangles. 

Open with DEXTER  
In the text 