Issue 
A&A
Volume 610, February 2018



Article Number  A86  
Number of page(s)  32  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201628870  
Published online  12 March 2018 
Investigating light curve modulation via kernel smoothing
I. Application to 53 fundamental mode and firstovertone Cepheids in the LMC
^{1}
Observatoire de Genève, Université de Genève, Ch. d’Ecogia,
1290
Versoix,
Switzerland
email: maria.suveges@unige.ch;sueveges@mpia.de
^{2}
Department of Physics & Astronomy, The Johns Hopkins University,
3400 N Charles St,
Baltimore,
MD
21218, USA
email: ria@jhu.edu
Received:
7
May
2016
Accepted:
3
July
2017
Context. Recent studies have revealed a hitherto unknown complexity of Cepheid pulsations by discovering irregular modulated variability using photometry, radial velocities, and interferometry.
Aims. We aim to perform a statistically rigorous search and characterization of such phenomena in continuous time, applying it to 53 classical Cepheids from the OGLEIII catalog.
Methods. We have used local kernel regression to search for both period and amplitude modulations simultaneously in continuous time and to investigate their detectability. We determined confidence intervals using parametric and nonparametric bootstrap sampling to estimate significance, and investigated multiperiodicity using a modified prewhitening approach that relies on timedependent light curve parameters.
Results. We find a wide variety of period and amplitude modulations and confirm that first overtone pulsators are less stable than fundamental mode Cepheids. Significant temporal variations in period are more frequently detected than those in amplitude. We find a range of modulation intensities, suggesting that both amplitude and period modulations are ubiquitous among Cepheids. Over the 12year baseline offered by OGLEIII, we find that period changes are often nonlinear, sometimes cyclic, suggesting physical origins beyond secular evolution. Our method detects modulations (period and amplitude) more efficiently than conventional methods that are reliant on certain features in the Fourier spectrum, and prewhitens time series more accurately than using constant light curve parameters, removing spurious secondary peaks effectively.
Conclusions. Period and amplitude modulations appear to be ubiquitous among Cepheids. Current detectability is limited by observational cadence and photometric precision: detection of amplitude modulation below 3 mmag requires spacebased facilities. Recent and ongoing space missions (K2, BRITE, MOST, CoRoT) as well as upcoming ones (TESS, PLATO) will significantly improve detectability of fast modulations, such as cycletocycle variations, by providing highcadence highprecision photometry. Highquality longterm groundbased photometric time series will remain crucial to study longerterm modulations and to disentangle random fluctuations from secular evolution.
Key words: methods: data analysis / methods: statistical / stars: variables: Cepheids / stars: oscillations / Magellanic Clouds
© ESO 2018
1 Introduction
Classical Cepheid variable stars (from hereon: Cepheids) have been the focus of a great deal of research since their discovery by Goodricke (1786), who suggested that their study “may probably lead to some better knowledge of the fixed stars”. Indeed, Cepheids have been of great historical importance for the understanding of stellar evolution and structure. Although Cepheids are usually considered to be highly regular variable stars, Cepheid pulsations were shown very early on to exhibit timedependencies (e.g., Eddington 1919).
In particular the changing periods of Cepheids have received much attention (e.g., Szabados 1983; Berdnikov & Ignatova 2000; Pietrukowicz 2001, 2002; Turner et al. 2006), since they offer an immense opportunity for studying the secular evolution of stars on human timescales (decades) and provide important tests of stellar evolution models (e.g., Fadeyev 2013; Anderson et al. 2016b). In addition, changing periods complicate phasefolding of timeseries data obtained over long temporal baselines and often have to be accounted for when determining the orbit of longperiod binary systems containing Cepheids (e.g., Szabados et al. 2013; Anderson et al. 2015). In addition, cyclic variations of pulsation periods exhibited by some Cepheids have been discussed in terms of the lighttime effect due to orbital motion (Szabados 1989), although only a few cases have been confirmed using radial velocities (Szabados 1991). Period changes may also be related to the muchdiscussed linearity of the periodluminosity relation, see GarcíaVarela et al. (2016) and references therein.
Recent detailed studies of both large samples of Cepheids in the LMC (Poleski 2008) and of individual stars in the Galaxy (e.g., Berdnikov et al. 2000; Kervella et al. 2017) have revealed intricate, possibly periodic period change patterns that are not necessarily consistent with the classical picture of secular evolution. One of the most notorious and intricate cases of period and amplitude variations is that of Polaris, the North Star (Arellano Ferro 1983), which is assumed to be crossing the classical instability strip for the first time (Turner et al. 2006). Polaris’ amplitude seemed to diminish to the point of disappearing, which had been interpreted as the Cepheid leaving the instability strip. However, more recent observations have shown that the pulsation amplitude has increased again, and Polaris remains a puzzle. Another unique wellknown example is that of V473 Lyrae (Burki & Mayor 1980; Burki et al. 1982). Combining many years worth of observations, Molnár et al. (2013) were able to trace this star’s amplitude modulation cycles and determined a modulation period of 1204 d. They discussed these modulations in the context of the Blažko (1907) effect, which is better known in RR Lyrae stars and has seen a boost in research thanks to the Kepler mission (e.g., Kolenberg et al. 2010; Benkő et al. 2011, 2014).
Percy & Kim (2014) presented evidence that some longperiod Cepheids exhibit amplitude changes of up to a few hundredths of a magnitude over timescales of a few hundred or thousands of days, potentially exhibiting cyclic behavior (e.g., for U Carinae). Such strong modulations are presumably not very common among classical Cepheids, or else they would likely be found more frequently in longterm photometric surveys such as the All Sky Automated Survey (Pojmanski 2002) or in other longterm Cepheid photometry (e.g., Berdnikov et al. 2014). Soszynski et al. (2008) mention that about 4% of fundamentalmode (FU) and 28% of firstovertone (FO) Cepheids are Blažko Cepheids, identified via secondary period peaks near the primary period (hereafter “twin peaks”), which are found after prewhitening the light curve using the primary period. Additionally, among the entire set of 3374 Cepheids, eight (all FO) were labeled as having variable amplitude. Soszyński et al. (2015a) have since provided additional targets of interest in this regard. Periodicity of such light curve modulations, if it can be firmly established, may be indicative of Cepheids pulsating in more than one mode (Moskalik & Kołaczkowski 2009). While this paper was under review, Smolec (2017) further reported light curve modulation in 51 fundamentalmode Cepheids located in the Small and Large Magellanic Clouds, none of which overlap with the stars discussed in the present work. Derekas et al. (2017) and Kanev et al. (2015) found a lowamplitude periodic modulation of the pulsation period analyzing the Kepler data of the (fundamentalmode) Cepheid V1154 Cygni. BRITE observations suggest a period modulation also for the fundamentalmode T Vulpeculae (Smolec et al. 2016). Evidence for nonradial modes in Cepheids has been found in a sample of 138 Small Magellanic Cloud (SMC) FO Cepheids that exhibit light curve modulation (Soszynski et al. 2010; Dziembowski 2016; Smolec & Śniegowska 2016).
Though difficult to detect with groundbased observatories, small amplitude light curve fluctuations appear to be rather common, if photometry is sufficiently precise and denselysampled. Derekas et al. (2012) first showed this for V1154 Cygni. Evans et al. (2015) used the MOST satellite to demonstrate the different types of irregularities seen in the fundamentalmode Cepheid RT Aurigae and the firstovertone Cepheid SZ Tauri. Such lowamplitude modulations and period “jitter” may be explained by convection and/or granulation (Neilson & Ignace 2014). Stothers (2009) furthermore proposed a model involving activity cycles to explain period and light amplitude changes in shortperiod Cepheids. However, light curve modulations remain difficult to detect, even with precise spacebased photometry (Poretti et al. 2015).
While most amplitude modulations in Cepheids are found using photometric measurements, the extreme precision afforded by modern planet hunting instruments has recently enabled the discovery of small amplitude spectral modulations in Cepheids (Anderson 2014, 2016). Furthermore, tentative evidence for modulated angular diameter variability in the longperiod Cepheid ℓ Carinae based on longbaseline interferometry has recently been presented by Anderson et al. (2016a).
In summary, recent advances in instrumentation have enabled the discovery that Cepheid variability is not as regular as often assumed. Whether or not irregularities are detected is dominated by observational precision and timesampling. Moreover, there is evidence that not all irregularities of Cepheid variability share the same origin; for instance, the timescales of radial velocity modulation in short and longperiod Cepheids are very different (Anderson 2014).
Given the patchy evidence for irregularities and modulations in Cepheid variability, it is important to characterize how and how often these phenomena occur. To this end, we have implemented a method based on local kernel estimation to detect irregularities in Cepheid pulsations. For the first time, this method allows us to search for smooth variations of light curve amplitudes and periods in continuous time and enables the quantification of the significance of the detected effects. We here describe our technique, which we apply to a total of 53 FU and FO Cepheids from the OGLEIII catalog (Soszynski et al. 2008). In a followup paper, we will then apply this technique to the full sample of OGLEIV classical Cepheids (Soszyński et al. 2015b) to investigate limits of detectability, the rate of occurrence of period and amplitude modulations in Cepheids, and to characterize them. This will be a crucial step toward a physical understanding of these phenomena.
This paper is structured as follows. We describe our method for analyzing light curves in Sect. 2, which is divided into target selection (Sect. 2.1) and a brief description of the slidingwindows based light curve modeling (Sect. 2.2). We present the results of this modeling in Sect. 3, which we divide into subsections dedicated to changing periods (Sect. 3.1), changing amplitudes (Sect. 3.2), and a discussion of how light curve shapes change with time (Sect. 3.3). We present the implications of the new method on results of a multiperiodicity analysis and on prewhitening artefacts (Sect. 3.4), compare the trends and fluctuations discovered among different groups of Cepheids (Sect. 4.1), investigate their relationships with physical parameters of the Cepheids (Sect. 4.2), and compare our results with the literature. We summarize and conclude in Sect. 5. We explain the statistical methodology in detail in Appendix A. Using simulations of modulated periodicity with parameters taken from real Cepheids, we benchmark the detectability of various modulation types and the performance of the kernel method in Appendices B and C. Figures illustrating the results for all 53 Cepheids and tables containing the numerical results from the fitting procedures are given in in Appendices D and E.
2 Data and methods
2.1 Data and target selection
We analyzed Iband photometric timeseries data (light curves) of classical Cepheids in the LMC published by the second and thirdgeneration Optical Graviational Lensing Experiment (Udalski et al. 1999; Soszynski et al. 2008, OGLEII and III). These data were taken with the 1.3 m Warsaw telescope at Las Campanas, Chile and reduced using difference imaging (Alard & Lupton 1998; Alard 2000; Wozniak 2000; Udalski et al. 2008). OGLE photometry is particularly wellsuited for the study of modulations in Cepheid variability due to its high quality (precision of up to 0.02 mag), long temporal baseline (spanning up to 12 yr), large number of observations (up to ≳ 1500 per target), excellent homogeneity, and large number of Cepheids available for study. The time sampling of the light curves from the two surveys is seasonal, with the longest gap (between the end of OGLEII and the start of OGLEIII) up to 300 dfor some stars. The instrumentation changed between the two phases of the survey, around HJD − 2 450 000 ≈ 2000. Iband photometry in the LMC from the third phase was carefully calibrated to seamlessly fit with the photometric data from the second based on more than 620 000 stars in 78 overlapping subfields (Udalski et al. 2008). All data were obtained from the OGLEIII server^{1}.
This work has the dual goals of assessing the performance of our methodology and of investigating the phenomenology of the modulated variability exhibited by Cepheids. To this end, we selected a sample of 53 Cepheids consisting both of ones likely to exhibit modulations and ones likely to be stable pulsators.
Whether or not a given Cepheid is likely to show the effect is difficult to determine a priori. A sign of modulations may be the presence of peaks in the secondary periodogram at a frequency very close to the primary one (a “twin” of the primary peak). The reason is that if the harmonic decomposition of the oscillation or its period varies with time, then prewhitening with a constant model will not lead to a perfect removal of the primary pulsation from the light curve, and the residual signal would appear in the secondary periodogram as a twin peak. These twin peaks may or may not be resolved by the Rayleigh criterion, depending on various factors such as the type of the modulation (trendlike, periodic or irregular), its typical timescale, its amplitude, and the distribution of the observational times. In our sample, the absolute difference between the twin peaks f_{1} − f_{2} was between 0.00017 and 0.00076 d^{−1} except for one Cepheid, CEP1564, for which this was 0.0176 d^{−1}). Among our targets, the frequency separations between primary and twin peak were scattered around the Rayleigh limit (approximately 0.00025 d^{−1} for our data), so we had both resolved and unresolved twin peaks.
We thus selected as likely irregular candidates those Cepheids for which a multiperiodicity analysis reveals a twin peak after prewhitening. We refer to these targets here as “twinpeak Cepheids” rather than adopting the terminology of Soszynski et al. (2008) who refer to them as “Blažko Cepheids”, since it is still unclear whether the origin of the irregularities of the Cepheid variability is the same as that of the amplitude modulation found in RR Lyrae stars. In order to assess the efficacy of the twin peaks phenomenon as a diagnostic for identifying pulsation irregularities and to provide a baseline for comparison, we also selected likely regular, or “control”, sample consisting of target stars that do not exhibit twin peaks. Since firstovertone (FO) Cepheids are considered to be more irregular (less stable) than fundamentalmode (FU) Cepheids, we treated these two groups separately. To create a basis for target selection, we first modeled all individual light curves of fundamental (FU) and firstovertone (FO) Cepheids that consisted of more 700 observations, and inspected their residual periodograms. For the purpose of sample selection, we modeled light curves using a nonperiodic fifthorder polynomial trend (to account for a possible temporal evolution of the instrument zeropoint as well as other spurious or true changes in mean magnitude) and a Fourier series with ten harmonics using the period from the OGLEIII catalog. Secondary (residual) periodograms were computed for each star using the method of Zechmeister & Kürster (2009). Including the polynomial trends proved efficient at removing artefacts from the residual periodograms, such as high peaks near 0, 1, 2, … d^{−1} that otherwise dominated. Although a precise light curve modeling could have called for more or less than ten harmonics, we found this a sufficient approach for sample selection. In all later stages of the analysis, in particular during the sliding window analysis, we used a more detailed modeling with different harmonic orders, which were determined separately for each Cepheid (cf. Sect. 2.2 and Appendix A).
Figure 1 shows examples of each of the twin peaks and control Cepheid groups. The left hand side exemplifies the twin peaks case with (OGLELMC)CEP1998, which has a high secondary peak in the residual periodogram. For this star, the appearance of the twin peak occurs together with a prominent separation of the folded residual light curves according to observation times: the early observations trace a different line from later ones. Such a visual separation is not observed in every twinpeaks Cepheid. Instead, various degrees of this pattern can be found, although its prominence does seem to correlate with the height of the peak in the residual periodogram. The fact that the earliest and the latest observations are close together, while the reddish dots of midsurvey observations are far from them, indicates a variation nonlinear in time, and thus excludes a misestimated period from the possible explanations. For comparison, the right hand panel shows CEP1543, for which the residuals are flat and no significant secondary peaks are found in the residual periodogram.
Since our method aims to smoothly and continuously trace temporal variations of pulsation periods and amplitudes, it is necessary for all objects studied to have a large number of sufficiently densely and uniformly distributed observations over the survey time span. We thus limited possible targets to those with at least 700 Iband observations, distributed in a way that ensured none of the time intervals used in our sliding windows fits contained less than 90 observations.
Of all the Cepheids that satisfied these conditions, we selected a sample of 12 FU and 12 FO Cepheids that exhibit twin peaks. We further included five additional FO Cepheids with clearly visible changing amplitudes, which we encountered during light curve inspections. These were included despite having a slightly lower limit of 70 observations for each data subinterval. This was deemed acceptable, since FO Cepheids have nearly sinusoidal light curve shapes that require fewer harmonics for an adequate light curve representation. All of these amplitudechanging Cepheids exhibit twin peaks, and therefore will be treated as part of the twinpeak FO group. Similarly, we selected 12 FU and 12 FO Cepheids for which prewhitening did not reveal twin peaks as the control samples.
2.2 Detecting period and amplitude modulation using sliding windows
Studying timedependent variability phenomena has been gaining traction for some time. A sotermed timedependent Fourier analysis has been applied to nonlinear pulsation models of RR Lyrae stars (Kovacs et al. 1987), RRc stars observed by the Kepler spacecraft (Moskalik et al. 2015), RR Lyrae stars in the Galactic Bulge observed by OGLE (Netzel et al. 2015), and 138 FO Cepheids in the SMC (Smolec & Śniegowska 2016). As an alternative, the analytic signal processing method has also been applied to hydrodynamical models (Kolláth et al. 2002) and to investigate the period doubling phenomenon in RR Lyrae stars using Kepler data (Szabó et al. 2010). Given the sensitivity of the method to (seasonal) gaps in the groundbased OGLE data, analytic signal processing is not a suitable choice for the present investigation.
Since previous studies of timedependent Fourier and analytic signal processing techniques adopted fixed oscillation frequencies, any fluctuations in pulsation period were absorbed as phaseshifts in the Fourier phase coefficient (Moskalik et al. 2015; Szabó et al. 2010). Here, we sought to develop a method capable of efficiently dealing with data gaps that simultaneously determines changes in Fourier coefficients and pulsation period. Moreover, our method makes no assumptions as to the periodicity or repeatability of any detected modulations.
We have adopted a highly flexible model to describe the potentially diverse and presumably small variations of Cepheid light curves^{2}. Physical causes for changes in period, for instance, may originate from secular evolution or binarity (lighttime effect). However, amplitude modulations are not so easily explained and modeled. In our sample, the separation of the residual light curves shown in the left panel of Fig. 1 suggests nonlinearity of the changes. Moreover, the diversity of twinpeak structures in the examined secondary periodograms, though the irregular OGLE time sampling and the photometric noise undoubtedly affect the shapes, also suggests that modulation patterns may not be strictly repetitive. Figure 2 gives a few examples of this diversity. While the present work was under review, Smolec (2017) presented an investigation of periodic light curve modulation in SMC and LMC Cepheids based on a systematic search for double modulation side peaks in periodograms following a standard prewhitening technique. Such an approach assumes periodicity of any detected modulations as well as only mild effects of the uneven time sampling and noise. Our visualizations of the residual light curves, with the observing time colorcoded, suggest that many kinds of secondary peak structures can be associated with visible modulations of the period and/or the harmonic parameters, and these are not necessarily strictly periodic or linear. Our goal is therefore to be open to all possibilities, and allow for nonlinear, cyclic, and noncyclic components in the modeling of the modulations. Local kernelmodeling (sliding window technique; e.g., Fan & Gijbels 1996) is a simple option that provides this flexibility.
We defined these sliding windows using a grid of times that covers the entire observational baseline. The first grid point τ_{1} is fixed to the time of the first observation of the star, and every subsequent grid point is defined as τ_{i} = τ_{1} + (i − 1) × 30 days. For most of our Cepheids, 137 or 128 grid points were thus defined. One target (CEP2580) was limited to 80, since observations of this target started roughly five years after the others. We wish to obtain a local estimate of pulsation period and harmonic amplitudes at each gridpoint τ_{i}. Therefore, at each τ_{i}, we selected the data within a threeyear window centered at that point, and fit a harmonic model with a thirdorder polynomial trend: (1)
where are assumed to be independent Gaussian errors. The harmonic order M, individually selected for each Cepheid, is based on the constant models using the whole time span: M = M_{0} + 2, where M_{0} is the order of the best constant model for all data by the Bayes Information Criterion (Schwarz 1978), and the two extra terms are added to allow for variations in the light curve shape. M was subsequently kept fixed at all gridpoints, though all the model parameters were refitted in each window. Thus, we obtain bestfit pulsation period, harmonic amplitudes, and polynomial coefficients at each gridpoint. Any changes in mean brightness are thus absorbed by the polynomial term in Eq. (1) and treated as nuisance parameters.
The parameter estimation is performed by a weighted nonlinear least squares procedure that optimizes both the harmonic parameters and the pulsation period separately in each window. The result is a time series of the bestfit model parameters θ(τ_{i} ) with corresponding pointwise confidence bands, where θ(τ_{i}) can stand for any of the pulsation period, harmonic parameters, or peaktopeak amplitude. We estimated both approximate theoretical pointwise confidence bands and ones based on a Monte Carlo experiment; the two were, in general, very close to each other.
The weighting scheme in this estimation procedure is chosen in a special way, to put emphasis on the process near the central gridpoint. To decrease bias there, we combined the usual error weighting with the kernel weights of local modeling: inverse variances (as given by the square of the photometric uncertainty) were multiplied with a factor derived from a normal density centered at τ_{i} with a standard deviation of 182.5 days (at the ends of the survey timespan, the definition remains the same; we did not lengthen the series by the addition of artificial points based on some rule of extrapolation from the observed values). Doing so attributes a higher influence to observations made closer to τ_{i} (important not to oversmooth, and to trace irregularities better), while also weighting data according to their reliability. The result is a weighting scheme that increases the impact of the most relevant and most reliable observations, while it reduces the high correlation observed between estimates in neighboring windows when using simple error weighting.
In addition to its flexibility, the local sliding window model has a few more advantages over complex global models with constant parameters. Since it is local, abrupt shifts in mean magnitude due to calibration effects between OGLEII and III data affect it only in windows including the time of the shift, and thanks to the weighting scheme, the effect is attenuated in windows centered relatively far from this time. Moreover, the thirdorder polynomial component in model (1) accounts for any further local variation in the mean magnitude that could bias the results of frequency analysis. In the absence of information in the sampling gaps, the estimates may of course be biased or have a higher than average variance.
The finite window size has two main effects on the estimates. First, the detectable timescales of period and amplitude modulations of the Cepheids are determined by the 12yr total survey timespan and the window size. Periods longer than approximately the full timespan cannot be distinguished from trends, so approximately 12 yr is the longperiod limit of modulation cycles we can identify. In the shortperiod limit, any fluctuations on timescales shorter than about two years are smoothed out due to our use of sliding windows giving high weight to observations within an approximately twoyear time interval. Second, the faster the modulation, the more downward biased (underestimated) the estimated modulation amplitude. This bias is due to the local estimate being a kind of average value within the temporal sliding window. It is proportional to the characteristic frequency of the modulation and can be estimated in a simple way if this latter is known or estimated. We discuss this and provide empirical bias correction formulae in Appendix C.
Bias is also present at the start and end of the observation period, if the fitted parameters are not approximately constant there, since we lack information in the subinterval of the window that stretches over the end. This bias is strongest at the endpoints, decreases as the window includes more data, and vanishes at 1.5 yr from the ends (the halfwidth of the window). We investigated the effect of the sampling gaps and the finite observation timespan using simulations, and found this to be a minor effect, although sampling gaps can indeed cause some systematic distortions in the estimates of an underlying trendlike and oscillatory pattern. Moreover, data gaps do not lead to false detections in the absence of modulation (see Appendix B for a detailed discussion).
To assess whether a constant model sufficiently explains the observed photometric time series, we repeated our sliding windows analysis on a simulated, perfectly repetitive stable reference model using the bestfit constant parameters from a global light curve modeling. Confidence bands were added to this curve by applying a nonparametric Monte Carlo resampling based on the residuals (cf. Appendix A).The estimated functions θ(τ_{i} ) were then compared to the obtained confidence bands. In order to assess significance of departures from the constant parameter values, we employed the multiple hypothesis testing procedure of Benjamini & Yekutieli (2001) to avoid spurious detections due to random fluctuations amplified by strong correlation between neighboring windows.
The fitted magnitudes from the sliding window method enabled an improved prewhitening. We approximated the noisefree magnitude of the star at time t_{i} by a weighted average of the fitted value from the models at the two closest windows, with the weights based on the differences between t_{i} and the window centers. We computed the improved residuals by subtracting these fitted values from the observed magnitudes. Then, we performed a secondary period search (Zechmeister & Kürster 2009) on these improved residuals, and we checked whether twin peaks were still apparent, or if other, weak secondary modes appeared.
In summary, the local kernel modeling method presented here does not impose assumptions on the type of period changes encountered and simultaneously traces temporal variations in pulsation period and light curve shape. This is an improvement over the O−C technique, which does not allow us to fully disentangle fluctuations due to period and Fourier amplitude changes. While the use of sliding windows represents a limitation on the time resolution, it provides tools to attribute statistical significance to the timevariation of signals, and to obtain a coherent, continuoustime picture of the simultaneous variations of the period and Fourier composition.
A detailed description of the fitted local model, the model selection, the stable reference model, the error analysis and the assessment of significance taking into account the correlation caused by the overlap of the windows are given in Appendix A. We complement this with realistic simulated examples illustrating the detection power of the model and the limitations imposed by the OGLE time sampling and the kernel size in Appendix B, using real OGLE observation times and light curve profiles, trends and fluctuations similar to those found in our Cepheid sample. Appendix C considers the bias of the method and provides a means of correcting for it. The analysis presented in the paper was performed using the statistical computing environment R (R Core Team 2015).
Fig. 1
Left: CEP1998 as an example of a twinpeak target, for which secondary periods are found after prewhitening the light curve using the primary period. Right: CEP1543 as an example from the control group (no twin peaks). Colors trace observation date, increasing from yellow to blue. Top panels: phasefolded (with primary period) OGLEIII Iband light curve. Center panels: residuals after subtracting Fourier series model. Bottom panels: periodogram of the residuals shown in the center panel. Pink vertical lines at 0, 1, … d^{−1} indicate parasite frequencies due to residual trends in the time series. The solid light blue vertical line indicates the OGLE catalog (primary) frequency. The dashed light blue lines correspond to its ± 1, 2, … aliases. 
Fig. 2
Enlarged environment of the removed OGLE catalog frequency and the secondary peak in the secondary periodogram for four twinpeak Cepheids. The OGLE catalog frequency is indicated by a vertical blue line. 
Basic parameters and their variation for the twinpeak and amplitudechanging Cepheid sample.
3 Results
3.1 Period changes
Tables 1 and 2 present the results of applying the kernel smoothing technique to the 53 selected Cepheids. Columns DP and T_{P} contain the maximal span of these changes, that is, max P(t) − min P(t), and the length of these intervals in percentage of the total observation span. Figures D.1 and D.2 provide a visualization of the variations in pulsation period compared to a stable reference Cepheid simulated using the best stable parameter estimates.
The figures suggest a continuous broad range of possible variations potentially extending to many or even all Cepheids, rather than separate subclasses of steady and nonsteady pulsators. Within the timescale limitations mentioned at the end of Sect. 2.2 and discussed in Appendix B, the range of the variations may include slow irregular or nearlinear trends, stochastic or multiscale oscillations, quasiregular changes, and any combinations of these. Figure 3 presents an example of each (CEP2132: nearlinear trend; CEP1621: irregular fluctuations, CEP1140: quasiperiodic changes; CEP1833: combination of damped quasiperiodic changes and nearlinear trend). Appendix B suggests that though the OGLE time sampling can cause systematic distortions on the shape of the estimate of a real modulation pattern (see Appendix B), it cannot cause artefacts similar to the observed phenomena. Neither the trends of CEP2132 and 1833, nor the large fluctuations of CEP1621 can be fully explained in this way. Considering the simulation results in Appendices B and C, the nonsignificant quasiperiodic modulation in CEP1140 is more likely to be attributable to a fast oscillation around or above the upper frequency detection limit of our window than to noise or effect of sampling gaps (the weak quadraticlike trend may be the consequence of the time sampling or possibly end effects, though the latter are more likely to cause estimates leveling out than a sharp increase or decrease as here). The exceptionally large scatter of the estimates of the modulation frequency and amplitude presented in Appendix C in the case of CEP1833 warrants prudence in accepting the shape of its modulation, but the existence of a trend seems to be certain and some additional instabilities very likely.
The wide variety of different combinations of these relatively pure types can be appreciated in Figs. D.1 and D.2 among our Cepheid sample. The intervals of deviation from a stable model are much more frequent and last much longer among twinpeak Cepheids (Fig. D.1) than in the control group (Fig. D.2), as emphasized by the extent of grayshaded areas on the figures, and shown by T_{P} in Tables 1 and 2. This is so for both FU and FO Cepheids.
Among the variety of modulation types found (trendlike, stochastic, oscillationlike or arbitrary combination), visual inspection suggests a relatively strong trendlike component for CEP1521, CEP1527, CEP1704, CEP2217 (all FOs), and CEP2132 (FU). Similar trends are visually less obvious, but potentially present for stars CEP1405 (FO), CEP1833 and CEP2470 (both FUs). For most stars, this trend appears nonlinear, or has added fluctuating components (stochastic or oscillationlike); almost pure nearlinear patterns are shown only by CEP2132 and CEP1521, the latter nevertheless with some increasing oscillations toward the end of the observation period (where end effects may also be present). The strongest trend is observed for CEP1833: its frequency change within the decadelong OGLEIII timespan is about 0.008 c/d. Given the several different types of period changes seen here, however, the observational baseline may not be sufficient to ascertain that these period changes are truly caused by secular evolution (cf. also Soszynski et al. 2008): trends observed on such short timescales may actually prove to be a portion of fluctuations with a long characteristic timescale. The frequent presence of relatively strong fluctuations on various timescales further strengthens this impression.
For four twin peak Cepheids (FU CEP1140, FO CEP1536, CEP1564 and CEP1693), the multiple testing procedure did not find significant period changes (T_{P} = 0 in Table 1), and thus, instability of the pulsation period is not confirmed in these stars. CEP1418 also has only a very short interval of period deviation. A strong instability on a shorter timescale than the sliding window length could cause remaining periodicity after removal of the (average) primary frequency. Although the frequency separation between the primary and the secondary periodogram peaks, which is commonly used as an indicator of the modulation’s typical timescale, does not suggest a highfrequency periodic modulation, the sliding window estimates in Fig. D.1 suggest indeed fast (lowlevel) variations for these stars. The strength of such fast variations are strongly underestimated by our threeyear window (see Appendix B), so it is possible that these stars do have a highamplitude, fast period oscillation producing perceivable traces in the secondary periodogram. In addition, CEP1140 and CEP1536 have long intervals when their amplitude deviates significantly from the mean value, which offers an alternative explanation. A look at the period changes of CEP1693 and CEP1418 in Fig. D.1 reveals a slow, weak nonlinear trend (together with some comparatively strong oscillations) in their pulsation period. Although this is not significant with the small number of observations used per window and with our very conservative error assessment procedure, it may be sufficient to give rise to a twin peak when trying to prewhiten by a constant model fitted to all data. For the fifth star, CEP1564, we find no significant deviations from stable pulsations using our multiple hypothesis testing procedure.
Within the control group, our results suggest identical variation types except for trends discernible by eye, and on average faster quasiperiodic or stochastic changes. The variations appear generally milder than those of the twinpeak group, as shown by the values of DP in Tables 1 and 2. The residualbased significance assessment finds these statistically nonsignificant, so we cannot exclude a noise origin of these modulations. However, as mentioned above, variations on timescales comparable to or shorter than the window length are generally underestimated by the kernel method. Fast and strong cyclic variations can cause overdispersion in the residuals at all phases in the secondary light curve, and the underfitted variation causes a strong general overdispersion of the residuals. As a consequence, our conservative residualbased significance assessment procedure (cf. Sect. 2.2 and Appendix A) yields a very broad confidence band around the P_{cat} value of the constant model, and thus, we find the modulation to be nonsignificant. It follows that the twinpeak phenomenon, though it seems to be a good indicator of changes that are trendlike on the survey timespan, can fail to indicate even strong oscillatory modulationsin the pulsation period, and thus miss many potentially scientifically interesting cases. Two such cases merit being mentioned, that of CEP0727 (FU) and CEP1638 (FO). Their secondary periodogram does not show a twin peak, however, their DP values in Table 2 are 0.002 d and 0.0017 d, respectively, both among the highest in all our sample.
Basic parameters and their variation for the control sample.
Fig. 3
Examples for the basic types of period modulation found among the Cepheids. The plots show the kernelestimated pulsation period (in days) versus Julian Date (in days), plotted as a solid thick black line, together with its bootstrapped pointwise CI (thin black lines; see Sect. 2.2). The heavy orange line denotes the catalog period, which was used as known (nonoptimized) value in the fitted stable reference model for the Cepheid. The orange band indicates the nonparametric bootstrap CI around this period, obtained from the procedure described in Sect. 2.2. The dotted black horizontal and vertical lines are aids to the eye to estimate the extent and time interval of the changes. The gray background highlights time intervals where the deviation from the stable reference model was found significant by a multiple testing procedure (Benjamini & Yekutieli 2001). Leftmost panel: nearlinear trend, middle left: stochastic fluctuations, middle right: quasiperiodic instabilities, rightmost: combination of a trend and weakening quasiperiodic changes. 
3.2 Amplitude changes
We find much fewer stars to exhibit significant variations in the amplitude of the leading harmonic term, (cf. Eq. (1)), than in period. Figures D.3 and D.4 present the estimated A_{1} (t) curves for all stars. Tables 1 and 2 summarize the maximum span DA_{1} = max A_{1}(t) − min A_{1}(t) of the first harmonic amplitude A_{1}(t) from the kernel fits, and the extent of the time intervals of significant deviations from the average A_{1} in their final columns DA_{1} and . It appears indeed that the variations in the pulsation period are easier to detect, and for this reason they have been discussed in the literature much more extensively. Nevertheless, we find significant, more or less cyclic, amplitude variations for the twin peak FO Cepheids CEP1527, CEP1535 and CEP1536, in addition to those five Cepheids included because of their clearly visible strong amplitude changes. Several other FO Cepheids (CEP1405, 1561 and 1605) exhibit shorter, erratic excursions from otherwise fairly stable pulsation amplitudes. With the exception of CEP1536, these stars also show changes in pulsation period. There are two FU stars as well with significant amplitude changes, CEP1748 and CEP1140, but contrary to the majority of the FO sample with amplitude modulations, these stars do not undergo significant period changes, as discussed in Sect. 3.1. Our method recovers all (3) Cepheids identified ashaving variable amplitudes in the OGLE catalog, and adds several more examples to this interesting class.
Longterm, trendlike amplitude changes (within the detectability limits of the 12yr survey) seem to be very rare. Two stars that show a pattern compatible with it are the FO Cepheids CEP1561 and 1693. However, even for those, slow fluctuations around a mean amplitude which is stable on the longterm cannot be excluded. One FO Cepheid, CEP2217, has a significantly higher first harmonic amplitude with the sliding window estimation than with the constant reference model. This is due to the fact that its strong trendlike period variation causes large phase shifts of temporally distant observations, and smears the light curve foldedwith a constant period.
Our results suggest that small amplitude variations may be present for many more Cepheids, although these amplitude fluctuations tend to stay below the 95% confidence level with OGLEIII time cadences and our conservative significance assessment procedure. Taking this into account, we interpret our results as an indication that amplitude fluctuations on millimagnitude level and below are a common, possibly ubiquitous, phenomenon whose detectability is currently limited by the availability of sufficient time resolution and photometric precision.
Fig. 4
Four examples of the effect of the sliding window smoothing on residual periodograms. The GLS periodogram of the residuals from the stable reference model is plotted in black. Twin peaks appear very close to the gray background vertical line that indicates the primary pulsation frequency. Daily aliases of the primary frequency are shown as dashed vertical gray lines. We superposed the periodogram of the residuals from the sliding window fit in red. 
3.3 Changes in light curve shape
The morphology of light curves is most commonly described using relative amplitudes and phases, Φ_{j1} and R_{j1} for the leading few harmonics, most importantly for j = 2, 3 (Simon & Lee 1981). However, the continuoustime changes of these parameters are too weak to consider using these quantities, since the variations in the second and third harmonic amplitudes are below the detection limit.
Nevertheless, it is possible to visually compare light curve shapes at different epochs, for example where the peaktopeak amplitudes of the Cepheids are very different. Figures D.5–D.9 show the light curves of the Cepheids in our sample at two such epochs selected individually for each star, such that furthermore the windows do not overlap. We plot the folded light curves so that maximum brightness occurs at phase 0.6 to facilitate visual comparison. Confidence bands are added based on bootstrap repetitions generated with resampled errors superposed to the reconstructed estimated local light curve.
The fundamentalmode Cepheids in our sample, both twinpeak and control, exhibit little scatter in their light curves. In quite a few cases, there are discernible but tinylooking differences in the minimum or maximum brightness or in the pattern of the brightening branch (in a few cases, these can be due to outliers, for example CEP2215 in Fig. D.8, or to unfortunate phase gaps in the data such as for CEP1932 in Fig. D.5). There are a few unusual cases, like that of CEP1833, which has many downward scattered observations in its more recent window, but none in its earlier window; or those of CEP1418 (Fig. D.5), CEP1543, CEP1711, and CEP2774 (Fig. D.8), which seem to have stronger overdispersion than other fundamentalmode Cepheids. These stars have smaller amplitude (< 0.3 mag) than those with less noisy light curve (0.4 mag or above), so this can also be due to the relatively smaller magnitude span of the plots, or can indicate possible further variations, not adequately modeled by the threeyear sliding window estimate.
The variations exhibit a larger range for the overtone stars. As well as a few cases that show light curve variations comparable to the fundamentalmode sample (often those with relatively high average amplitude), and CEP1536 for which our procedure selected apparently an unnecessarily high harmonic order resulting in a wiggly fit (cf. Sect. 2.2), there are many that show high dispersion of observed magnitudes together with strong light curve shape variations. The dispersion of the observationsaffects the quality of the estimated light curve shapes, as is shown by the broad confidence bands around the estimates. In many cases, the observations with high residuals are far off from the center of the window in real observation time. This suggests that relatively strong changes might occur on a shorter timescale than the kernel length.
3.4 Residual periodograms
Our study presented so far led to the conclusion that for the Cepheids in our sample, the presence of a “twin peak" in the secondary periodogram after prewhitening is related to period and/or light curve shape modulations. Prewhitening a light curve that intrinsically contains such modulations with a model of constant period and harmonic amplitudes will be obviously only approximate. Using instead the local estimates of period and Fourier amplitudes yields a more precise magnitude estimate at every observation time, and hence helps to remove the remnants of the main oscillation (the twin peak) from the secondary periodogram. Moreover, a secondary period search on the residuals of a local prewhitening can be expected to detect weak secondary modes more efficiently and more precisely than a constant model. We thus constructed residual periodograms for all Cepheids in the sample by subtracting the timedependent bestfit models.
We compared the two types of residual periodograms, once computed based on the stable reference fit with constant parameters over the full OGLEIII timespan, and once computed using the sliding window fit. Four examples are shown in Fig. 4. In all but two of the 29 twinpeak cases, the sliding window fit removed perfectly the twin peak from the residual periodogram. Moreover, no new spurious peaks were found in the residual periodograms of the 24 control Cepheids, indicating that the method does not introduce any artefacts into the results.
Residual periodograms obtained by timedependent prewhitening are generally uniform and flat for twinpeak FU, control FU, and control FO Cepheids. The top left panel of Fig. 4 shows the secondary periodograms after constant prewhitening (black) and after prewhitening with local estimates (red) in such a case, for the FU twinpeak Cepheid CEP2191. In some cases, a slight, albeit most likely nonsignificant maximum appears at approximately 0, 1, 2, …c∕d. It cannot be decided whether these indicate a parasite frequency of terrestrial origin or a remaining mean trend which is detected by frequency analysis as a lowfrequency signal.
Twinpeak FO Cepheids exhibit more complex patterns in their residual periodograms. Only five of 17 show flat residual periodogram following prewhitening with our timedependent bestfit models. Four of these belong to the group of five FO Cepheids that we selected for their obviously changing amplitudes (CEP0916, CEP1275, CEP1955, and CEP2820), while CEP1536 was an FO twinpeak group member. CEP1119, the 5th Cepheid selected for obvious amplitude changes, still exhibits a twin peak after prewhitening (this can be seen in the top right panel in Fig. 4). This might be caused by imperfect fits in some windows: the sliding window period estimates in the bottom row of Fig. D.1 suggest very sharp period changes that were not precisely fitted. Three other overtone Cepheids exhibit mostly flat residuals similar to those of FU Cepheids with (likely spurious) peaks near 0, 1, 2, …c∕d. For two of them, this peak is weak, for the third (CEP1535, one of the newly discovered amplitudechanging Cepheids), it is strong.
CEP1564, shown in the bottom right panel of Fig. 4, represents an exception. This is the only FO twinpeak Cepheid in our analysis for which we find no significant variability of pulsation period (P_{1} ~ 2.063 d) and A_{1} (A_{1} ~ 0.074 mag), despite a very prominent secondary peak at P_{twinpeak} = 2.141 d after prewhitening with a constant model (see also the bottom right panel Fig. 2). Prewhitening with a timedependent model yields a secondary periodogram and a secondary frequency almost identical to the former one. Folding the residuals with the secondary period yields a scattered, diffuse sinusoidal light curve with a peaktopeak amplitude of about 7.5 mag shown in the bottom panel of Fig. 5. The residuals folded by the primary pulsation frequency lack any trace of the pattern which is typical of other twinpeak Cepheids (compare the upper panel of Fig. 5to the middle panel of Fig. 1). Moreover, the kernel fits using a broad search interval for the local pulsation frequency failed to find any strong modulation, either trendlike or fluctuating. The fits shown in Fig. D.1 remained stable for a wide range of search intervals (including also the frequency of the twin peak). Possible explanations for the persistence of the secondary peak can be that it results from fast fluctuations around or above the detectability limit of our kernel, or that both frequencies are physical, meaning that the star is blended or is a physical binary with another variable star, or that this star is indeed pulsating in two closeby modes.
The residual periodograms of seven remaining FO twinpeak Cepheids exhibit weak secondary peaks similar to those presented in Fig. 4, in the bottom left panel. Their strength ranges from probably nonsignificant to probably significant. All but two of them are at lower frequencies, that is, longer periods than the primary pulsation frequency, and the ratio of the secondary to the primary period P_{2} ∕P_{1} ∈ [1.11, 1.36] (a weak fundamental mode would be between 1.35 and 1.4). For the two at shorter period, P_{2} ∕P_{1} = 0.86 and 0.38. The second of these may be aliased; for the frequency corresponding to the alias, the ratio is 1.46. Only four LMC Cepheids exhibiting secondary periods longer than the primary are mentioned in the OGLEIII catalog, all having primary periods longer than three days. Neither is a member of our sample. Three of the four have P_{2} ∕P_{1} ∈ [1.12, 1.26], while the fourth has P_{2}∕P_{1} = 1.46, similar to our stars. Soszynski et al. (2008) identified a relatively frequent, probably nonradial mode with period ratios around 0.6 in the Magellanic Cloud classical Cepheids; Soszyński et al. (2015a) reported 206 FO stars exhibiting such secondary frequencies among a total of 3530 FO (5.8%) Cepheids in both Magellanic clouds (of which 82 in the LMC), based on joined OGLEIII and OGLEIV data. Despite the relatively frequent occurrence of this mode, none of our stars shows this. In conclusion, our study finds seven out of 29 FO Cepheids (24%) to show some indication of secondary peaks and spanning a broad interval of frequency ratios, which suggests that secondary frequencies, though weak, may be even more common than previously thought in overtone pulsators.
Fig. 5
Top panel: residuals of CEP1564 from the constant model, phasefolded with the average primary period. Bottom panel: residuals of CEP1564 from the kernel fits, phasefolded with the secondary period. The typical error bar is shown in the top right corner. 
4 Discussion
4.1 Separation of trends and oscillatory terms
4.1.1 Modeling timedependent light curve parameters
The complex variations seen in the figures of Appendix D do not suggest a simple statistical model. Aiming in this paper only at a rough separation and quantification of trendlike and fluctuationlike components, we modeled the time dependence of the pulsation period, amplitude of the first harmonic term and the peaktopeak amplitude with a linear model consisting ofa trend and an oscillatory component. This is purely heuristic and intentionally avoids interpretation of period changes as evidence for secular evolution (applies only to linear trends), the lighttime effect (cyclic changes) or other physical origin. This is also warranted since changes in amplitude have no clear theoretical explanation. The benefit of using this model is its simplicity and ability to roughly capture the characteristic size of longterm trends and shortterm fluctuations.A physical interpretation of the observed modulations can later be based on the variations described heuristically.
where θ(t_{i}) can represent timedependent pulsation periods P or first harmonic amplitudes A_{1} at time t_{i}. f_{θ} is the frequency of the oscillatory term, and the error η_{θ,i} is assumed to follow a normal distribution with the locally estimated error σ_{θ,i} on the parameter estimate θ(t_{i}). The strong correlations introduced by the overlapping windows make it necessary to include an autoregressive structure between the consecutive estimates, represented by the correlation coefficient ρ. The model is fitted for a given frequency f_{θ} by generalized least squares (Draper & Smith 1998).
Searching for the best approximation, we fit this model at a series of test frequencies f_{θ} in the range between the minimal and maximal reasonable frequencies (as constrained by the width of the sliding windows and the full timespan of the observations). Figure 6 shows an example of these fits. The left panel shows loglikelihoods of the fitted models for pulsation period and peaktopeak amplitude as a function of frequency f_{θ} for the FU twinpeak CEP1748. The best fits corresponding to the highest peak of the loglikelihood profiles are shown in the center and the right panels. The obtained fits appear to capture the trend and, in most cases, the dominant quasiregular oscillatory variations, although the model is clearly only a rough approximation.
In order to quantitatively characterize the period and amplitude modulations in Cepheids, we extract the key parameters of the model (2): (a) the slope of the trends β_{P} and , see Eq. (2); (b) the amplitudes of the bestfit oscillatory components ΔP and ΔA_{1}, which are defined as and , respectively; and (c) the frequencies f_{P} and of the dominant oscillatory components. The estimated parameters for A_{1} and P are given (with estimated standard errors) in Appendix E. However, there are several important facts to keep in mind. As explained in Appendix C, the estimated frequencies can be “aliased”, that is, frequencies above the kernel method’s upper detection limit can be perceived as lower frequencies. Additionally, the amplitudes ΔP and ΔA_{1} of fast modulations are biased downwards. This bias depends on the modulation frequency, and thus can be corrected through an empirical estimated relationship between the two, if the modulation frequency is well known (see Appendix C). However, since it cannot be decided based on our data whether the estimated frequency is correct or not, the correction may be flawed, particularly for the period modulations in our control sample, for which Fig. D.2 suggests the possibility of relatively fast modulations.
We explore the distribution of the estimated fluctuation parameters among the different Cepheid groups and with respect to the average period in the next two subsections (with the above caveats kept in mind). For amplitudes of the modulations, we give both biascorrected and uncorrected versions; however, our conclusions do not differ in the two cases.
Fig. 6
CEP1748 shown as example of the heuristic model (2) of the changes of the pulsation parameters. Left panel: the loglikelihoods for the pulsation period (in black) and the peaktopeak amplitude (in blue; the loglikelihood functions were shifted to have common minimum at zero). The other two panels show the fits using the frequency yielding the highest likelihood (blue dots) against the sliding window estimates (thick black line), their 95% confidence interval (thin black lines), the stable value from the reference fit (thick orange line) and its confidence band (orange band). The middle panel presents the variations of the pulsation period, the right panel those of the peaktopeak amplitude. 
4.1.2 Trends and fluctuations of pulsation periods
Figure 7 shows the parameters of the trend + oscillation model fitted to the variations of the pulsation period in the four Cepheid groups. As a colorcoding of twin peak frequency separations shows, there does not appear to be a clear relation between modulation frequency and twin peak frequency separations, despite this being commonly assumed. Any tentative indications of such a relation are strongest for the T.FO group.
The approximate trend component (leftmost panel) is on average higher in the twinpeak groups than in the control groups, both for fundamentalmode and firstovertone Cepheids. Inspecting the two twinpeak overtone stars with the lowest trend values (CEP1536 and CEP1561), the time series of periods from the sliding window fits confirms the absence of an overall trend, and the reason for the appearance of twin peaks in their secondary periodogram may be due to the comparatively strong fluctuating component in their pulsation period, and to some very significant changes in A_{1} for CEP1536.
The relative amplitude ΔP∕P_{cat} of the oscillatory component in the third panel of Fig. 7 is higher in the overtone groups than in the fundamental mode groups. The shift between different modes is larger than the difference between the control and the twinpeak groups both within the fundamental and the overtone Cepheid sample. This supports that twin peaks are more likely to be the result of longterm, trendlike period instability rather than of fluctuations. The fact that twinpeak groups have, on average, lower characteristic frequencies of period fluctuations than control groups, regardless of pulsation mode, agrees with this: slower semiregular variations are more likely than fast oscillating ones to induce twin peaks in residual periodograms when prewhitening is carried out with constant light curve parameters.
4.1.3 Trends and fluctuations in brightness amplitude
Figure 8 shows the temporal variations of the first harmonic amplitude A_{1} in the four Cepheid groups. The trends in A_{1} (leftmost panel) show a pattern similar to the trends in the pulsation period: they tend to be on average higher in the twinpeak groups than in the control groups, though the difference is less prominent than with the pulsation period. The behavior of frequencies of the oscillatory term in the first harmonic amplitude shows no difference across the groups. The relative amplitudes of the oscillatory term in A_{1} are on average higher among first overtone Cepheids than among fundamental mode objects, especially when we consider the biascorrected amplitudes.
In summary, our results suggest that the twinpeak phenomenon is indicative of trends as well as slow, relatively highamplitude fluctuations of primarily the pulsation period. Particularly strong changes in the amplitude may also be identified by twin peaks, although this ability is limited to only the strongest or most trendlike cases in OGLE data. Soszynski et al. (2008) found 28% of the FO and 4% of the FU sample to show the twinpeak phenomenon. However, our study suggests also that potentially interesting cases where the modulations have oscillatory or stochastic character on significantly shorter timescales than the window length can be missed, as was found in the case of CEP0727 and CEP1638 (see Sect. 3.1). Using the twinpeaks phenomenon only to identify cases of modulated pulsation might therefore miss scientifically interesting subclasses, and can give biased estimates about the occurrence and typology of the modulations.
Fig. 7
Distribution of the logarithm of the absolute trend β_{P} (leftmost panel), the frequency of the oscillatory term f_{P} (second panel) and the logarithm of the uncorrected and corrected relative amplitude of the modulation, that is, the amplitudes ΔP and ΔP_{c} of the oscillatory term divided by the catalog period P_{cat} (third and rightmost panels) in model (2), for the changes in the pulsation period of all Cepheids in our sample. Thecolors in the twinpeak groups indicate the absolute value of the separation of the twin peaks, ranging from light blue for small separations to dark blue for large ones. The label T.FU and C.FU corresponds to fundamentalmode twinpeak and control Cepheids, respectively; similarly, the labels T.FO and C.FO indicate first overtone twinpeak and control cepheids, respectively. The T.FO group includes the 5 cepheids selected because of their strong amplitude changes. 
Fig. 8
Distribution of the logarithm of the absolute trend (leftmost panel), the frequency of the oscillatory term (second panel) and the logarithm of its relative size, that is, the (uncorrected and corrected) size of the changes ΔA_{1} and ΔA_{1c} divided by Ā_{1}, the first harmonic amplitude of the average stablemodel (third and rightmost panels, respectively), estimated from model (2). Groups and their labels are the same as in Fig. 7. 
4.2 Light curve variations versus physical parameters
In Fig. 9 we plot the (tenbase) logarithm of absolute values of period trends and biascorrected amplitude of period fluctuations, log ΔP_{c}, against average logarithmic pulsation period log P_{cat}. The left panel shows the wellknown monotonically increasing relationship between log P and log β_{P} with large scatter. This is in agreement with the trends observed for period changes by means of OC diagrams (e.g., Szabados 1983; Pietrukowicz 2001; Turner et al. 2006) and is a consequence of evolutionary timescales, which are shorter for the highermass longperiod Cepheids (e.g., Bono et al. 2000; Fadeyev 2013; Anderson et al. 2016b). The large scatter can be in part due to the 12year timespan of OGLE: this may be too short to ascertain whether a slow, trendlike change is indeed a portion of an evolutionary longterm period change, or only a fluctuation on timescales longer than the OGLE timespan.
Our selection procedure, namely random choice from all Cepheids based only on the presence or absence of the twin peak, implied that our control FO sample have on average shorter periods than stars in the twinpeak FO group, as is discernible from Fig. 9. A similar, though smaller and less clean separation is also apparent for FU Cepheids.
The right panel of Fig. 9 suggests that period fluctuations become stronger with increasing average period. Such a trend is not consistent with an interpretation in terms of the lighttime effect. We further find no candidates for binarity based on the lighttime effect within this sample: the values of ΔP ∝ a_{1} sin i determined are too low. Oscillatory period fluctuations have previously been reported for individual longperiod Cepheids (e.g., RS Puppis, see Berdnikov et al. 2009). However, to the best of our knowledge no such relationship has yet been firmly established. We are currently extending this study using the full sample of OGLEIII and IV Cepheids in the Magellanic System (Soszyński et al. 2015b), which will enable a more detailed investigation of this relation. Despite the different observational cadences over the LMC and SMC through OGLEIII and IV, this will provide a more comprehensive picture of the distribution of modulation parameters using a much larger sample in part with a longer timespan, and is aimed at enabling a more indepth physical interpretation of these phenomena.
As a consequence of the periodluminosity relation of Cepheids, similar correlations exist between these parameters and mean I or V magnitude. However, no other relation with color or position on the colormagnitude diagram (CMD) was found. Figure 10 shows that Cepheids found to be variable in period or amplitude occupy all regions of the log P − log A parameter space covered by our sample. This further corroborates the ubiquity of light curve modulation among Cepheids. Smolec (2017) investigated all fundamentalmode Cepheids in the OGLE Magellanic Cloud collection (Soszynski et al. 2008, 2010; Soszyński et al. 2015b). However, in that study, detection was based on the presence of a particular symmetric doublet structure in the residual periodogram after a standard constantmodel prewhitening, which is aimed primarily to detect strong, strictly periodic variations, and does not capture all the possible manifestations of a multiscale or not necessarily strictly repetitive modulation. However, it is capable to identify highfrequency periodic modulations, which would be missed by the threeyear sliding window. Notably, none of the 53 Cepheids discussed in the present work is mentioned in Smolec (2017). Our ongoing analysis of the full OGLE Cepheid collection using the kernel method will enable further insights and a more complete comparison. Within its range of detectable modulation frequencies, the sliding window method is much more sensitive to pick out a variety of complex modulation patterns than classical prewhitening analyses that implicitly assume periodicity of modulations.
Fig. 9
Period trend ∣ β_{P} ∣ (left panel) and size of period fluctuations ΔP (right panel) as estimated from model (2) for the 53 Cepheids. Downward pointing triangles: control overtone Cepheids, upward pointing triangles: overtone twinpeak Cepheids, diamonds: fundamentalmode control Cepheids, squares: fundamentalmode twinpeak stars. In the left plot, negative trends are indicated in blue, positive ones in orange. 
Fig. 10
Doublelogarithmic periodamplitude diagram for program Cepheids using the original OGLEIII periods and peaktopeak amplitudes. Constant objects are shown as open circles, variable objects as filled circles. Left panel: changing periods. The color corresponds to column T_{P} of Tables 1 and 2, from cyan (0%) to magenta (100%). Right panel: variable amplitudes. The color corresponds to column of Tables 1 and 2, from cyan (0%) to magenta (100%). 
5 Conclusions
In this paper, we apply local kernel modeling, a wellknown method of statistics to investigate periodic and nonperiodic temporal variations of Cepheid light curve parameters. We apply this method to 53 classical Cepheids from the OGLEIII catalog of variable stars (Soszynski et al. 2008), selected according to the presence or absence of a secondary oscillation frequency close to the primary (referred to as “twin peak”), or because of very strong, visually obvious amplitude changes. We compare on the one hand the behavior of fundamentalmode with firstovertone Cepheids, and on the other, the behavior of Cepheids that exhibit twin peaks in secondary periodograms with those that do not. Our method yields estimates of the pulsation period and light curve parameters as smooth functions of time; we estimate the significance of the found deviations from the bestfit constant parameters by bootstrapping residuals (Monte Carlo methods) and by multiple hypothesis testing procedures with respect to the bestfit stable model.
We find period modulations to be probably a very frequent, possibly ubiquitous phenomenon among Cepheids. Changes in light curve amplitudes also seem to occur frequently, although they are harder to identify reliably. Our results suggest that twin peaks are related to instabilities on longer timescales in the period or in the light curve parameters. The characteristic size of these instabilities for Cepheids are such that with a given time sampling pattern and photometricprecision, period changes are easier to detect than variations in the harmonic content of the light curve.
Over the sample of stars considered, we find a wide range of degrees to which amplitudes and periods can change, and timescales of the variations ranging from the shortest to the longest detectable. This suggests an extension of the range of both amplitude and period changes occurring in Cepheids to beyond our detection limits, which are imposed by the time sampling and precision of the photometry used. Applying our method to more densely sampled and more precise photometry from space (K2, CoRoT) should allow the detection of smaller amplitude changes on complementary timescales.
We detect several different types of behavior (nearlinear, oscillatory, and stochastic) among period and amplitude variations. Specifically, we find that the majority of Cepheids exhibit period changes beyond or different from ones expected from secular evolution (see also Poleski 2008). These may rather be related to other timedependent phenomena such as convection, granulation, rotation (spots), (episodic?) massloss, secondary modes at frequencies close to the primary, or other forms of modulation such as the Blažko effect seen in RR Lyrae stars.
For Cepheids exhibiting “twin peaks”, we find that prewhitening with our timedependent bestfit light curve estimates generally removes the “twin peak” from the secondary periodogram. After a timedependent prewhitening, several overtone twinpeak stars exhibited weak secondary peaks, in majority at longer periods than the primary. No such secondary periods were detected in fundamentalmode stars.
As a next step, we will apply this method to a much larger sample of stars, such as the OGLEIV catalog of Cepheid variable stars, and datasets with higher photometric precision. In doing so, we will explore the overall characteristics and occurrence rates of period and amplitude variations as well as their dependence on the stellar properties. Revealing the timedependence of Cepheid light curves will thus yield new constraints on the structure of the outer envelopes of Cepheids and serve to understand the interaction between the pulsations and the medium they propagate inside of. Upcoming spacemissions such as TESS and PLATO willprovide highquality photometry of many Cepheids, allowing us to explore a larger range of amplitude variations and improve population statistics.
Acknowledgements
We thank the anonymous referee for valuable comments that have helped to improve the quality of this manuscript. R.I.A. acknowledges financial support from the Swiss National Science Foundation. The research made use of the public OGLEIII database (http://ogledb.astrouw.edu.pl/~ogle/CVS/) and the NASA’s ADS bibliographic services.
Appendix A Detailed statistical methodology
A.1 Local estimation
Windows. Corresponding to the aim of the study, we performed nonlinear harmonic series fitting, optimizing also over period, in threeyear windows of the time series. We fixed the centers of the windows at a time grid with a 30 day separation, which allowed us to follow the period, amplitude and shape variations with a reasonable temporal resolution, thus fitting about 130 windows for each of our targets.
Weights. In order to obtain a smooth picture of the local changes in the light curve shape with most emphasis on the observations close to the center of the window, we used a particular weighting scheme, which combined kernel weighting with the usual inverse squared error weights (Fan & Gijbels 1996). Suppose that we have N observations Y _{i} with errors σ_{i} at times t_{1}, …, t_{N}. For a window centered at time τ_{k}, we defined the first component of the weight of an observation at t_{i} as (A.1)
with a Gaussian kernel , and with bandwidth h = 182.5. This component ensured that the observations close to the window center contribute more to the fit than observations farther off, and cut the effect of observations outside a threeyear interval. As an example, Fig. A.1 shows a window centered at JD = 2850 as a grayshaded area of a part of the time series of OGLELMCCEP1621 (top panel). The Gaussian kernel, computed at the times of the observations, is shown in the middle panel in black. We used the inverse squared errors; (A.2)
as the other component of the weights, presented in the middle panel of Fig. A.1 as red spikes. The final weights, shown in the bottom panel, were determined by (A.3)
This scheme ensured both that we obtain an estimate based on the most relevant observations, and that data with comparatively large errors have a weaker influence on the estimate than data with smaller errors.
Model formula. Within each window, and using the above described weighting procedure, we fitted a harmonic + thirdorder polynomial model of the form (A.4)
where are assumedto be independent Gaussian errors.
Polynomial order. Since a neglected nonlinear trend can cause bias in the frequency estimate, we included a thirdorder polynomial trend into the local model. Inspection of the fits for the Cepheids suggests that the mean magnitude can occasionally vary rapidly, and in such periods, the effect influences visibly the frequency estimate, as can be seen in the top panel of Fig. A.2. Mean magnitude variations within the window were thus accounted for by the polynomial term in the local model.
Fig. A.1
Weighting scheme. The top panel shows a portion of the observed time series of OGLELMCCEP1621, with the fitted window highlighted by a gray background. The middle panel shows the two weight sequences: in black, the kernel weights, in red, the inverse error weights. The bottom panel exhibits the combined weights used in the fit. 
Fig. A.2
Difference of fits with locally constant (red), locally linear (blue) and locally cubic (black) polynomial terms on the example of OGLELMCCEP1833. The top panel shows the estimates of frequencies as a function of time. The middle panel shows the temporal dependence of the estimated mean magnitude. The bottom panel presents the coefficient of the cubic term in the locally cubic model versus time. 
Harmonic order M. For a sliding window fit, the choice of M is crucial. It must be fixed and the same for all windows, otherwise the appearance of new significant terms or the dropout of formerly significant ones may result in unreasonable jumps of model parameter estimates. As the required M can be very high for some stars such as bump Cepheids and low for others with sinelike light curve, we have to choose M individually.The choice should be such that it allows for some coefficients becoming temporarily significant and falling out again. It is also desirable to allow for a somewhat higher order than that of a stable model for the complete observation span, since if there are in truth temporal variations, a stable model sees a sort of average, and we can get a too simple model. However, a too high M adds many nonsignificant parameters to the model, which in turn can change the values estimated for the significant terms. Thus, the individual values were determined by fitting a stable model with M_{init} = 15, performing a statistical model selection procedure based on the Bayes Information Criterion (Schwarz 1978). After finding a maximal order M_{final}, we finally fixed the harmonic order at M = M_{final} + 2.
Model fitting. According to our aim of investigating the temporal behavior of the frequency and the light curve parameters, we treated the problem as a nonlinear model, optimizing over the linear parameters a_{0}, …, a_{3}, s_{1}, c_{1}, …, s_{M}, c_{M} as well as over the frequency f. Optimization was initialized with the stable model parameters, and was restricted to an adapted frequency interval around the catalog frequency.
A.2 Error analysis
We conducted a careful error analysis to assess the significance of the presumably small temporal variations in variability parameters. We therefore employed several different methods of determining confidence intervals for the estimated timevarying parameters, and to compare the timevarying model to the temporally stable alternative reference model for each Cepheid.
Uncertainty of the estimates. For local linear regression models with normally distributed errors, uncertainty of the estimated parameters can be obtained by closed formulae (see Eq. (3.6) in Fan & Gijbels 1996), if the errors are known. In that case, the error on the estimates follows a multivariate normal distribution. However, we treated the pulsation period of the Cepheid as a parameter to be fitted. Hence, our model became a nonlinear regression model, so this formula and the multivariate normality of the estimates are only approximate (and asymptotic). We therefore used a double procedure to obtain uncertainty of the parameters:
 1.
The above cited formula from Fan & Gijbels (1996), extended for the nonlinear model.
 2.
Parametric bootstrap, using the given observational error bars as the standard deviation of the generating normal distribution at each epoch, assuming independent errors, and using the fitted nonstable models interpolated to each observing time as the errorfree light curve. The sliding window estimation was then performed on all repetitions. Taking the quantiles 0.025 and 0.975 of the obtained estimates yielded the 95 % bootstrap confidence intervals (CIs), plotted around the estimates in Figs. 3, 6 and D.1−D.4.
In all windows of all Cepheids, these two estimates were very close together, so in our nonlinear model, nonlinearity does not cause large deviations in the statistical distribution of the parameter estimates from approximations based on theory.
Stable reference model. We would like to know whether a stable model, with parameters constant over time, can yield a plausibleexplanation for the estimates from the local sliding window model. For this, we simulated the bestfit stable model with added noise and repeated the sliding window estimation 250 times for each Cepheid.
However, in the case of the stable model, the added noise cannot be simulated from the error bars using a Gaussian model, since the residuals from the stable model are strongly overdispersed with respect to the given error bars. It was necessary to allow for an underestimation of the errors, if we want to have a plausible zero hypothesis. Therefore, we applied a nonparametric bootstrap of the standardized residuals. We scaled all the residuals r_{1}, r_{2}, … at times t_{1}, t_{2}, … with the error bars σ_{1}, σ_{2}, … at that time (the mean of the residuals is zero by construction, so we don’t need to subtract it). If the error bars are underestimated by a common factor with a good approximation, and the error distribution is a locationscale model (not necessarily Gaussian, but admitting also heavytailed distributions such as the Lorentzian/Cauchy), this procedure obtains a homoschedastic error sample s_{1} = r_{1}∕σ_{1}, s_{2} = r_{2}∕σ_{2}, … We resampled these with repetition, obtaining a (scaled) sample (each of thestarred residuals is equal to an arbitrary one of the original scaled residuals). At each time t_{1}, t_{2}, …, we computed a simulated noise value by scaling back the scaled residual with the local standard error on the observation: This simulated noise value was added to the fitted magnitudes of the stable reference model, obtaining simulated “observed” magnitude values of a noisy Cepheid with stable parameters and errors consistent with the assumption of stability.
For the detection of such possibly very fine effects as shortterm period and amplitude changes in Cepheids, we considered also the effect of the precision of times, although this may be a minor effect. The OGLEIII database gives the HJD dates in days with a precision of five digits. In order to assess the effect of rounding, we replaced the times t_{1}, t_{2}, … by values jittered in their sixth digit.
Finally, the time series was fitted by the local kernel model, obtaining an estimate for every parameter value at every window center τ_{i} (, denoted generally by θ^{*}(τ_{i})). This procedure was repeated R = 250 times, and at each window center, the quantiles 0.025 and 0.975 of the obtained estimates were taken to get the 95% pointwise bootstrap confidence intervals on our hypothetic stable Cepheid. We performed the bootstrap procedure also with unchanged times. We found that the jittering has a negligible effect on the resulting confidence bands, and gives visibly larger CIs only rarely in short time intervals. Thus, our results are robust against the finite precision of the times. We indicate the CIs obtained with jittered times as orange band around the estimated stable bands in Figs. 3, 6, and D.1−D.4.
The procedure is approximate for many reasons (the fitted Cepheid is itself an estimation, the residuals follow a correlated joint distribution with variances somewhat different from , the procedure is based on the assumption of homogeneous underestimation of errors, and although the observational error distribution is not restricted to Gaussian, it must nevertheless be a locationscale distribution), but it accounts for the two dominant effects, namely, the overdispersion of the residuals and the implications of the irregular sparse time sampling on the sliding window estimates.
A.3 Attribution of significance: multiple hypothesis testing
The null hypothesis to be tested for each model parameter (pulsation period, amplitudes) at each window center τ_{i} is that the true local value of θ is in fact equal to the stable value . We computed a probability of seeing a discrepancy between the local and the stable value equal to or higher than the one determined under this null hypothesis, then compare this probability to a predefined confidence level α (0.05 in our study): if the probability of the found discrepancy is higher than this level, we cannot reject the null hypothesis of the local estimate being equal to the stable one.
To compute the probability, we supposed that the discrepancy follows a Gaussian distribution with mean 0 and variance computed from the empirical distribution of the repetitions obtained with the bootstrap procedure described above. This distribution represents the null distribution, that is, when thepulsation parameters of the Cepheid are in fact stable. We visually confirmed that the repetitions do follow a Gaussian distribution using quantilequantile plots (for a short description, see e.g., Süveges 2014). Using this null distribution, we could then compute a pvalue, namely, the probability that such a discrepancy is given by the kernel method at window center τ_{i} when the Cepheid is a steadily pulsating one (with the bestfit stable reference parameters).
These pvalues are pointwise and have been computed for each window center (about 130 over the total time span of OGLEII and III). Their number itself raises a problem. Generally, if we have a true zero hypothesis, and we assess significance of an alternative at level α separately K times using completely independent data, even in the case of a true zero hypothesis we can expect approximately Kα significant values (that is, false positives), simply by randomness. Thus, finding a few significant pointwise pvalues should not necessarily be considered as having found a significant global discrepancy from stability. Multiple hypothesis testing methods aredeveloped to deal with this situation.
Moreover, pvalues computed for neighboring windows are strongly correlated, due to the overlap in the data used for each local estimate. If by sheer randomness, at τ_{i} there was a pvalue smaller than our confidence limit α, there is an increased probability that we will have a similarly small pvalue at the next or next few window centers, creating a false effect of a longer period where the pulsation parameters were significantly different from the stable one. Thus, particular versions of multiple hypothesis testing procedure must be used, which in addition takes into account the dependence between neighboring pvalues. Such procedures have been developed (among others) by Hommel (1988); Hochberg (1988); Benjamini & Hochberg (1995); Benjamini & Yekutieli (2001) and Meskaldji et al. (2013). Here, we applied the procedure by Benjamini & Yekutieli (2001) that imposes only few restrictions on the dependence structure and is one of the most powerful among the alternatives. The time intervals during which this method indicates significance are highlighted in Figs. 3, 6, D.1−D.4 in gray.
Appendix B Simulated instability
We investigate how the detectability of the minute effects reported in this work are affected by data structure, and quantify the detection limits of the kernel smoothing method to trace modulations of the pulsations using relatively sparse and unevenly sampled data. To this end, we have simulated trendlike, as well as combined trend and oscillationlike variations of pulsation periods and amplitudes, for the range of parameters representative of our sample.
Noisefree light curves were generated as follows. First, we determined the instantaneous phase by (numerical) integration of the timevarying frequency: (B.1)
where we replaced ϕ_{0} = 0 without loss of generality. We substituted f(s) by combinations of a linear trend and periodic fluctuations in pulsation period: (B.2)
Varying first harmonic amplitudes were produced by
where the notation agrees with that of Eq. (2), and and were taken from similar linear+oscillatory model fits to the timevarying harmonic coefficients, such that .
Next, ϕ(t), s_{1} (t) and c_{1} (t) were inserted into the model formula (A.4) to obtain the pure noisefree light curves at times t_{i}: (B.5)
where a_{0} and the m ≥ 2 harmonic coefficients were kept constant. We did not add longterm polynomial components in these simulations.
To create realistic simulations, we adopted real time samplings and our estimates of fitted trend+oscillation model parameters for three real Cepheids: CEP1405 (FO twin peak featuring oscillations in A_{1}, as well as a trend and oscillation in P), CEP1748 (FU twin peak with significant variations only in A_{1}) and CEP2470 (FU twin peak with significant variations only in P). We tested the performance of the kernel method in several different modulation types from simple to complex: (1) using period trends of increasing slope up to the trend found in OGLELMCCEP1405, (2) superposing oscillations of varying frequency and amplitude to the trend of OGLELMCCEP2470 and (3) simulate the full estimated model of OGLELMCCEP1748 and OGLELMCCEP1405. To these noisefree light curves, we added independent Gaussian error terms using the OGLE magnitude error estimates as the standard deviation of the Gaussian, generating 250 independent repetitions of noisy light curves for each investigated pure light curve. Finally, we repeated the sliding window estimation on each of the repetitions, using the same tuning parameters as in the investigation of the observed Cepheids.
B.1 Trends in the pulsation period
To investigate the method’s power of detecting trends in pulsation period, we simulated a stable light curve with parameters identical to those of CEP1405, adding linear trends in pulsation period as β_{P} = aβ_{P,1405}, where β_{P,1405} = −1.8 × 10^{−7}d∕d (see Table E.1), and a ∈{0, 0.1, 0.25, 0.5, 1, 2}. Figure B.1 shows the median estimates and the 0.025 and 0.975quantiles for a = 0, 0.25 and 1. The local kernel estimation is unbiased apart from end effects at the beginning and end of the observation span, both for the period trend and the constant amplitude. The characteristic OGLE time sampling does not leave a strong systematic imprint on the estimates. The simulation with a = 0 (i.e., a perfectly stable light curve) reveals no tendency to incorrectly suggest a nonexisting trend or fluctuation.
Fig. B.1
Simulated trends in pulsation period (top panel) combined with a constant amplitude (bottom panel) as estimated by the localkernel method. Thin red, green and blue solid lines in the top panel: the simulated values using slopes 0, 0.25 and 1 times that of OGLELMCCEP1405 (such that the blue line coincides with the trend component found for this Cepheid). The true amplitude, which had the same constant value in all simulations, is represented by a blue line in the bottom panel. Dotted red, green and blue lines: the pointwise median of the kernel estimates on the 250 repetitions. The colors show a 95% symmetric confidence band around the median. 
Inserting nonzero trends into simulations of CEP1405, we find that our method would have detected trends as small as a quarter of the observed trend, which can be seen as a green line in Fig. B.1. Both the trend in the period and the stability of the amplitude is well estimated, which suggests that the model is able to correctly disentangle effects in pulsation period and amplitude.
B.2 Combined oscillation and trend in the pulsation period
Trends are smooth long timescale variations and are usually wellestimated by local models. However, modulations on shorter timescales can also occur in the Cepheids, and the length of the temporal (sliding) window crucially affects what timescales can be detected.
To simulate the case of short timescale modulations, we adopted the parameters for CEP2470, whose period and amplitude fluctuations are near the detectability limit of our threeyearwide sliding window. We superposed five different oscillations of the pulsation period to a period trend of β_{P} = β_{P,2470} = −5.983 × 10^{−8}d∕d. Four of these had the same oscillation amplitude Δ_{P} = Δ_{P,2470} as the bestfit estimate of CEP2470, with frequencies f_{P} = a f_{P,2470} with a ∈{0.25, 0.5, 0.75, 1}. The fifth simulation used Δ_{P} = 4Δ_{P,2470} and f_{P} = f_{P,2470}. For all simulations, we kept all other parameters at the values of the stable reference model of CEP2470 (no variations in the amplitude were added).
Fig. B.2
Simulated oscillations of the period, superposed to the period trend estimated for CEP2470. The parameters of the simulated modulations are shown in the upper right corner of the panels. The true simulated modulations of the period are shown with solid red lines. The 0.025, 0.975 quantiles and the median of the kernel estimates are indicated by the colored band and the dotted thick line, similarly to Fig. B.1. The thin black horizontal line in the top panel shows the full width of the sliding window, to ease comparison with the period of the fluctuations. The rugplot added to the xaxis indicates the observation times. 
Figure B.2 illustrates the results of this investigation. As expected, slow period fluctuations that are on the order of or longer than the sliding window duration are wellestimated by the local kernel method (two top panels). The size of faster period fluctuations (f_{P} ≥ 0.75 f_{P,2470}) is increasingly underestimated (middle panel), and the kernel procedure eventually fails to detect fluctuations on timescales that approach the time interval within which our weights are nonnegligible. The kernel procedure does, however, allow us to recover the correct fluctuation frequency if the intensity of such fast fluctuations is increased (Δ_{P} = 4Δ_{P,2470}, bottom panel), even though the fluctuation intensity is then underestimated by about a factor of four. Thus, longtimescale fluctuations (relative to sliding window size) are well estimated, whereas the detectability decreases for shorter fluctuation timescales and depends on the intensity of the fluctuations. The characteristic size of the fluctuations is underestimated. This bias is discussed in detail in Appendix C. The harmonic amplitudes A_{1}, which were kept constant in these simulations, were correctly reproduced by all five simulations and are not shown here for brevity.
Fig. B.3
Local kernel estimates of the modulation parameters of a Cepheid simulated using the estimated parameters of CEP1748. The notation and the symbols are the same as in Fig. B.2. 
B.3 Adding amplitude variations
We have verified the ability of the local kernel method to correctly detect amplitude modulations in the case of stable pulsation periods. To this end, we simulated variations in A_{1} using CEP1748’s estimated A_{1} fluctuations, setting both β_{P,1748} = 0 and Δ_{P,1748} = 0. As Fig. B.3 shows, amplitude modulations on timescales longer than the window length are wellestimated (bottom panel). However, based on the bias in the amplitude estimates of the period fluctuations, we expect similar detection bias regarding fluctuations of A_{1}; we will discuss these biases in detail in Appendix C.
Combining amplitude and period fluctuations, we have simulated a case in analogy to CEP1405. Here we adopted a fluctuation timescale that would be readily detectable using a 3year sliding window (cf. Appendix B.1). As Fig. B.4 shows, the kernel method provides unbiased estimates of both the variations in A_{1} and P, and reveals no issues related to the separation of the two phenomena.
Appendix C Reliability and bias of the modulation parameter estimates
Appendices B.2 and B.3 used simulated amplitude and period fluctuations to show that both types of fluctuations (separate or combined) are accurately recovered if the timescale of the fluctuations is at least halfof the sliding window’s timespan. However, they also showed that the sliding window introduces bias in the estimated intensity of period and/or amplitude fluctuations, which implies that estimates obtained via a heuristic model (2) are also potentially biased. Here, we investigate how the limitations of the kernel method influence our results, specifically regarding the frequency and amplitude of any detected fluctuations in A_{1} and P.
Fig. B.4
Local kernel estimates of the modulation parameters of a Cepheid simulated using the estimated parameters of CEP1405. The notation and the symbols are the same as in Fig. B.2. 
We selected five Cepheids spanning a range of different catalog pulsation periods (CEP1527, P_{cat} = 1.49 d, CEP2217, P_{cat} = 2.31 d, CEP2191, P_{cat} = 4.21d, CEP1140, P_{cat} = 8.19 d, and CEP1833, P_{cat} = 19.16 d) and different time samplings. We then simulated ten different periodic period modulations for each Cepheid that cover the frequency range detectable by our sliding window using modulation frequencies f_{i} = i∕3652.5, i = 1, 2, …, 10. The basic model in each case was the catalog pulsation period and the harmonic coefficients determined using the stable reference model. Here we assumed high modulation amplitudes (0.002 c∕d) in order to clearly illustrate the relationship of the detection bias as a function of fluctuation frequency (the expected ratio between true and estimated modulation amplitude does not depend on the signaltonoise ratio). We estimated each of the simulated light curves (with 200 different added Gaussian white noise sequence for each) with our sliding window method, and fitted model(2) to each of the estimated time series of periods.
Figure C.1 presents the results of this bias estimation. Its top panel shows the distribution of fitted modulation frequencies as a function of the true modulation frequency. For modulation frequencies ≲ 0.0015 c∕d, that is, a twoyear modulation period, our method achieves essentially unbiased frequency estimates: in most cases the estimates for all stars are very close to the true modulation frequency. Specifically, we do not find a dependence of the fluctuation estimates on the catalog pulsation period. The simulations involving the longestperiod Cepheid (CEP1833) however exhibit significantly higher variance, in particular at the lowest frequencies. This may be due to its long period: fewer full light curve cycles are observed in one window, and thus the change in period cannot be estimated as precisely as for a shorterperiod Cepheid. Above the limit modulation frequency of 0.0015 c∕d, the estimated frequencies for all five simulated Cepheids are much more scattered, and interestingly, their median often falls not on the true value, but on its yearly alias.
Fig. C.1
Modulation frequency (top panel) and the ratio between the modulation amplitude and the true amplitude (bottom panel) estimated by model (2) applied to local sliding window estimates of simulated period modulations. Colored dots show the median of the estimates based on 200 different Gaussian noise realizations added to each of five simulated modulated light curves, and vertical segments indicate the range between the 0.025 and 0.975 quantiles. Cyan: using parameters of CEP1527, green: CEP2217, blue: CEP2191, dark red: CEP1140, orange: CEP1833. All simulations used the same sequence of true modulation frequencies. For visibility, the dots based on different Cepheids are horizontally offset from the true value. In the top panel, the solid black line indicates x = y, i.e., the location of accurate estimates, with short horizontal segments aiding the eye as to where horizontally offset points should lie. The dashed line indicates the yearly alias of the x = y line. In the bottom panel, the horizontal gray line indicates y = 1, the location of the ratios in case of unbiased estimates. The solid black line is the best linear fit using the five medians at each modulation frequency corresponding to the seven lowest true modulation frequencies. 
Repeating the above procedure for modulations in the amplitude instead of period with a chosen modulation amplitude of 2 mmag, we obtain a very similar bias pattern, shown in Fig. C.2. To assess the bias of the modulation amplitude estimates from model (2) for both period and amplitude modulations, we computed the ratios between the estimated and the true modulation amplitude, as shown the bottom panels of Figs. C.1 and C.2. The plots suggest that this ratio is nearly linear as a function of the true modulation frequency in the frequency range that can be reliably detected. Fitting this relationship (between 0 < f_{P} < 0.0018 for period modulations and between for amplitude modulations), we obtain the empirical biascorrection formulae
which we use to estimate the true underlying modulation in Table E.1 and in Figs. 7−9.
Fig. C.2
The modulation frequency (top panel) and the ratio between the modulation amplitude and the true amplitude (bottom panel) estimated by model (2) applied to local sliding window estimates of simulated amplitude modulations. The notation is the same as in Fig. C.1 The solid black line is the best linear fit using the five mediansat each modulation frequency corresponding to the five lowest true modulation frequencies. 
Furthermore, the above simulations can be used to form an idea how much the particular time sampling influences locally the shape of the estimated pattern. Figure C.3 presents the estimated modulations of the period in a case where the frequency of the modulation was fixed around the upper detection limit of our threeyear window (corresponding to a period of two years). The simulation in the uppermost panel uses the time sampling, catalog period and stable harmonic parameters of CEP2217 (FO, harmonic order H = 3, P_{cat} = 2.31 d), the middle one, that of CEP1140 (FU, harmonic order H = 12, P_{cat} = 8.19 d), and the bottom panel, that of CEP1833 (FU, harmonic order H = 12, P_{cat} = 19.16 d). The systematic distortions caused by the different time samplings are noticeable, as well as the increasing statistical uncertainty due partly to the need of estimating more parameters from a similar amount of data. The shape of the estimated signal for the characteristics of CEP2217 and 1140 is remarkably stable, the signal, even at this relatively high frequency, is doubtless present, albeit underestimated, and the deformations depend predominantly on the time sampling, but only little or very little on the noise. In our other two simulations, using CEP1527 and 2191, the estimated signal pattern is even less variable than using CEP2217. The larger variance of the estimates in the case of CEP1833 implies the higher uncertainty of the estimates of the frequency and the amplitude.
Fig. C.3
Local kernel estimates of the period modulation in simulations of a modulation of f_{P} = 0.00137 c∕d, using the stable parameters of CEP2217 (top panel), CEP1140 (middle panel), and CEP1833 (bottom panel). Solid dark red line: the true signal, dotted dark red line: median of the local kernel estimates on 200 random realizations, orange band: the stripe between the pointwise 0.025 and 0.975 quantiles of the estimates. 
Appendix D Figures of the temporal behavior of a few variability parameters
D.1 Pulsation period
Fig. D.1
Temporal variation of the pulsation period of the twinpeak Cepheids (solid thick black line), together with its bootstrapped pointwise CI (thin black lines). Heavy orange line: catalog period, orange band: 95% nonparametric bootstrap CI around this period based on the bootstrap procedureof Sect. A.2. The dotted black vertical lines indicate years after HJD − 2 450 000. The gray background highlights time intervals where the deviation from the stable reference model was found significant (cf. Sect. A.3). The times of the observations are indicated with a rugplot at the bottom of the panels. 
Fig. D.2
Temporal variation of the pulsation period of the control Cepheids. For the notation, see the caption of Fig. D.1. 
D.2 Amplitude of first harmonic term A_{1}
Fig. D.3
Temporal variation of the leading harmonic amplitude A_{1} of the twinpeak Cepheids. For the notation, see the caption of Fig. D.2. 
Fig. D.4
Temporal variation of the leading harmonic amplitude A_{1} of the control Cepheids. For the notation, see the caption of Fig. D.2. 
D.3 Light curve shape
Fig. D.5
Kernelestimated light curve shapes in two little or nonoverlapping windows for twinpeak fundamental Cepheids, inwhich the peaktopeak amplitude or the amplitude of the leading harmonic term was very different. Fitted curves are indicated by the lines of different colors, and the observations in the different windows, by the different colors and shapes of the symbols. The times of the window centers are given in the legend. The observations close to the window center and therefore more influential in the fit have darker red or blue colors, whereas those in the wings are shown in lighter shades. Approximate bootstrap confidence intervals for the curves are given as light red and light blue stripes around the estimated lines (Sect. 2.2 and Appendix A). 
Fig. D.6
Light curve shapes in two different windows for twinpeak overtone Cepheids. For the notation, see the caption of Fig. D.5. 
Fig. D.7
Light curve shapes in two different windows for overtone Cepheids showing strong amplitude changes. For the notation, see the caption of Fig. D.5. 
Fig. D.8
Light curve shapes in two different windows for control fundamentalmode Cepheids. For the notation, see the caption of Fig. D.5. 
Fig. D.9
Light curve shapes in two different windows for control overtone Cepheids. For the notation, see the caption of Fig. D.5. 
Appendix E Table of trends and fluctuations in the pulsation period and the first harmonic amplitude
References
 Alard, C. 2000, A&AS, 144, 363 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alard, C., & Lupton, R. H. 1998, ApJ, 503, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Anderson, R. I. 2014, A&A, 566, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Anderson, R. I. 2016, MNRAS, 463, 1707 [NASA ADS] [CrossRef] [Google Scholar]
 Anderson, R. I., Sahlmann, J., Holl, B., et al. 2015, ApJ, 804, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Anderson, R. I., Mérand, A., Kervella, P., et al. 2016a, MNRAS, 455, 4231 [NASA ADS] [CrossRef] [Google Scholar]
 Anderson, R. I., Saio, H., Ekström, S., Georgy, C., & Meynet, G. 2016b, A&A, 591, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Arellano Ferro A. 1983, ApJ, 274, 755 [NASA ADS] [CrossRef] [Google Scholar]
 Benjamini, Y., & Hochberg, Y. 1995, J. Roy. Stat. Soc., Ser. B, 57, 289 [Google Scholar]
 Benjamini, Y., & Yekutieli, D. 2001, Ann. Stat., 29, 1165 [CrossRef] [Google Scholar]
 Benkő, J. M., Szabó, R., & Paparó, M. 2011, MNRAS, 417, 974 [NASA ADS] [CrossRef] [Google Scholar]
 Benkő, J. M., Plachy, E., Szabó, R., Molnár, L., & Kolláth, Z. 2014, ApJS, 213, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Berdnikov, L. N., & Ignatova, V. V. 2000, in IAU Colloq. 176, The Impact of LargeScale Surveys on Pulsating Sta r Research, eds. L. Szabados, & D. Kurtz, ASP Conf. Ser., 203, 244 [NASA ADS] [Google Scholar]
 Berdnikov, L. N., Ignatova, V. V., Caldwell, J. A. R., & Koen, C. 2000, New Astron., 4, 625 [NASA ADS] [CrossRef] [Google Scholar]
 Berdnikov, L. N., Henden, A. A., Turner, D. G., & Pastukhova, E. N. 2009, Astron. Lett., 35, 406 [NASA ADS] [CrossRef] [Google Scholar]
 Berdnikov, L. N., Kniazev, A. Y., Sefako, R., Kravtsov, V. V., & Zhujko, S. V. 2014, Astron. Lett., 40, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Blažko S. 1907, Astron. Nachr., 175, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Bono, G., Caputo, F., Cassisi, S., et al. 2000, ApJ, 543, 955 [NASA ADS] [CrossRef] [Google Scholar]
 Burki, G., & Mayor, M. 1980, A&A, 91, 115 [NASA ADS] [Google Scholar]
 Burki, G., Mayor, M., & Benz, W. 1982, A&A, 109, 258 [NASA ADS] [Google Scholar]
 Derekas, A., Szabó, G. M., Berdnikov, L., et al. 2012, MNRAS, 425, 1312 [NASA ADS] [CrossRef] [Google Scholar]
 Derekas, A., Plachy, E., Molnár, L., et al. 2017, MNRAS, 464, 1553 [NASA ADS] [CrossRef] [Google Scholar]
 Draper, N. R., & Smith, H. 1998, Applied Regression Analysis (Wiley), 2nd ed. [Google Scholar]
 Dziembowski, W. A. 2016, Commmunications of the Konkoly Observatory Hungary, 105, 23 [NASA ADS] [Google Scholar]
 Eddington, A. S. 1919, MNRAS, 79, 177 [NASA ADS] [Google Scholar]
 Evans, N. R., Szabó, R., Derekas, A., et al. 2015, MNRAS, 446, 4008 [NASA ADS] [CrossRef] [Google Scholar]
 Fadeyev, Y. A. 2013, Astron. Lett., 39, 746 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, J., & Gijbels, I. 1996, Local Polynomial Modelling and Its Applications (London: Chapman and Hall) [Google Scholar]
 GarcíaVarela, A., Muñoz, J. R., Sabogal, B. E., Vargas Domínguez, S., & Martínez, J. 2016, ApJ, 824, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Goodricke, J. 1786, Roy. Soc. Lond. Philos. Trans. Ser. I, 76, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Hochberg, Y. 1988, Biometrika, 75, 800 [CrossRef] [MathSciNet] [Google Scholar]
 Hommel, G. 1988, Biometrika, 75, 383 [CrossRef] [Google Scholar]
 Kanev, E., Savanov, I., & Sachkov, M. 2015, in Eur. Phys. J. Web Conf., 101, 06036 [CrossRef] [Google Scholar]
 Kervella, P., Trahin, B., Bond, H. E., et al. 2017, A&A, 600, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kolenberg, K., Szabó, R., Kurtz, D. W., et al. 2010, ApJ, 713, L198 [NASA ADS] [CrossRef] [Google Scholar]
 Kolláth, Z., Buchler, J. R., Szabó, R., & Csubry, Z. 2002, A&A, 385, 932 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kovacs, K., Buchler, J. R., & Davis, C. G. 1987, ApJ, 319, 247 [NASA ADS] [CrossRef] [Google Scholar]
 Meskaldji, D. E., Thiran, J.P., & Morgenthaler, S. 2013, ArXiv eprints [arXiv:1112.4519] [Google Scholar]
 Molnár, L., Szabados, L., Dukes, et al. 2013, Astron. Nachr., 334, 980 [NASA ADS] [CrossRef] [Google Scholar]
 Moskalik, P., & Kołaczkowski, Z. 2009, MNRAS, 394, 1649 [NASA ADS] [CrossRef] [Google Scholar]
 Moskalik, P., Smolec, R., Kolenberg, K., et al. 2015, MNRAS, 447, 2348 [NASA ADS] [CrossRef] [Google Scholar]
 Neilson, H. R., & Ignace, R. 2014, A&A, 563, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Netzel, H., Smolec, R., & Moskalik, P. 2015, MNRAS, 447, 1173 [NASA ADS] [CrossRef] [Google Scholar]
 Percy, J. R., & Kim, R. Y. H. 2014, J. Am. Association of Variable Star Observers (JAAVSO), 42, 267 [Google Scholar]
 Pietrukowicz, P. 2001, Acta Astron., 51, 247 [NASA ADS] [Google Scholar]
 Pietrukowicz, P. 2002, Acta Astron., 52, 177 [NASA ADS] [Google Scholar]
 Pojmanski, G. 2002, Acta Astron., 52, 397 [NASA ADS] [Google Scholar]
 Poleski, R. 2008, Acta Astron., 58, 313 [NASA ADS] [Google Scholar]
 Poretti, E., Le Borgne, J. F., Rainer, M., et al. 2015, MNRAS, 454, 849 [NASA ADS] [CrossRef] [Google Scholar]
 R Core Team 2015, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria [Google Scholar]
 Schwarz, G. 1978, Ann. Stat., 6, 461 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Simon, N. R., & Lee, A. S. 1981, ApJ, 248, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Smolec, R. 2017, MNRAS, 468, 4299 [NASA ADS] [CrossRef] [Google Scholar]
 Smolec, R., & Śniegowska, M. 2016, MNRAS, 458, 3561 [NASA ADS] [CrossRef] [Google Scholar]
 Smolec, R., Moskalik, P., Evans, N. R., Moffat, A. F. J., & Wade, G. A. 2016, ArXiv eprints [arXiv:1612.02970] [Google Scholar]
 Soszynski, I., Poleski, R., Udalski, A., et al. 2008, Acta Astron., 58, 163 [NASA ADS] [Google Scholar]
 Soszynski, I., Poleski, R., Udalski, A., et al. 2010, Acta Astron., 60, 17 [NASA ADS] [Google Scholar]
 Soszyński, I., Udalski, A., Szymański, M. K., et al. 2015a, Acta Astron., 65, 329 [NASA ADS] [Google Scholar]
 Soszyński, I., Udalski, A., Szymański, M. K., et al. 2015b, Acta Astron., 65, 297 [NASA ADS] [Google Scholar]
 Stothers, R. B. 2009, ApJ, 696, L37 [NASA ADS] [CrossRef] [Google Scholar]
 Süveges, M. 2014, MNRAS, 440, 2099 [NASA ADS] [CrossRef] [Google Scholar]
 Süveges, M., Anderson, R. I., & Eyer, L. 2016, IAU Focus Meeting, 29, 505 [Google Scholar]
 Szabados, L. 1983, Ap&SS, 96, 185 [NASA ADS] [CrossRef] [Google Scholar]
 Szabados, L. 1989, Commmunications of the Konkoly Observatory Hungary, 94, 1 [Google Scholar]
 Szabados, L. 1991, Commmunications of the Konkoly Observatory Hungary, 96, 123 [Google Scholar]
 Szabados, L., Anderson, R. I., Derekas, A., et al. 2013, MNRAS, 434, 870 [NASA ADS] [CrossRef] [Google Scholar]
 Szabó, R., Kolláth, Z., Molnár, L., et al. 2010, MNRAS, 409, 1244 [NASA ADS] [CrossRef] [Google Scholar]
 Turner, D. G., AbdelSabour AbdelLatif, M., & Berdnikov, L. N. 2006, PASP, 118, 410 [NASA ADS] [CrossRef] [Google Scholar]
 Udalski, A., Soszynski, I., Szymanski, M., et al. 1999, Acta Astron., 49, 223 [Google Scholar]
 Udalski, A., Szymanski, M. K., Soszynski, I., & Poleski, R. 2008, Acta Astron., 58, 69 [NASA ADS] [Google Scholar]
 Wozniak, P. R. 2000, Acta Astron., 50, 421 [NASA ADS] [Google Scholar]
 Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
We previously reported on our work in progress at the XXIX IAU General Assembly (Süveges et al. 2016).
All Tables
Basic parameters and their variation for the twinpeak and amplitudechanging Cepheid sample.
All Figures
Fig. 1
Left: CEP1998 as an example of a twinpeak target, for which secondary periods are found after prewhitening the light curve using the primary period. Right: CEP1543 as an example from the control group (no twin peaks). Colors trace observation date, increasing from yellow to blue. Top panels: phasefolded (with primary period) OGLEIII Iband light curve. Center panels: residuals after subtracting Fourier series model. Bottom panels: periodogram of the residuals shown in the center panel. Pink vertical lines at 0, 1, … d^{−1} indicate parasite frequencies due to residual trends in the time series. The solid light blue vertical line indicates the OGLE catalog (primary) frequency. The dashed light blue lines correspond to its ± 1, 2, … aliases. 

In the text 
Fig. 2
Enlarged environment of the removed OGLE catalog frequency and the secondary peak in the secondary periodogram for four twinpeak Cepheids. The OGLE catalog frequency is indicated by a vertical blue line. 

In the text 
Fig. 3
Examples for the basic types of period modulation found among the Cepheids. The plots show the kernelestimated pulsation period (in days) versus Julian Date (in days), plotted as a solid thick black line, together with its bootstrapped pointwise CI (thin black lines; see Sect. 2.2). The heavy orange line denotes the catalog period, which was used as known (nonoptimized) value in the fitted stable reference model for the Cepheid. The orange band indicates the nonparametric bootstrap CI around this period, obtained from the procedure described in Sect. 2.2. The dotted black horizontal and vertical lines are aids to the eye to estimate the extent and time interval of the changes. The gray background highlights time intervals where the deviation from the stable reference model was found significant by a multiple testing procedure (Benjamini & Yekutieli 2001). Leftmost panel: nearlinear trend, middle left: stochastic fluctuations, middle right: quasiperiodic instabilities, rightmost: combination of a trend and weakening quasiperiodic changes. 

In the text 
Fig. 4
Four examples of the effect of the sliding window smoothing on residual periodograms. The GLS periodogram of the residuals from the stable reference model is plotted in black. Twin peaks appear very close to the gray background vertical line that indicates the primary pulsation frequency. Daily aliases of the primary frequency are shown as dashed vertical gray lines. We superposed the periodogram of the residuals from the sliding window fit in red. 

In the text 
Fig. 5
Top panel: residuals of CEP1564 from the constant model, phasefolded with the average primary period. Bottom panel: residuals of CEP1564 from the kernel fits, phasefolded with the secondary period. The typical error bar is shown in the top right corner. 

In the text 
Fig. 6
CEP1748 shown as example of the heuristic model (2) of the changes of the pulsation parameters. Left panel: the loglikelihoods for the pulsation period (in black) and the peaktopeak amplitude (in blue; the loglikelihood functions were shifted to have common minimum at zero). The other two panels show the fits using the frequency yielding the highest likelihood (blue dots) against the sliding window estimates (thick black line), their 95% confidence interval (thin black lines), the stable value from the reference fit (thick orange line) and its confidence band (orange band). The middle panel presents the variations of the pulsation period, the right panel those of the peaktopeak amplitude. 

In the text 
Fig. 7
Distribution of the logarithm of the absolute trend β_{P} (leftmost panel), the frequency of the oscillatory term f_{P} (second panel) and the logarithm of the uncorrected and corrected relative amplitude of the modulation, that is, the amplitudes ΔP and ΔP_{c} of the oscillatory term divided by the catalog period P_{cat} (third and rightmost panels) in model (2), for the changes in the pulsation period of all Cepheids in our sample. Thecolors in the twinpeak groups indicate the absolute value of the separation of the twin peaks, ranging from light blue for small separations to dark blue for large ones. The label T.FU and C.FU corresponds to fundamentalmode twinpeak and control Cepheids, respectively; similarly, the labels T.FO and C.FO indicate first overtone twinpeak and control cepheids, respectively. The T.FO group includes the 5 cepheids selected because of their strong amplitude changes. 

In the text 
Fig. 8
Distribution of the logarithm of the absolute trend (leftmost panel), the frequency of the oscillatory term (second panel) and the logarithm of its relative size, that is, the (uncorrected and corrected) size of the changes ΔA_{1} and ΔA_{1c} divided by Ā_{1}, the first harmonic amplitude of the average stablemodel (third and rightmost panels, respectively), estimated from model (2). Groups and their labels are the same as in Fig. 7. 

In the text 
Fig. 9
Period trend ∣ β_{P} ∣ (left panel) and size of period fluctuations ΔP (right panel) as estimated from model (2) for the 53 Cepheids. Downward pointing triangles: control overtone Cepheids, upward pointing triangles: overtone twinpeak Cepheids, diamonds: fundamentalmode control Cepheids, squares: fundamentalmode twinpeak stars. In the left plot, negative trends are indicated in blue, positive ones in orange. 

In the text 
Fig. 10
Doublelogarithmic periodamplitude diagram for program Cepheids using the original OGLEIII periods and peaktopeak amplitudes. Constant objects are shown as open circles, variable objects as filled circles. Left panel: changing periods. The color corresponds to column T_{P} of Tables 1 and 2, from cyan (0%) to magenta (100%). Right panel: variable amplitudes. The color corresponds to column of Tables 1 and 2, from cyan (0%) to magenta (100%). 

In the text 
Fig. A.1
Weighting scheme. The top panel shows a portion of the observed time series of OGLELMCCEP1621, with the fitted window highlighted by a gray background. The middle panel shows the two weight sequences: in black, the kernel weights, in red, the inverse error weights. The bottom panel exhibits the combined weights used in the fit. 

In the text 
Fig. A.2
Difference of fits with locally constant (red), locally linear (blue) and locally cubic (black) polynomial terms on the example of OGLELMCCEP1833. The top panel shows the estimates of frequencies as a function of time. The middle panel shows the temporal dependence of the estimated mean magnitude. The bottom panel presents the coefficient of the cubic term in the locally cubic model versus time. 

In the text 
Fig. B.1
Simulated trends in pulsation period (top panel) combined with a constant amplitude (bottom panel) as estimated by the localkernel method. Thin red, green and blue solid lines in the top panel: the simulated values using slopes 0, 0.25 and 1 times that of OGLELMCCEP1405 (such that the blue line coincides with the trend component found for this Cepheid). The true amplitude, which had the same constant value in all simulations, is represented by a blue line in the bottom panel. Dotted red, green and blue lines: the pointwise median of the kernel estimates on the 250 repetitions. The colors show a 95% symmetric confidence band around the median. 

In the text 
Fig. B.2
Simulated oscillations of the period, superposed to the period trend estimated for CEP2470. The parameters of the simulated modulations are shown in the upper right corner of the panels. The true simulated modulations of the period are shown with solid red lines. The 0.025, 0.975 quantiles and the median of the kernel estimates are indicated by the colored band and the dotted thick line, similarly to Fig. B.1. The thin black horizontal line in the top panel shows the full width of the sliding window, to ease comparison with the period of the fluctuations. The rugplot added to the xaxis indicates the observation times. 

In the text 
Fig. B.3
Local kernel estimates of the modulation parameters of a Cepheid simulated using the estimated parameters of CEP1748. The notation and the symbols are the same as in Fig. B.2. 

In the text 
Fig. B.4
Local kernel estimates of the modulation parameters of a Cepheid simulated using the estimated parameters of CEP1405. The notation and the symbols are the same as in Fig. B.2. 

In the text 
Fig. C.1
Modulation frequency (top panel) and the ratio between the modulation amplitude and the true amplitude (bottom panel) estimated by model (2) applied to local sliding window estimates of simulated period modulations. Colored dots show the median of the estimates based on 200 different Gaussian noise realizations added to each of five simulated modulated light curves, and vertical segments indicate the range between the 0.025 and 0.975 quantiles. Cyan: using parameters of CEP1527, green: CEP2217, blue: CEP2191, dark red: CEP1140, orange: CEP1833. All simulations used the same sequence of true modulation frequencies. For visibility, the dots based on different Cepheids are horizontally offset from the true value. In the top panel, the solid black line indicates x = y, i.e., the location of accurate estimates, with short horizontal segments aiding the eye as to where horizontally offset points should lie. The dashed line indicates the yearly alias of the x = y line. In the bottom panel, the horizontal gray line indicates y = 1, the location of the ratios in case of unbiased estimates. The solid black line is the best linear fit using the five medians at each modulation frequency corresponding to the seven lowest true modulation frequencies. 

In the text 
Fig. C.2
The modulation frequency (top panel) and the ratio between the modulation amplitude and the true amplitude (bottom panel) estimated by model (2) applied to local sliding window estimates of simulated amplitude modulations. The notation is the same as in Fig. C.1 The solid black line is the best linear fit using the five mediansat each modulation frequency corresponding to the five lowest true modulation frequencies. 

In the text 
Fig. C.3
Local kernel estimates of the period modulation in simulations of a modulation of f_{P} = 0.00137 c∕d, using the stable parameters of CEP2217 (top panel), CEP1140 (middle panel), and CEP1833 (bottom panel). Solid dark red line: the true signal, dotted dark red line: median of the local kernel estimates on 200 random realizations, orange band: the stripe between the pointwise 0.025 and 0.975 quantiles of the estimates. 

In the text 
Fig. D.1
Temporal variation of the pulsation period of the twinpeak Cepheids (solid thick black line), together with its bootstrapped pointwise CI (thin black lines). Heavy orange line: catalog period, orange band: 95% nonparametric bootstrap CI around this period based on the bootstrap procedureof Sect. A.2. The dotted black vertical lines indicate years after HJD − 2 450 000. The gray background highlights time intervals where the deviation from the stable reference model was found significant (cf. Sect. A.3). The times of the observations are indicated with a rugplot at the bottom of the panels. 

In the text 
Fig. D.2
Temporal variation of the pulsation period of the control Cepheids. For the notation, see the caption of Fig. D.1. 

In the text 
Fig. D.3
Temporal variation of the leading harmonic amplitude A_{1} of the twinpeak Cepheids. For the notation, see the caption of Fig. D.2. 

In the text 
Fig. D.4
Temporal variation of the leading harmonic amplitude A_{1} of the control Cepheids. For the notation, see the caption of Fig. D.2. 

In the text 
Fig. D.5
Kernelestimated light curve shapes in two little or nonoverlapping windows for twinpeak fundamental Cepheids, inwhich the peaktopeak amplitude or the amplitude of the leading harmonic term was very different. Fitted curves are indicated by the lines of different colors, and the observations in the different windows, by the different colors and shapes of the symbols. The times of the window centers are given in the legend. The observations close to the window center and therefore more influential in the fit have darker red or blue colors, whereas those in the wings are shown in lighter shades. Approximate bootstrap confidence intervals for the curves are given as light red and light blue stripes around the estimated lines (Sect. 2.2 and Appendix A). 

In the text 
Fig. D.6
Light curve shapes in two different windows for twinpeak overtone Cepheids. For the notation, see the caption of Fig. D.5. 

In the text 
Fig. D.7
Light curve shapes in two different windows for overtone Cepheids showing strong amplitude changes. For the notation, see the caption of Fig. D.5. 

In the text 
Fig. D.8
Light curve shapes in two different windows for control fundamentalmode Cepheids. For the notation, see the caption of Fig. D.5. 

In the text 
Fig. D.9
Light curve shapes in two different windows for control overtone Cepheids. For the notation, see the caption of Fig. D.5. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.