Testing the Stothers model

Variable turbulent convection as the cause of the Blazhko effect – testing the Stothers model


The amplitude and phase modulation observed in a significant fraction of the RR Lyrae variables – the Blazhko effect – represents a long-standing enigma in stellar pulsation theory. No satisfactory explanation for the Blazhko effect has been proposed so far. In this paper we focus on the Stothers (2006) idea, in which modulation is caused by changes in the structure of the outer convective zone, caused by a quasi-periodically changing magnetic field. However, up to this date no quantitative estimates were made to investigate whether such a mechanism can be operational and whether it is capable of reproducing the light variation we observe in Blazhko variables. We address the latter problem. We use a simplified model, in which the variation of turbulent convection is introduced into the non-linear hydrodynamic models in an ad hoc way, neglecting interaction with the magnetic field. We study the light curve variation through the modulation cycle and properties of the resulting frequency spectra. Our results are compared with Kepler observations of RR Lyr. We find that reproducing the light curve variation, as is observed in RR Lyr, requires a huge modulation of the mixing length, of the order of 50 per cent, on a relatively short time-scale of less than 40 days. Even then, we are not able to reproduce neither all the observed relations between modulation components present in the frequency spectrum, nor the relations between Fourier parameters describing the shape of the instantaneous light curves.

convection – hydrodynamics – methods: numerical – stars: oscillations – stars: variables: RR Lyrae – stars: individual: RR Lyr

1 Introduction

The periodic amplitude and phase modulation observed in many RR Lyrae variables – the Blazhko effect – is one of the most important, still unsolved problems in classical pulsation theory. In the last few years, thanks to extensive ground-based observation campaigns (e.g., Kolenberg & Guggenberger, 2007; Jurcsik et al., 2009) and satellite missions, CoRoT (Chadid et al., 2010) and Kepler (e.g., Benkő et al., 2010), we finally obtained nearly continuous data, allowing for detailed study of the light variation associated with the Blazhko cycle. One of the most exciting new findings is the period doubling phenomenon (Kolenberg et al., 2010a), never detected from ground-based data. It occurs at some phases of the Blazhko cycle and manifests in an alternating shape of the light variation in consecutive pulsation cycles (see also Szabó et al., 2010). It was not detected in any non-modulated RR Lyrae star so far. Satellite data allow to study the light curve changes over the Blazhko cycle in detail. Clearly, consecutive Blazhko cycles are not exactly repetitive (Kolenberg et al., 2011). Also, the period of the Blazhko modulation differs from cycle to cycle. In the frequency spectra, the Blazhko phenomenon manifests itself as equidistantly spaced multiplets around main pulsation frequency and its harmonics. Triplets and quintuplets are detected in ground-based observations (see e.g., Jurcsik et al., 2008; Kolenberg et al., 2009). In satellite data, not only triplets and quintuplets are visible, but also tenth-order side peaks (Chadid et al., 2010).

All the newly discovered features of the Blazhko effect put stringent constraints on the models proposed to explain this longstanding enigma. In fact, the two most elaborated models, the Magnetic Oblique Rotator/Pulsator model (MORP, Shibahashi, 2000), and the Non-radial Resonant Rotator/Pulsator model (NRRP, e.g., Van Hoolst et al., 1998; Nowakowski & Dziembowski, 2001; Dziembowski & Mizerski, 2004) are ruled out by observations for several reasons. First, strong dipole magnetic fields, necessary for the MORP model to work, were not detected in RR Lyrae stars (Chadid et al., 2004; Kolenberg & Bagnulo, 2009). Secondly, both the MORP and NRRP models predict that the observed light curve modulation should manifest through specific features in the frequency spectra. The dipole geometry of the magnetic field, considered in the MORP model, leads to a quintuplet structure in the frequency spectra. A more complicated field geometry, that could possibly lead to higher order modulation components was not considered so far. Within the NRRP model it is shown that the modes of are most easily excited, giving rise to a triplet structure in the frequency spectra. The excitation of higher order nonradial modes is less likely. In addition, due to cancellation effects, the expected amplitudes would be very low for high-order modes. To the contrary, in satellite data we clearly detect higher-order modulation components. Chadid et al. (2010) detected modulation components up to the tenth order. In both the MORP and NRRP models, strong asymmetries in the amplitudes of the side peaks, as observed in the majority of Blazhko stars, are not trivial to reproduce. Considering the triplets, in about 75 per cent of the Blazhko variables, the higher frequency side peak has a higher amplitude than the lower frequency side peak (Alcock et al., 2003). Third, both models propose mechanisms that imply a clockwork regularity, predicting robust modulation periods, e.g., equal to the rotation period of the star, and repeatable Blazhko cycles. Neither is seen in the data. The Blazhko variation can change considerably even from cycle to cycle (Kolenberg et al., 2011). In several stars two modulation periods are apparent (e.g., Detre & Szeidl, 1973; LaCluyzé et al., 2004; Szczygiel & Fabrycky, 2007). Hence, the connection between the period of the Blazhko modulation and the rotation period is highly questionable. Also, in several cases consecutive Blazhko cycles differ significantly, which manifests, e.g., in systematic growth or decay of the modulation amplitude.

We also mention the recent studies of the Blazhko phenomenon by Buchler & Kolláth (2011) and Jurcsik et al. (2011). Using the amplitude equation formalism (see e.g., Buchler & Goupil, 1984), Buchler & Kolláth (2011) showed that the high order half-integer resonance, proposed by Szabó et al. (2010) to explain the period doubling, can also cause the phase and amplitude modulation as observed in Blazhko variables. The result is very exciting and a more detailed analysis and confirmation through the hydrodynamic modelling would be of great value. Jurcsik et al. (2011) studied the Blazhko variables in the globular cluster M5. They showed that the Blazhko stars tend to be situated on the zero-age horizontal branch and at the blue edge of the fundamental mode instability strip. They speculate that this location hints that the Blazhko effect may be connected to the mode switch from the fundamental mode to the first overtone pulsation during the evolution.

Very recently, an idea proposed by Stothers (2006, 2010), which connects the Blazhko effect with variable magnetic stellar cycles, has gained popularity. In this idea, the Blazhko modulation is connected to the cyclic strengthening and weakening of turbulent convection in the outer stellar layers, caused by the postulated transient magnetic field. The field decays cyclically and is subsequently built up anew by the turbulent/rotational dynamo. Details of the process, particularly the interaction between magnetic filed, turbulent convection and pulsation were not discussed. Also, numerical estimates, e.g., on the expected strength of modulation or the necessary strength of the variable magnetic field, were not presented. For a critical analysis of the Stothers idea we refer the reader to Kovács (2009).

In this paper we use our pulsation hydrocodes to investigate whether variable turbulent convection may cause such a light variation as is associated with the Blazhko effect. The variation of turbulent convection is put into the hydrodynamical models in an ad hoc way – interaction with the magnetic field is neglected. Currently, there are no models capable of reproducing a variable magnetic field due to dynamo mechanisms and of describing its dynamical coupling with turbulent convection and highly non-linear pulsation, features we deal with in the case of RR Lyrae variables. As a consequence, our model is strongly simplified and our results will serve for qualitative comparison with the observations. Having this in mind, we state the premise of our paper. If we can reproduce the most important observational constraints, the Stothers idea is certainly worth further, more detailed investigation. But if we do not succeed it needs revision.

The structure of the paper is the following. In Section 2 we present the details of our numerical model testing the idea proposed by Stothers (2006). In Section 3 we describe the properties of our models and compare our results with observations, focusing on overall properties of Blazhko variables and on recent Kepler observations of the famous Blazhko variable, RR Lyr (KIC 7198959) (Kolenberg et al., 2011). Some additional models are discussed in Section 4 and conclusions are presented in Section 5.

2 Numerical methods

The basic tool in our modelling is the Warsaw non-linear convective pulsation hydrocode (Smolec & Moskalik, 2008a), which we briefly describe in Section 2.1. The code is slightly modified in order to model the effects of variable turbulent convection. Details, as well as the modelling procedure, are outlined in Section 2.2. In Section 2.3 we present the computed sequences of models to be analysed in this paper.

2.1 Non-linear convective hydrocodes

In all our computations we use pulsation codes described by Smolec & Moskalik (2008a). These are a static model-builder, linear non-adiabatic code and a direct time integration non-linear hydrocode. For convective energy transfer we use the Kuhfuß (1986) convection model reformulated for the use in stellar pulsation codes. Radiation is described in the diffusion approximation. The codes use a simple Lagrangian mesh.

For an extensive description and details of numerical implementation we refer the reader to Smolec & Moskalik (2008a). Here we note that the model equations contain eight order-of-unity scaling parameters that describe the turbulent convection model. These are mixing-length, , and parameters multiplying the turbulent fluxes and terms that drive/damp the turbulent energy, (turbulent pressure), (eddy-viscous dissipation), (convective heat flux), (kinetic turbulent energy flux), (buoyant driving), (turbulent dissipation) and (radiative cooling). Theory provides no guidance for their values. However, some standard values are in use, which result from comparison of a static, time-independent version of the model with the standard mixing-length theory (see Wuchterl & Feuchtinger, 1998; Smolec & Moskalik, 2008a). In practice, values of these parameters should be such that models satisfy as many observational constraints as possible.

We also stress that we use the original Kuhfuß (1986) prescription, in which essential buoyancy effects are included in convectively stable regions of the model (negative buoyancy). We note that the neglect of negative buoyancy (as is done, e.g., in the Florida-Budapest code, Kolláth et al., 2002) can lead to unphysical effects in the computed models. For an extensive discussion on the subject we refer the reader to, e.g., Smolec & Moskalik (2008b) or Smolec (2009).

2.2 Exploring Stothers’ idea

Stothers (2006) proposed that the Blazhko modulation is connected to the postulated cyclic strengthening and weakening of turbulent convection in the outer stellar layers, caused by a transient magnetic field. Turbulent convection becomes more vigorous during the decay of the magnetic field and it is quenched as the magnetic field builds up anew. Details of the underlying processes, particularly the interaction between magnetic field, turbulent convection and pulsation were not elaborated by Stothers.

Our code neglects the effects of magnetic fields and as such is not suitable for modelling the dynamical coupling of pulsation and turbulent convection on the one hand, and the transient magnetic fields postulated by Stothers (2006) on the other hand. We only assume that the strength of the turbulent convection varies in time and do not elaborate on the physical mechanism behind such changes. The strength of turbulent convection is changed by cyclic variation of one of the -parameters entering our model (see Section 2.1). The mixing length parameter, , is our first choice, as it regulates the overall strength of convection and affects all the turbulent quantities entering the model (see equations in e.g., Smolec & Moskalik, 2008a). We assume that the mixing length is either a sinusoidal or a triangular function of time, as illustrated in Fig. 1.

Figure 1: Possible shapes of mixing-length modulation considered in this study.

The initial steps in our modelling correspond to a standard hydrodynamic model computation, without modulation of the turbulent convection (constant mixing-length). First, we construct a static equilibrium model. Next, we compute the non-linear model. The static model is perturbed with a scaled linear velocity eigenfunction followed by time evolution of the model. The model integration is stopped, when the full-amplitude single-periodic pulsation (the limit cycle) is reached. Then we start to modulate the turbulent convection. The non-linear model integration is restarted with the mixing-length parameter varying in time according to the chosen functional form (see Fig. 1). The model is integrated for several Blazhko cycles. After a few initial cycles, the consecutive ones are almost indistinguishable from one another, indicating that the Blazhko limit cycle is reached. Then, the model integration is stopped and the resulting light variation is subjected to a detailed analysis (Section 3). 3

In our computations we start to modulate the turbulent convection at the phase of maximum radius – the initial phase of modulation relative to the unperturbed mono-mode pulsation is equal to zero. To check whether the choice of the initial phase affects the results, we have computed several additional models with different initial phases: 0.25, 0.50 and 0.75, and starting the modulation of turbulent convection with an either increasing or decreasing mixing-length parameter. The final pulsation state (the Blazhko limit cycle) is insensitive to the initial phase. The computed light curves are almost indistinguishable for the models with different initial phases (the trajectories overlap in the Fourier parameter plots, Figs. 35). The amplitudes of the peaks in the frequency spectra (see Section 3.4) are also identical to within a fraction of a per cent.

In Fig. 2 we present the light variation for a typical model over slightly more than one Blazhko cycle. The modulation of the pulsation amplitude is clearly visible. In the lower panel we present the close-up of maximum amplitude phases at which the period doubling occurs.4

Figure 2: Bolometric magnitude versus time for a particular model of set M6 (). In the lower panel we show the close-up of phases at which period doubling occurs.

2.3 Computed model sequences

In our modelling we focus on the fundamental mode models only, as most Blazhko variables are fundamental mode pulsators. The chemical composition of all our models is the same, and , which is typical for RR Lyrae stars. We note that the metallicities of the Blazhko variables do not differ from those of the non-modulated RR Lyrae pulsators (Smolec, 2005). In all model computations we use OP opacities (Seaton, 2005) at low temperatures supplemented with the Alexander & Ferguson (1994) opacity data. Opacities were computed for the solar mixture of Asplund et al. (2004). Other physical parameters of the initial, non-modulated models (masses and luminosities) are collected in Table 1. Parameters of set A, which is the most explored set in this study, are close to those quoted for RR Lyr (Kolenberg et al., 2010b). In sets B, C and D the luminosity is varied, which allows to study the effects of a different -ratio.

A 0.65 50.0
B 0.65 40.0
C 0.65 60.0
D 0.65 70.0
Table 1: Physical parameters of the computed sequences of initial models. The chemical composition of all our models is the same, and . Effective temperatures are in the range of  K.

The convective parameters of the initial non-modulated models are collected in Table 2. We adopt only one set of convective parameters, very similar to the one we have adopted in Baranowski et al. (2009) and Smolec (2009). With these convective parameters, overall properties of the Galactic Cepheids, both fundamental mode and first overtone pulsators, were modelled successfully. Here, we only increased the eddy-viscous dissipation ( parameter) in order to match the typical pulsation amplitudes of the RR Lyrae stars. We note that the modelling of RR Lyrae light curves is a difficult problem and some discrepancies with the observations still remain (see e.g., the detailed comparison done by Kovács & Kanbur (1998) or Feuchtinger (1999)). Nevertheless, the general properties of the RR Lyrae light curves are quite well reproduced with the current convective hydrocodes. In particular, the light curves of individual stars can be nicely matched with the hydrodynamic models (see e.g., Marconi, 2009; Marconi & Clementini, 2005).

1.5 0.65 1.0 1.0 1.0 0.0 0.0 1.0
Table 2: Convective parameters considered for non-modulated initial models. Parameters , , , and are given in the units of standard values (, , and ; see Smolec & Moskalik, 2008a, for details).

The static models, the first step in our modelling, consist of 150 mass zones extending down to . Fifty outer zones have equal mass, down to the anchor zone in which the temperature is set to (hydrogen ionisation). The mass of the remaining zones increases geometrically inward. In the second step of our modelling procedure, we have computed non-linear full amplitude models. To reach the limit-cycle pulsation, we have computed 2000 pulsation cycles by default. These non-linear non-modulated models were subsequently used as initial models for model integrations with variable turbulent convection.

Several sequences of modulated models were computed. The properties of these models are summarised in Table 3. We investigated the effects of different amplitudes of the turbulent convection modulation, different modulation periods, and different modulation shapes (see Fig. 1). In each set, several models of different effective temperatures, extending through the whole instability strip, were computed. For most of the computed models, the physical parameters of set A were used (Table 1). We note that the mean physical parameters of the modulated models, such as mean radius or mean effective temperature, are affected by the modulation of turbulent convection and differ from the static equilibrium values. In the following, we will refer to the particular models by providing the adequate set name from Table 3 and the effective temperature of the underlying static model, .

Set physical function
M1 A sine 10% 60 d -
M2 A sine 20% 60 d -
M3 A sine 30% 60 d -
M4 A sine 20% 40 d -
M5 A sine 40% 60 d -
M6 A sine 50% 60 d -
M7 A sine 50% 40 d -
M8 A sine 50% 120 d -
M7L4 B sine 50% 40 d -
M7L6 C sine 50% 40 d -
M7L7 D sine 50% 40 d -
MT1 A triang. 50% 60 d 45 d
MT2 A triang. 50% 60 d 15 d
Table 3: Parameters of turbulent convection modulation for the model sequences considered in this study. For definitions of , and see Fig. 1.

3 Results

In this section we analyse the computed model sequences, described in Section 2.3, focusing on the light curve variation through the Blazhko cycle, and on the properties of the frequency spectra. First, in Section 3.1, we describe the overall properties of the light curve modulation in terms of the Fourier decomposition parameters (e.g., Simon & Lee, 1981). The variety of behaviours and trends we compute can serve for comparison with continuous observations of Blazhko variables by satellite missions. In Section 3.2 we compare our results with Kepler observations of RR Lyr, for which a detailed analysis of the light variation was already done (Kolenberg et al., 2011). In Section 3.3 we briefly analyse the pulsation period changes in our models and in Section 3.4 we present the results of the frequency analysis of our models, focusing on the overall properties of the frequency spectra and how these compare with properties of Blazhko variables and RR Lyr in particular. Finally, in Section 3.5 we provide a more detailed discussion on the period doubling phenomenon.

3.1 Light curve variation through the Blazhko cycle

To study the light curve variation through the Blazhko cycle we use the time dependent Fourier analysis (Kovács et al., 1987). We fit the Fourier series to the light variation of the consecutive pulsation cycles,


and analyse time variation of the Fourier decomposition parameters of low order, the amplitude, , amplitude ratio, , and phase difference, ,


We note that the derivation of instantaneous amplitudes and phases through the time-dependent Fourier analysis is justified, as the modulation is slow compared with the period of oscillation.

For the analysed smooth hydrodynamical data, the Fourier parameters also vary smoothly with time. If period doubling occurs, the alternating shapes of the consecutive pulsation cycles manifest in series of wiggles that appear in the run of the Fourier parameters. In Figs. 35 we plot different relations between the Fourier parameters for different models. We present the effects of different amplitudes of the mixing-length modulation (Fig. 3), different modulation periods (Fig. 4) and different modulation shapes (Fig. 5) on the light curve variation. The instantaneous mixing length, , and , are plotted versus the amplitude, , in the top, middle and bottom rows of these figures, respectively. As in our models consecutive Blazhko cycles are repetitive, relations take the form of closed trajectories. In order to visualise the variation across the instability strip, we display the models with different effective temperatures, for which we choose =6300 K (left columns), =6600 K (middle columns; models to be compared later with RR Lyr) and =6800 K (right columns; models with period doubling). All models are characterised by the physical parameters of set A (Table 1). On each trajectory in these figures, a plus sign indicates the phase of maximum mixing length, and the direction of the time flow along the trajectory is depicted schematically. In addition, in each panel crosses connected with dashed line indicate the location of several non-modulated models with fixed mixing-length values, in a range covered by our modulated models. Below we describe the properties of the light curve variation connected with the location within the instability strip and different parameters of the modulation of the turbulent convection.

Figure 3: Light curve variation through the Blazhko cycle for the models with different amplitude of mixing-length modulation, 10 per cent (set M1), 30 per cent (set M3) and 50 per cent (set M6). The instantaneous mixing-length, and are plotted versus in the top, middle and bottom rows of the Figure. In consecutive columns we plot the results for models of different (6300 K, 6600 K and 6800 K). A plus sign on each trajectory indicates the phase of maximum mixing length. The direction of the time flow along each trajectory is shown on a schematic circle. Crosses connected with a long-dashed line correspond to the non-modulated models with different values of the mixing length. For =6300 K and the fundamental mode is linearly stable and the models were not computed. For =6800 K and the models switch into first overtone (1O) pulsation (two models in the Figure).
Figure 4: The same as Fig. 3, but for the models with different periods of mixing-length modulation; 60 days (M6), 40 days (M7) and 120 days (M8).
Figure 5: The same as Fig. 3, but for the models with different shapes of mixing-length modulation; sinusoidal (set M6) and asymmetric triangular (MT1 and MT2).

General properties of the light curve modulation. Variation through the instability strip.

Analysis of Figs. 35 reveals some general properties of the computed trajectories, independent of the parameters of turbulent convection modulation (amplitude, period and shape). A basic observation is that the higher the effective temperature is, , the higher is the mean pulsation amplitude, . Also, the higher the mean pulsation amplitude, the smaller the range of variation of and . For the models with =6300 K we are close to the linear red edge of the fundamental mode. In fact, for the models with a strong modulation of the turbulent convection (e.g., of the sets M3 and M6 in Fig. 3), phases of large instantaneous mixing-length values () correspond to the pulsationally stable equilibrium models. On the other hand, at higher temperatures (=6800 K), small values of the mixing length () at some modulation phases, if adopted in non-modulated models, would lead to first overtone (1O) pulsation (see the long-dashed lines in the figures). It is clear that the temporary light curve of the model in which turbulent convection is modulated is significantly different from that of the non-modulated model adopting the same value of the mixing length. In addition, hysteresis is clearly visible for the modulated models. Very different shapes of the light curve are possible at two phases characterised by the same instantaneous mixing length. The described features are in qualitative agreement with the observations. Jurcsik, Benkö & Szeidl (2002) noted that the light curves of Blazhko stars are never (at any phase of the Blazhko cycle) like those of non-Blazhko stars, thus always distorted. A detailed analysis of the light curve variation in Kepler RR Lyr data (Kolenberg et al., 2011) reveals a similar behaviour of the Fourier parameters as plotted in the figures. A more quantitative comparison is postponed to the next section.

Another interesting feature is visible in the vs. plots (top rows of Figs. 35). We note that the minimum value of the amplitude does not coincide with the maximum value of the mixing-length parameter, as one might naively expect, but it is delayed. The reason for such a delay is not clear. It shows the inertia of the dynamical system, which does not adjust instantaneously to the changing dissipation.

Considering the relations vs. and vs. (middle and bottom rows of Figs. 35) we see that for low temperature models (of small mean pulsation amplitudes) trajectories have the shape of double loops, while at higher temperatures (and high mean pulsation amplitudes) single loops are present. The direction of trajectories (in case of double loops, the direction of the larger loop) is always clockwise for vs. trajectory (bottom rows of Figs. 35). For vs. trajectories (middle rows of Figs. 35) the direction is either counter-clockwise (for cooler models, =6300 K and =6600 K) or clockwise (for the hottest models, =6800 K). The analysis of additional models across the instability strip reveals that the transition occurs at around =6700 K. For these models we deal with an almost single-valued dependence of on . We note that the direction of the vs. loop is related to a particular asymmetry of the triplet components in the frequency spectra (Szeidl & Jurcsik, 2009). We discuss this problem in more detail in Section 3.4.2.

Considering the vs. relation, we first note that intuitively we expect that should correlate with . The higher the amplitude, , the larger the non-linearity and consequently larger the contribution of the harmonic terms (and thus, higher ). In our models, in which both and are functions of time, the relation between and is complicated. For the hottest models, of higher mean pulsation amplitude, correlates with . For the lower temperature models, of lower mean pulsation amplitudes, the trajectories are more extended (”blown-up”) and the relation is not obvious. At some phases in the Blazhko cycle, a high amplitude is accompanied by lower values of .

Period doubling is clearly visible for the highest temperature models ( K) displayed in Figs. 35, in which the mixing-length modulation is strong (Fig. 3). During the phases with period doubling, the trajectories are no longer smooth. The series of wiggles is present, which indicates that the consecutive pulsation cycles alternate. Period doubling will be discussed in more detail in Section 3.5.

In addition to the models running horizontally in the instability strip, we have computed some models with different ratio. Four sequences of models were computed, M7L4, M7, M7L6 and M7L7 (see Table 3), in which the luminosities vary from to . The mass is the same in all these model sequences () and was varied. In all these sequences the amplitude of the mixing-length modulation was 50 per cent, the modulation period 40 days, and the shape of the modulation was sinusoidal. Except for the shifts in the computed trajectories, connected with the different luminosities of the initial models, the trajectories are very similar and hence not presented in a separate figure. We have noted that the direction of the vs. trajectories for high temperature models depends on the luminosity and is clockwise for the low-luminosity models ( and ) and counter-clockwise for the high-luminosity models ( and ), in which displays a flat dependence on . Consequently, a clockwise direction of the vs. trajectory is present only in higher-temperature and lower-luminosity models.

Different strength of the turbulent convection modulation.

In Fig. 3 we plot the results for the models with a different amplitude of the turbulent convection modulation, per cent (M1), per cent (M3) and per cent (M6). The period and shape of the modulation are the same for these three sets (, sine). As expected, the stronger the modulation, the more complicated (‘non-linear’) trajectories we get and higher the ranges of variation of the Fourier parameters. This property of the computed trajectories will be used in the next section to estimate the required amplitude of the mixing-length modulation to reproduce the RR Lyr observations. We also note that, depending on the strength of the modulation, the mean values of Fourier parameters vary. This concerns the amplitude, , and the amplitude ratio, . The mean phase difference is not very sensitive to the strength of the modulation. For the most convective models of =6300 K, the stronger the modulation, the smaller the amplitude and the amplitude ratio.

Different period of the turbulent convection modulation.

In Fig. 4 we plot the results for models with different period of turbulent convection modulation, (M7), (M6) and (M8). The amplitude and shape of the modulation are the same for these three sets ( per cent, sine). The basic observation is that for longer modulation period, the computed loops are larger (more ”blown-up”), and thus, the larger the ranges of variation of Fourier parameters, but the effect is not as strong as in case of different amplitudes of turbulent convection modulation. We also note that the range of effective temperatures in which the period doubling occurs, depends on the modulation period, but we postpone the detailed discussion to Section 3.5.

Different shape of the turbulent convection modulation.

In Fig. 5 we plot the results for models of sets M6 (sinusoidal modulation of the mixing length), MT1 and MT2 (asymmetric triangular modulation). The period and amplitude of mixing-length modulation are the same for these three sets (, per cent). We conclude that the differences between the trajectories are rather minute and the light curve variation is very similar for the different shapes of modulation.

3.2 Comparison with RR Lyr

RR Lyr is not only the prototype of the RR Lyrae variables but also a famous Blazhko variable with a strongly modulated light curve. Its Blazhko period is subject to small variations, its present value being  days (Kolenberg et al., 2011). Thanks to its location in the Kepler field of view, we now have excellent, nearly continuous photometric data covering slightly more than three Blazhko periods – seasons Q1+Q2 of the Kepler long cadence data (Jenkins et al., 2010a, b) of which Q1 is already public5.

In this section we compare the light curve variation through the Blazhko cycle, as it is observed in RR Lyr, with our models. The Kepler photometric passband is rather wide, covering the combined range of the standard V and R passbands (Koch et al., 2010). It justifies the direct comparison of our model bolometric light curves with the Kepler RR Lyr light curve. This simplification does not affect the estimates presented in this section. We note that Nemec et al. (2011) derived the transformation between the Fourier parameters in Kp and in V passbands based on the photometry of three non-modulated RR Lyrae stars in both bands. The differences are systematic but very small.

In Fig. 6 we plot the vs. and vs. relations as observed for RR Lyr (season Q2 of Kepler data). The ranges of variation of the Fourier parameters are rather large. In our models, as noted in the previous section, the larger the amplitude of turbulent convection modulation is, the larger is the range of variation of the Fourier parameters. Now, we can estimate the strength of mixing-length modulation necessary to explain the RR Lyr observations, through comparing the model and observed ranges of variation of , and . To this aim, for each Fourier parameter , , we define the parameter ,


constructed using the minimum, , and maximum, , values of during the modulation cycle. We note that the value of is independent of simple scaling of the Kepler data by a constant factor. For RR Lyr we have , and (Kolenberg et al., 2011, Fig. 6). The computed values for the models with a different strength of the turbulent convection modulation, per cent (set M1), per cent (set M3), and per cent (sets M6 and M7) are collected in Table 4. For each set, the values for seven models of different are computed. The values that agree within 20 per cent with RR Lyr values are marked with an asterisk. Note that the sets M6 and M7 have the same amplitude of the mixing-length modulation ( per cent), but for set M7 the modulation period is shorter (40 days, very close to RR Lyr’s modulation period).

Figure 6: Fourier parameters, and , plotted vs. for Kepler RR Lyr data (season Q2) and for two models of set M3.
M1 (A=10%) M3 (A=30%) M6 (A=50%) M7 (A=50%)
6300K 0.16 0.26 *0.25 *0.51 1.05 0.53 0.89 1.47 0.64 0.79 1.41 0.65
6400K 0.16 0.27 0.17 *0.53 0.78 0.39 0.84 1.08 0.51 *0.74 0.99 0.50
6500K 0.17 0.19 0.10 0.48 0.48 *0.30 0.76 0.73 0.41 *0.66 0.65 0.41
6600K 0.14 0.06 0.07 0.42 0.19 *0.23 *0.64 *0.35 0.37 *0.57 *0.30 0.38
6700K 0.11 0.02 0.07 0.33 0.08 *0.22 *0.53 0.19 0.39 0.46 0.17 0.40
6800K 0.09 0.04 0.08 0.27 0.16 *0.23 0.43 0.28 0.37 0.37 0.27 0.38
6900K 0.07 0.06 0.07 0.21 0.17 *0.22 0.34 0.28 0.36 0.32 0.28 0.36
RR Lyr 0.62 0.38 0.27 0.62 0.38 0.27 0.62 0.38 0.27 0.62 0.38 0.27
Table 4: Ranges of variation of the Fourier parameters for the different models of sets M1, M3, M6 and M7 compared with the RR Lyr Kepler data. In the first column effective temperature of the initial non-modulated model, , is given. In the consecutive columns the ranges of variation of Fourier parameters, , and (see eq. 3) are given for different sets of modulated models. In the last row, the data for RR Lyr are provided for reference. Model values that agree with the corresponding RR Lyr value within 20 per cent are marked with an asterisk.

To reproduce the ranges of the Fourier parameter variation in RR Lyr, a mixing-length modulation with an amplitude equal to at least 30 per cent is necessary. The best match is found for an amplitude of the mixing-length modulation equal to 50 per cent (sets M6 and M7). Then, both and can be matched for the model with =6600 K. For models with =6600 K, the mean pulsation period ( d) matches RR Lyr’s period (0.567 d, Kolenberg et al., 2011) almost exactly. The range of the phase variation, , however, is slightly higher for these models than is observed in RR Lyr.

The most stringent constraints on the Stothers model come from highly modulated stars, like RR Lyr. Tiny modulations as observed for many Blazhko stars can be easily reproduced assuming a small amplitude of turbulent convection modulation in our models. It is now evident that, if the mechanism proposed by Stothers is responsible for the Blazhko effect, the observed large modulations of the light curves (as we observe in RR Lyr) require significant changes of the mixing-length over the Blazhko cycle. Note that the modulation amplitude, , as defined in Fig. 1, is actually a semi-amplitude. For the considered models, per cent means that the mixing length, , varies in a huge range, from 0.75 to 2.25. In our opinion, such large changes in the effectiveness of the turbulent convection on relatively short time-scales of typically tens to hundreds of days (typical Blazhko periods) are highly unlikely (see discussion in Section 5).

Now we focus on the shape of the Fourier parameter trajectories. A close inspection of Figs. 35, particularly the models with =6600 K (for which the period of the fundamental mode matches that of RR Lyr) shows that the model vs. loops are quite similar to the RR Lyr loop (top panel of Fig. 6). What we typically see in our models (around =6600 K) is a larger (”blown-up”) part of the loops on the left side and cusps on the right side. This is what we observe in RR Lyr, however, the cusp in our models is located at higher values of compared to RR Lyr. A comparison of Kepler RR Lyr data with two models of set M3 (for which we get the best match) is shown in the upper panel of Fig. 6. We note that also the direction of the time flow is the same for both RR Lyr and for the models (clockwise).

Considering the vs. loops (bottom panel of Fig. 6) we cannot reproduce the behaviour that we observe in RR Lyr. In RR Lyr anticorrelates with through the majority of the Blazhko cycle. This is not reproduced by our models, although at some phases of the modulation increases as decreases. The bottom panel of Fig. 6 provides a direct comparison of the model (set M3) and RR Lyr loops. Also here, the direction of the time flow is the same for our models and for RR Lyr (bigger loop, counter-clockwise).

3.3 Changes of the pulsation period

In this section we analyse the changes of the pulsation period. In the mathematical description, a period change is equivalent to a change of the pulsation phase. From the observational point of view, the two are indistinguishable. It is the physical interpretation where the difference occurs. Period changes are caused by the overall changes of the stellar structure, while phase changes may result, e.g., from the nonlinear interaction of pulsation modes (which does not affect the stellar structure).

In our models, due to the changes of the convective structure, the pulsation period changes during the modulation. The amplitude of the period variation, , is most sensitive to the amplitude of turbulent convection modulation, and hence can be used to estimate the required strength of mixing-length modulation in a similar way we have done in the previous section. In Table 5 we list the period changes for the models with different amplitudes of the turbulent convection modulation, per cent (set M1), per cent (set M3) and per cent (set M6). The period changes were estimated based on the radius variation of our hydrodynamic models, using the time difference between consecutive radius maxima as an estimate of the instantaneous period. Then, the full amplitude of the period variation was used to derive values. It is evident that the larger the amplitude of mixing-length modulation is, the larger the amplitude of period variation. Also, the trend of decreasing amplitude of the period variation with increasing effective temperature is clear. In Blazhko RR Lyrae variables, the pulsation period changes are clearly detected, with a typical amplitude, , between 0.2 and 1.4 per cent (Molnár & Kolláth, 2010). For RR Lyr, the period change is 0.83 per cent (Kolenberg et al., 2011). Comparison with values collected in Table 5 indicates that in order to reproduce such period variations, large amplitudes of the mixing-length modulation are necessary, at least of the order of 50 per cent, in agreement with the value derived in the previous section and in agreement with estimate of Molnár & Kolláth (2010), based on the analysis of the periods of linear equilibrium models.

[per cent]
M1 M3 M6
6300K 0.12 0.39 0.97
6400K 0.11 0.34 0.88
6500K 0.092 0.31 0.77
6600K 0.095 0.30 0.64
6700K 0.051 0.28 0.51
6800K 0.11 0.27 0.44
6900K 0.10 0.20 0.49
Table 5: Amplitudes of the period change, , for several models of different and different amplitudes of turbulent convection modulation, 10 per cent (M1), 30 per cent (set M3) and 50 per cent (set M6).

3.4 Analysis of frequency spectra

In this section we focus on the frequency analysis of our models, for which we use the Period04 package (Lenz & Breger, 2005). Differently from Section 3.1, we now analyse the light variation over many pulsation periods and many Blazhko cycles. At this long time-scale the system is stationary and we can fit the data with a Fourier sum of the following form:


In the above formula, the first line corresponds to the main pulsation frequency, , and its harmonics, , the second and third lines correspond to the components of the triplets, , the next two lines describe other higher-order multiplet components and in the last line we account for the modulation frequency, , and its harmonics, . We assume that the components of the multiplets are equidistant. Also the modulation frequency, , is known, as it is equal to the inverse of the assumed period of mixing-length modulation (see Table 3). Consequently, only one frequency is determined from the data, the main pulsation frequency, which is not known a priori (the non-linear period differs from the linear one). By default we use , and for all the models (note that signal at is strong in our models).

The length of the hydrodynamic data we analyse corresponds to four full Blazhko cycles. To speed up the computations, the sampling of the model light curve is degraded to roughly 60 points per pulsation period (twice the Kepler resolution for RR Lyr).

Below, we present some details of our analysis for one particular model, and later (Section 3.4.2) we discuss the general properties of the frequency spectra of our models and compare them with the observations.

Frequency spectra analysis – particular case

We analyse one particular model from set M6 with . The model displays a clear period doubling in the computed light curve (see, e.g., Fig. 3). The mean pulsation period for this model is , which is not far from the RR Lyr period ( d, Kolenberg et al., 2011). As noted before, the best match with RR Lyr’s period is obtained for the model with =6600 K, but this model does not display the period doubling.

For the discussed model, the variation of the bolometric magnitude with time for slightly more than one Blazhko cycle (70 days) is shown in the upper panel of Fig. 2. Period doubling is clearly visible during the phases of maximum pulsation amplitude. The lower panel of Fig. 2 shows these phases. The period doubling is obvious and lasts for several pulsation cycles. However, its amplitude is not very large.

In Fig. 7 we plot some results of the frequency analysis. The upper panel of Fig. 7 shows the close-up around the fundamental mode frequency, , and its harmonic, , after prewhitening the data with the fundamental mode frequency and its harmonics (dashed lines), and five consecutive components of the multiplet structure. In the middle panel of Fig. 7, we show the vicinity of the fundamental mode frequency. All removed frequencies are marked with dashed lines. Clearly, many more higher-order side peaks are visible. Residual power at the position of corresponds to the long-term evolution of the model toward the final limiting cycle pulsation (which is of exponential character). The bottom panel of Fig. 7 shows the vicinity of the half-integer frequency (HIF), . The signal at this frequency is a signature of the period doubling (see Szabó et al., 2010), which appears in the discussed model during a short fraction of the modulation cycle. Consequently, the signal is strongly modulated with the Blazhko period and hence many equally-spaced peaks are clearly visible, with the separation corresponding to the modulation frequency. The envelope of the signal located at the HIF’s is very regular and resembles a Gaussian. The width of this Gaussian corresponds to the duration of the period doubling behaviour observed in the model ( days).

Figure 7: Frequency spectrum for the particular model of set M6 () after prewhitening with the Fourier sum described by Eq. (4) (upper panel). Two close-ups are shown: at the location of fundamental mode frequency (middle panel) and at its subharmonic, (bottom panel). Removed frequencies are marked with dashed lines.

Our analysis is focused on the first-order side peaks, the triplets. The relations between their amplitudes, , and (see Eq. 4) were studied for several Blazhko stars observed from the ground (e.g., Jurcsik et al., 2006) as well as for the Kepler RR Lyr observations (Kolenberg et al., 2011). For our models, the relations between the amplitudes of the signals at different frequencies are shown in Fig. 8. This is the analog of fig. 7 from Kolenberg et al. (2011), where the results for RR Lyr are presented. The model we discuss now has physical parameters and a mean period close to RR Lyr, but the modulation period is slightly longer (60 days). However, the described features of the frequency spectrum do not depend strongly on the modulation period and the analysis below is comparative.

Figure 8: Amplitude ratios, , , and amplitude of the half-integer frequencies, plotted versus the frequency (harmonic orders indicated by vertical dashed lines). Results for the model of set M6 with .

For we observe an exponential decrease with harmonic order . At the amplitude ratio drops to very small values (just like in RR Lyr). For the modulation components the decrease of amplitude ratios ( and ) is less steep. It is also less steep for than for – all in agreement with RR Lyr. For we observe a tail at high also in agreement with RR Lyr. For low , increases (above 1) in the case of RR Lyr. Here it is not the case, but such behaviour can be obtained for other models (see next section).

For the amplitudes of the HIF’s resulting from period doubling, in the discussed model the highest peak is located at followed by and . For RR Lyr the highest peak is observed at followed by and . We discuss this discrepancy in more detail in Section 3.5.

Frequency spectra analysis – model sequences

A detailed frequency analysis was conducted for all the model sequences discussed in this paper. Typical results for models with =6800 K and different parameters of turbulent convection modulation are shown in Fig. 9. These models are similar to that of set M6 discussed in the previous section, have the same amplitude of mixing-length modulation (50 per cent), but a longer modulation period (set M8) or a different shape of the modulation (triangular for MT1 and MT2). We stress that the results discussed below are characteristic for all our model sequences.

Figure 9: Relations between the amplitudes of the triplet components for the models of sets M6, M8, MT1 and MT2 (all with =6800 K). Panel (a): vs. in logarithmic scale; panel (b): vs. ; panel (c): vs. , and panel (d): vs. . The distribution of values from the MACHO database peaks at 0.3, which is marked with the horizontal dashed line for reference.

For the amplitudes of the harmonic frequencies, (panel a of Fig. 9), an interesting behaviour is observed at harmonic order around 7. The amplitudes are already very small at this order so the plots are in logarithmic scale to reveal the details. Local minima are clearly visible. They fall at , except for the model with a longer modulation period (set M8, ). In all these models period doubling is present, however the origin of the discussed minima and their possible connection with period doubling is not clear (see also Section 3.5).

Now we discuss the amplitudes of the modulation components, and (panels b and c of Fig. 9). For we observe more ‘erratic’ behaviour than for . At low orders, an increase of the amplitude ratio (above 1) is possible, just as it is observed in RR Lyr. decreases with increasing order. A plateau is visible for harmonic orders between 3 and 5.

Panel d of Fig. 9 shows the variation of the parameter defined as (Alcock et al., 2003),


with . The line at is marked for reference (the distribution of from the MACHO database peaks at this value, Alcock et al., 2003; Kolenberg et al., 2011). It is clear that in all our models is positive at low harmonic orders and hence, the amplitudes of the higher frequency triplet components are higher. In the intermediate range of harmonic orders the situation is reverted and at largest the higher frequency components are higher again.

We focus our attention on the triplet components around the main frequency . The amplitudes of these triplet components are the most robust outcome of our model and should be compared with observation first. In about 75 per cent per cent of Blazhko RR Lyrae stars, the higher frequency side peak has a larger amplitude than the lower frequency side peak and thus is positive (Alcock et al., 2003). This is in agreement with our models. However, in 25 per cent of the observed Blazhko variables is negative, as the lower frequency side peak has a larger amplitude. Unfortunately, in none of our models this is the case. In Fig. 10 we plot the values of vs. for the models of sets M7L4, M7, M7L6, M7L7, with different luminosities of the initial non-modulated models. The trend of the decreasing values with decreasing luminosity as well as with decreasing temperature is clearly visible. However, in no case is negative. This is true for all model sequences we have investigated, regardless of the properties of turbulent convection modulation (period, amplitude and shape). As in a quarter of known Blazhko RR Lyrae stars is negative, we regard the failure to reproduce negative values in our models as another significant challenge for the Stothers mechanism.

Figure 10: plotted vs. the for models with different luminosities of the initial non-modulated models (sets M7L4, M7, M7L6, M7L7).

3.5 Period doubling phenomenon

The period doubling phenomenon, which manifests in alternating shapes of the light curves of the consecutive pulsation cycles, was discovered very recently in Kepler RR Lyr data (Kolenberg et al., 2010a). The effect is also clearly visible in two other Kepler targets, V808 Cyg (KIC 4484128) and V355 Lyr (KIC 7505345), both being Blazhko variables (Szabó et al., 2010). In some other Blazhko stars, the confirmation of period doubling requires longer observation. The strength of period doubling depends on the phase in the Blazhko cycle. Period doubling was not found in any non-modulated RR Lyrae variable so far (Szabó et al., 2010).

The phenomenon of period doubling is not new. Period doubling is clearly present in RV Tauri stars which show alternating deep and shallow minima in their light and radial velocity curves. It was also found in radiative hydrodynamic models of W Vir stars (Buchler & Kovács, 1987; Kovács & Buchler, 1988), Cepheids and BL Her variables (Moskalik & Buchler, 1990, 1991; Buchler & Moskalik, 1992). Moskalik & Buchler (1990) traced the origin of the period doubling in hydrodynamic models to the destabilising role of the half integer resonance between the fundamental mode and a higher order overtone ( equal to 1 or 2). They showed that the period doubling can occur close to the resonance centre when the resonant overtone mode is either excited or only weakly damped. Moskalik & Buchler (1990) also showed that, depending on the amount of dissipation in the model, the half-integer resonance leads either to single period-doubling or to the period-doubling cascade (Feigenbaum cascade) and chaos.

The period doubling phenomenon in Blazhko RR Lyrae stars observed by Kepler was analysed in detail by Szabó et al. (2010), who also proposed the underlying mechanism. Their convective hydrodynamic models (with fixed convective parameters) display period doubling, the origin of which was traced to the 9:2 resonance between the fundamental mode, and the ninth order overtone, . In these models, the ninth overtone is a trapped envelope mode (see e.g., Buchler & Kolláth, 2001) with a much higher growth rate than that computed for the neighbouring overtone modes. As noted above, the weak damping of the ninth overtone favours the occurrence of period doubling. In their detailed analysis Kolláth et al. (2011) used the relaxation technique (Stellingwerf, 1974) to determine the stability of the fundamental mode pulsation through the Floquet stability coefficients. They showed that indeed the fundamental mode is destabilized through the resonance and excluded all other possible half-integer resonances. In addition, they showed that the resonance can lead not only to a single period doubling, but also to the period-four and period-eight limit cycles (the Feigenbaum cascade).

Our models also display period doubling, clearly visible in the light curves (Fig. 2, Figs. 35) and in the frequency spectra (Figs. 7 and 8). It occurs only in models of higher temperatures, with in a range 6700 K–6900 K (depending on the model sequence) and always at phases around maximum pulsation amplitude and minimum values of the mixing length (see right columns of Figs. 35). Our linear analysis confirms the findings of Szabó et al. (2010) and Kolláth et al. (2011). During the phases of minimum mixing length, in which period doubling occurs, the equilibrium models computed assuming corresponding, fixed values of the mixing length are very close to the resonance centre. We note that the mixing length has to be sufficiently small. Period doubling does not occur in models with small amplitude of the mixing-length modulation (Fig. 3). Also, the ninth overtone in our models is trapped in the outer envelope and is close to being unstable.

For the period doubling to occur, the model has to be close to the resonance condition. As discussed by Szabó et al. (2010), if the Stothers mechanism is indeed operational in Blazhko variables and the convective structure of the star varies during the Blazhko cycle (as it does in our models), it is natural that period doubling occurs only at certain Blazhko phases, at which the physical conditions are favourable. Qualitatively, this is what we see in our models. Favourable conditions (proximity to the 9:2 resonance) occur during the phases of low mixing length. These are the phases at which period doubling occurs. We note that in our models period doubling is strictly repetitive and always appears at the same phase of the Blazhko cycle. In fact, at other phases the period doubling does not vanish entirely. Its remnants serve as a seed for the fast growth of alternations during the next Blazhko cycle. We would like to stress that turbulent convection itself does not cause period doubling, which is a resonant phenomenon and may occur in purely radiative models as well. The modulation of turbulent convection in our models just changes the structure of the outer stellar layers, bringing the instantaneous model periods close to the resonance condition.

As discussed above, the transient occurrence of the period doubling during the Blazhko cycle can be naturally interpreted within the Stothers model. Nevertheless, a more quantitative comparison with observations reveals some discrepancies. In our models period doubling is limited to the phases of maximum pulsation amplitude. For the three Kepler Blazhko RR Lyrae stars it is most prominent in the phases preceding the maximum pulsation amplitudes, but not only, it is also visible close to the phases of minimum pulsation amplitude (Szabó et al., 2010). Certainly it is present in a much wider range of Blazhko phases than it is in our models.

As for the frequency spectra, the amplitudes of the half integer components in all three stars observed by Kepler follow the same pattern. The highest peak is observed at followed by and . This is not what we compute in our models, as already mentioned in Section 3.4.1. In Fig. 11 we show the amplitudes of HIF’s, for models with different patterns of turbulent convection modulation ( is the same in these models, equal to 6800 K). These are the same models as displayed in Fig. 9. The highest peak is always located at followed by and . For the latter two peaks the amplitudes are comparable. A local maximum is visible at for our models of sets M6, MT1 and MT2. Such a maximum is expected if the 9:2 resonance is indeed responsible for the period doubling. The fact that it is very weak and not present in all our models is a puzzle.

Figure 11: Amplitudes of half integer frequencies for models of sets M6, M8, MT1 and MT2 (initial non-modulated model was the same for all four models displayed, and its effective temperature was 6800 K).

4 Additional models

Figure 12: Light curve variation through the Blazhko cycle. vs. and vs for models in which different -parameters of the convective model are modulated.

The results presented so far reveal some disagreements between the computed models and the observations. The frequency analysis shows that in all our models is positive, so the higher frequency side peak around the fundamental mode frequency has a higher amplitude than the lower frequency side peak. In about 25 per cent of Blazhko variables is negative. Also, we cannot reproduce the details of the light curve variation. In particular the anticorrelation between amplitude ratio, , and amplitude, , observed in RR Lyr through the majority of the Blazhko cycle cannot be reproduced. Probably, the most important challenge for the Stothers mechanism is a huge range of variation of the mixing length. This large variation is required to reproduce the ranges of light curve variation in strongly modulated stars like RR Lyr, as well as the pulsation period changes we observe in Blazhko variables.

To get an idea whether and how these difficulties may be overcome, we have computed additional sequences of models in which we vary other -parameters of the turbulent convection model than the mixing length (see Section 2.1). So far, we modulated the mixing length only, which is the natural choice, as the mixing-length parameter controls the overall strength of convection. Now, as an exercise, we vary the other parameters entering the turbulent convection model. We modulate the strength of eddy-viscous dissipation (), the strength of the source function (), the strength of the convective heat flux () and the strength of the turbulent dissipation (, turbulent-cascade term). During the mixing-length modulation, all the above terms were modulated simultaneously, as these -parameters enter the model in pair with the mixing-length parameter (, , and ; see e.g., Smolec & Moskalik, 2008a). We note that there is no physical justification behind the modulation of any particular term listed above. Only a detailed 3D magnetohydrodynamical computation could clarify how a variable magnetic field could affect particular phenomena connected with the turbulent convection. Unfortunately, such computations and even appropriate models do not exist currently.

In all computed model sequences, the shape of modulation is sinusoidal, its period is 40 days and the amplitude of its modulation is 50 per cent. Particular -parameters vary around the values defined in Table 2.

In Fig. 12 we plot the vs. and vs. relations for models with three different initial temperatures (=6300 K, 6600 K and 6800 K), just like in Figs. 35. The cross in each panel corresponds to the initial non-modulated model. Thick solid trajectories correspond to the model of set M7 in which the mixing length is modulated (default in this paper). As expected, the ranges of variation of the Fourier parameters are smaller if we vary other -parameters of the model than the mixing length, as in each case only one term of the convective model is modulated. We focus on vs. relation (upper panels of Fig. 12). Observed trends are not always simple (particularly for models with =6300 K) and they depend on the temperature of the model. The increase of with decreasing , as we observe in RR Lyr, can be reproduced most easily when only the strength of the convective heat flux is modulated. Then, however, the range of variation of the Fourier parameters is small, despite a rather huge modulation of the convective heat flux. Therefore, it is difficult to propose a modulation to get a better model for RR Lyr. We have to vary more than one convective parameter to get a large range of variation of the Fourier parameters. For the model with =6600 K one could modulate only , and and keep the eddy-viscous dissipation fixed. On the other hand, Fig. 12 suggests that the same modulation adopted in the hotter model (=6800 K) would lead to increasing with increasing .

We have also investigated whether the negative values can be obtained through modulating particular terms in the convective model. In Fig. 13 we plot the vs. for the discussed models. Negative, albeit very close to zero, values of are obtained only for the hottest models for which either only convective heat flux is modulated (=6600 K and =6800 K) or only eddy-viscous dissipation is modulated (=6800 K). Modulation of other components of the convective model always leads to positive .

In all our models, the convective parameters vary around the values defined in Table 2. As noted in Section 2.3, with such values we can reproduce the observational properties of the non-Blazhko RR Lyrae variables and also classical Cepheids (with smaller eddy-viscous dissipation). With these parameters we neglect the effects of turbulent pressure and overshooting from the convective zone. The computation of models including these effects is very time-consuming. In order to check whether the inclusion of turbulent pressure and overshooting can change our results we have computed one additional model sequence in which we set and . We adopted the same modulation parameters as for set M6 (Table 3) i.e., sinusoidal modulation with a period equal to 60 d and an amplitude of per cent. The results are qualitatively the same as described in the previous sections. All models are characterised by positive values of . For the lower temperatures, the light curve variation is qualitatively the same as presented in, e.g., Fig. 3 (set M6). At higher temperatures, the range of variation of the Fourier parameters becomes smaller for models including turbulent pressure. Consequently, for such models, an even larger amplitude of turbulent convection modulation would be necessary to reproduce the strongly modulated Blazhko variables.

Figure 13: plotted vs. the for models, in which particular terms of the convective model were modulated.

5 Summary and Conclusions

In this paper we have investigated whether the mechanism proposed by Stothers (2006) is capable of reproducing the light curve variation and the properties of the frequency spectra we observe in RR Lyrae Blazhko variables. The mechanism proposed by Stothers is extremely complicated, as it assumes the coupling between a variable magnetic field and turbulent convection in high amplitude, strongly non-linear pulsators. To get a variable magnetic field a rotational/turbulent dynamo is postulated. There are no comprehensive models dealing with all these processes. Even the current models to deal with turbulent convection in non-linear pulsation are rather simple, 1D, one-equation formulas.

Recent observational progress poses serious problems for the two most popular models to explain the Blazhko effect, the Magnetic Oblique Rotator/Pulsator model and Non-radial Resonant Rotator/Pulsator model (see Introduction). The Stothers mechanism remains as a scenario that has not been confronted with concrete challenges from observations. The variable, tangled magnetic fields postulated by Stothers (2006) cannot be easily detected, making the idea hard to verify on purely observational grounds (Kolenberg & Bagnulo, 2009). Its stochastic nature makes it attractive in the light of irregularities commonly detected in the Blazhko cycles of many variables. However, this is only a general idea, not supported by any detailed calculations or modelling. To advance with the theory we proposed a simple model to check whether the variation of convective structure of the star can lead to the modulation we observe in Blazhko stars. The modulation of turbulent convection is introduced into our models in an ad hoc way, neglecting the dynamical coupling with the postulated magnetic field.

Comparison of our models with overall properties of Blazhko variables, as well as a detailed comparison with the strongly modulated prototype RR Lyr, observed by Kepler, reveals several discrepancies with the observations and challenges for the Stothers mechanism. In our opinion, the most important objection is the required strength of the turbulent convection modulation. In order to reproduce the ranges of the light curve variation observed in strongly modulated stars like RR Lyr, we have to modulate the mixing-length by up to 50 per cent on a time-scale of several tens of days. The same estimate results from the analysis of pulsation period changes (Section 3.3, see also Molnár & Kolláth, 2010). The physical reality of such strong modulation is, in our opinion, questionable, although definite claims require detailed magnetohydrodynamic modelling of the problem, which is beyond the scope of this paper. We note here that it is not possible to infer the mixing-length variation in the subatmospheric layers directly from observations. In a recent paper, Preston (2011) analysed the variation of FWHMs (full width at half maximum) of spectral lines in several RR Lyrae-type stars. This parameter is a measure of the atmospheric turbulence. Preston shows that maximum value of FWHM varies with the Blazhko cycle. However, it is not clear if changing intensity of the atmospheric turbulence relates in any way to changes of the mixing-length in the subatmospheric layers. A strong variation of the FWHM occurs also in non-modulated RR Lyrae stars (see Preston’s fig. 3), and, we may add, also in non-modulated hydrodynamical models with constant mixing-length (see, e.g., Benz & Stellingwerf, 1985). In our opinion, the variation of the maximum FWHM reflects the behaviour of the turbulence generated by the velocity gradients in the atmosphere, rather than the possible variation of the mixing length in the deeper layers. Preston (2011) shows that the maximum value of the FWHM is strongly correlated with the pulsation amplitude, which varies during the Blazhko cycle (his fig. 6). The larger the pulsation amplitude, the larger the maximum value of the FWHM. On the other hand, we note that the maximum of the FWHM occurs at pulsation phases close to the minimum radius. These are the phases at which the velocity gradients in the atmosphere are the strongest. Strong velocity gradients generate strong turbulence. Thus, the larger the pulsation amplitude, the stronger the atmospheric velocity gradients and, consequently, the stronger the turbulence and the higher the FWHM maximum.6

The light variation we compute resembles that observed in Blazhko variables, however, the details cannot be reproduced, at least for RR Lyr. We note that our results can be used in the future for comparison with other Blazhko variables for which satellite data will be soon available. This would clarify how severe the discrepancies between the model and the observed light variation are. As for the frequency spectra, we cannot reproduce the asymmetry of the modulation side peaks around the main pulsation frequency. In our models, the higher frequency side peak always has a higher amplitude, which is not the case in about a quarter of the Blazhko variables (Alcock et al., 2003).

The critical analysis of the Stothers idea by Kovács (2009) is also worth mentioning. He points out that the predictions by Stothers (2006, 2010) of expected period changes through the instability strip, that agree with values observed in some Blazhko stars, are not correct, as they result from a misinterpretation of the linear and non-linear model periods.

In the light of our results, the idea proposed by Stothers (2006) faces difficulties, and should be treated with caution. It needs refinement and further detailed studies in order to put it on solid physical grounds. Certainly, it is worth further investigation, but other possibilities to explain the Blazhko phenomenon should be explored as well.


Funding for this Discovery mission is provided by NASA’s Science Mission Directorate. The authors gratefully acknowledge the entire Kepler team, whose outstanding efforts have made these results possible.

Model computations presented in this paper were conducted on the psk computer cluster in the Copernicus Centre, Warsaw, Poland. We are grateful to James Nemec for comments on this manuscript. RS and KK are supported by the Austrian Science Fund (FWF projects AP 21205-N16 and T359/P19962, respectively).


  1. pagerange: Variable turbulent convection as the cause of the Blazhko effect – testing the Stothers modelVariable turbulent convection as the cause of the Blazhko effect – testing the Stothers model
  2. pubyear: 2010
  3. We note that our model predicts strictly periodic Blazhko cycles, which result from our assumption of a strictly periodic modulation of the mixing-length. This is not a necessary feature of the model. We could easily modulate the mixing-length in a quasi-periodic manner (as proposed in Stothers’ paper), obtaining a quasi-periodic modulation. We have chosen periodic modulation for simplicity.
  4. See also the animated gif available in the electronic version of the Journal. It shows the light curve variation through the Blazhko cycle, including the period doubling phenomenon, as well as the variation of the Fourier parameters.
  5. http://archive.stsci.edu/kepler/
  6. The question whether variations of the FWHM in the Blazhko stars are predominantly caused by changing velocity gradients in the atmosphere (as we expect, and which is the case for non-modulated stars) or whether possible modulation of the sub-atmospheric convection can also play a role, can be examined in more detail by numerical computations. This requires repeating the analysis of Benz & Stellingwerf (1985) for both the modulated models and for non-modulated models of different pulsation amplitudes. Such work however, is beyond the scope of this paper.


  1. Alcock C. et al., 2003, ApJ, 598, 597
  2. Alexander D.R., Ferguson J.W., 1994, ApJ, 437, 879
  3. Asplund M., Grevesse N., Sauval A.J., Allende Prieto C., Kiselman, D., 2004, A&A, 417, 751
  4. Baranowski R. et al., 2009, MNRAS 396, 2194
  5. Benkő J. et al., 2010, MNRAS, 409, 1585
  6. Benz W., Stellingwerf R.F., 1985, ApJ, 297, 686
  7. Buchler J.R., Goupil M.-J., 1984, ApJ, 279, 394
  8. Buchler J.R., Kolláth Z., 2001, ApJ, 555, 961
  9. Buchler J.R., Kolláth Z., 2011, ApJ submitted; arXiv:1101.1502
  10. Buchler J.R., Kovács G., 1987, ApJ, 320, L57
  11. Buchler J.R., Moskalik P., 1992, ApJ, 391, 736
  12. Chadid M., Wade G.A., Shorlin S.L.S., Landstreet, J.D., 2004, A&A, 413, 1087
  13. Chadid M. et al., 2010, A&A, 510, 39
  14. Detre L., Szeidl B., 1973, in Fernie J.D., ed., Variable Stars in Globular Clusters and Related Systems, p. 31
  15. Dziembowski W.A., Mizerski T., 2004, AcA, 54, 363
  16. Feuchtinger M., 1999, A&A, 351, 103
  17. Jenkins J. et al., 2010a, ApJ, 713, L87
  18. Jenkins J. et al., 2010b, ApJ, 713, L120
  19. Jurcsik J., Benkö J., Szeidl B., 2002, A&A, 390, 133
  20. Jurcsik J. et al., 2006, AJ, 132, 61
  21. Jurcsik J. et al., 2008, MNRAS, 391, 164
  22. Jurcsik J. et al., 2009, MNRAS, 400, 1006
  23. Jurcsik J., Szeidl B., Clement C., Hurta Zs., Lovas M., 2011, MNRAS, 411, 1763
  24. Koch D.G. et al., 2010, ApJ, 713, L79
  25. Kolenberg K., 2008, JPhCS, 118, 2060
  26. Kolenberg K., Bagnulo S., 2009, A&A, 498, 543
  27. Kolenberg K., Guggenberger E., 2007, CoAst, 150, 381
  28. Kolenberg K. et al., 2009, MNRAS, 396, 263
  29. Kolenberg K. et al., 2010a, ApJ, 713, L198
  30. Kolenberg K., Fossati L., Shulyak D., Pikall H., Barnes T.G., Kochukhov O., Tsymbal V., 2010b, A&A, 519, A64
  31. Kolenberg K. et al., 2011, MNRAS, 411, 878
  32. Kolláth Z., Buchler J.R., Szabó R., Csubry Z., 2002, A&A, 385, 932
  33. Kolláth Z., Molnar L., Szabó R., 2011, MNRAS in the press; arXiv:1102.0157
  34. Kovács G., 2009, AIPC, 1170, 261
  35. Kovács G., Buchler J.R., 1988, ApJ, 334, 971
  36. Kovács G., Kanbur S., 1998, MNRAS, 295, 834
  37. Kovács G., Buchler J.R., Davis C.G., 1987, ApJ, 319, 247
  38. Kuhfuß R., 1986, A&A, 160, 116
  39. LaCluyzé A. et al., 2004, AJ, 127, 1653
  40. Lenz P., Breger M., 2005, Communications in Asteroseismology, 146, 53
  41. Marconi M. 2009, AIP Conf. Ser., 1170, 223
  42. Marconi M., Clementini G., 2005, AJ, 129, 2257
  43. Molnár L., Kolláth Z., 2010, JPhCS, 218, 2027
  44. Moskalik P., Bucher J.R., 1990, ApJ, 355, 590
  45. Moskalik P., Bucher J.R., 1991, ApJ, 366, 300
  46. Nemec J. et al., 2011, MNRAS submitted
  47. Nowakowski R.M., Dziembowski, W.A., 2001, AcA, 51, 5
  48. Preston G.W., 2011, AJ, 141, 6
  49. Seaton M., 2005, MNRAS, 362, L1
  50. Shibahashi H., 2000, ASPC, 203, 299
  51. Simon N.R., Lee A.S., 1981, ApJ, 248, 291
  52. Smolec R., 2005, AcA, 55, 59
  53. Smolec R., 2009, PhD Thesis, N. Copernicus Astronomical Centre, Warsaw
  54. Smolec R., Moskalik P., 2008a, AcA, 58, 193
  55. Smolec R., Moskalik P., 2008b, AcA, 58, 233
  56. Stellingwerf R.F., 1974, ApJ, 192, 139
  57. Stothers R.B., 2006, ApJ, 652, 643
  58. Stothers R.B., 2010, PASP, 122, 536
  59. Szabó R. et al., 2010, MNRAS, 409, 1244
  60. Szczygiel D.M., Fabrycky D.C., 2007, MNRAS, 377, 1263
  61. Szeidl B., Jurcsik J., 2009, Communications in Asteroseismology, 160, 17
  62. Van Hoolst T., Dziembowski W.A., Kawaler S.D., 1998, MNRAS, 297, 536
  63. Wuchterl G., Feuchtinger M.U., 1998, A&A, 340, 419
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description