IS THE BLAZHKO EFFECT THE BEATING OF A NEAR-RESONANT DOUBLE-MODE PULSATION?

Variable stars of the RR Lyrae type often exhibit a slow modulation of their pulsations known as the Blazhko effect. First reported in 1907 by Russian astronomer Sergey Nikolaevich Blazhko (Blazhko 1907), the effect remains mysterious, with none of the numerous proposed explanations achieving widespread acceptance. See the recent reviews by Kovacs (2009) and Kolenberg (2012). Recent attempts to explain the Blazhko effect include the following.

• 1.  
  "The ninth overtone model" which involves an unstable resonance interaction between the fundamental mode and the ninth radial overtone, the oscillations of which lead to modulation of the fundamental (Szabo et al. 2010; Buchler & Kollath 2011; Kollath et al. 2011).
• 2.  
  "The Stothers model" in which an oscillatory magnetic field induces variations in the strength of turbulent convection which in turn modulates the fundamental Stothers (2006, 2010).
• 3.  
  "The shock model" in which a shock wave periodically generates turbulence in the atmosphere which in turn modulates the fundamental (Gillet 2013).

All of these models have significant problems which are discussed in Section 4.

In this paper a simple explanation is put forth for the Blazhko effect which involves (at least) two periodic oscillations, PO1 and PO2, of slightly different frequency. Typically one of these oscillations is strongly non-sinusoidal. The slow phase slip between these oscillations results in a combined waveform that displays a beat frequency or apparent modulation at the Blazhko frequency. We refer to this as "the beating-modes model." It is the only model where there is no actual modulation of the amplitude of the primary mode, rather the apparent modulation is an illusion produced by the partial cancellation that results from adding to it the other mode (or modes) which have nearly the same frequency but differing phases. This hypothesis is supported by a simple (non-hydrocode) model of the dynamics that accurately reproduces Kepler data for RR Lyr, as shown in Figure 1. In addition, the beating-modes model provides an explanation for the mysterious "double-maxima" waveform of V445 Lyr (See Section 3.1), providing additional strong evidence in its favor.

Zoom In Zoom Out Reset image size

Referring to Figure 1, note the many features of the data that are reproduced by the model including the shape of the pulsation waveform, the phase offset between the upper and lower modulation envelope functions, the size of the bump feature, and the behavior of the bump: hesitating near the bottom of the pulsation waveform, moving upwards for a while, fading out and disappearing for a while, and then reappearing near the bottom to complete the cycle. The disappearance of the bump lasts for about a third of the Blazhko cycle starting near the Blazhko minimum. The bump behavior of the model is a little surprising since it is entirely brought about by the addition of PO2, a pure sinusoid. PO1 by itself displays a bump at a fixed location about one quarter up from the bottom.

PO1 and PO2 may be generated by two unstable modes of the star. It is generally believed that the pulsation frequency for Blazhko stars corresponds to the fundamental mode, so this would be the likely choice for PO1, while PO2 may be generated by a nonradial mode of nearly the same frequency. An alternate possibility is that PO2 is a stable nonradial mode which becomes unstable by acquiring additional energy from the fundamental. (This is only possible if a mechanism exists that can accomplish this without phase-locking.) Since nonradial modes are involved, the usual 1D radial hydrocode software will not be able to test this model. Multi-D codes are now being developed which may be able to accomplish this at some future date (see e.g., Mundprecht et al. 2013; Geroux & Deupree 2014).

There have been several previous papers in support of the beating-modes model, but the current work goes much further in terms of finding a simple model that actually provides a remarkably good fit to the observed Blazhko data, including the mysterious motion and disappearance of the so-called bump feature. Guggenberger et al. (2008), building on earlier work by Breger & Kolenberg (2006), proposed a test of the phasing behavior in an effort to determine if beating was responsible for the Blazhko effect. They applied this to the two well-observed Blazhko RR Lyrae stars SS For and RR Lyr, with results consistent with the beating of two oscillations. Cox (1993), based on a combination of nonradial modal stability analysis combined with radial hydrocode analysis, proposed a near-resonant beating process between the fundamental and a nonradial mode, most likely the g₄, l = 1 mode. (Note that the stability analysis was for a static star and so might be significantly impacted by the presence of the high amplitude fundamental.) Sixteen years later, Cox (2009), working with improved code, finds that nonradial modes near the fundamental are always just slightly stable, and proposes that the g₄, l = 1 mode becomes unstable through nonlinear interaction with the fundamental. He still refers to this as a beating process not as a resonant mode interaction. Borkowski (1980) proposed a near-resonant double-mode pulsation involving the fundamental and either the second or third radial overtone, where the overtone frequency is higher than the second harmonic of the fundamental by the Blazhko frequency. Because the overtone is near the second harmonic instead of the fundamental itself, the Borkowski model is not validated by the results of this paper. Also since it only involves radial modes, one might have expected it to have been observed in 1D hydrocode models, but apparently it has not.

A variant of the beating-modes model is found in the work of Nowakowski & Dziembowski (2001). Their approach involves a resonance between the fundamental mode and a pair of nonradial modes with opposite orders ±m. Modes with nonzero m have oscillations as a function of the azimuthal angle. These can be thought of as traveling wave modes that move around the star either with or against the rotational motion of the star depending on the sign of m. As a result, these are split in frequency by the rotation. These modes can derive energy through nonlinear interaction with the fundamental if they satisfy the resonance condition that the sum of their frequencies must be exactly twice the fundamental frequency. As a result they generate a triplet spectrum with one peak on either side of the fundamental. According to Equation (77) of Nowakowski & Dziembowski (2001), the resulting Blazhko frequency is some fraction of the rotational frequency, most likely around 1/2. The pair are assumed initially stable and made unstable through the resonance. The resonant interaction between the modes is presumed stable and therefore the mode amplitudes and frequencies do not vary in time. Since all three peaks of the triplet correspond to actual modes, this is another type of beating-modes model, where the combined waveform appears modulated even though no actual modulation taking place. The main problem with this model is that there seems to be no good explanation for the amplitude disparity that is often observed for the opposite side peaks. However it seems reasonable to consider, due to the wide variation of behavior observed among Blazhko stars, that this model could apply to some stars which have a relatively balanced triplet.

The spectra for Blazhko stars includes progressions of side peaks to the left and right of the fundamental and each of its harmonics, all evenly spaced apart by the Blazhko frequency. These tend to diminish in amplitude when moving away to the left or right and often only the first side-peak is visible above the noise on each side. Most side peaks are nonlinear mixing product peaks and do not correspond to actual excited modes of the star. Usually these side peaks are attributed to modulation of the fundamental mode. However, if at least one of the side peaks immediately adjacent to the fundamental corresponds to an actual mode, this is sufficient through nonlinear processes (as discussed below) to generate the entire set of side peaks. The other frequencies are integer linear combinations of these two frequencies:

$$\begin{eqnarray}&&{{\omega }_{nm}}=n{{\omega }_{1}}+m{{\omega }_{2}},\end{eqnarray} \tag{ 1 }$$

where ω_nm is a side-peak frequency, n and m are integers (positive or negative), ω₁ is the frequency of PO1 (the fundamental), and ω₂ is the frequency of PO2. So, for example, the side-peak that is opposite to ω₂ (on the other side of the fundamental) has n = 2 and m = −1. If one side-peak is notably larger than the others then it is the likely choice for PO2. In the case of RR Lyr, there is a significantly larger peak on the higher frequency side. A spectral analysis (using the Period04 software package) of the Kepler data used in Figure 1(a) gives the amplitude (in units of 10⁷) of the upper peak as 0.0487 and of the lower peak as 0.0097, giving an amplitude ratio of about 5.02 and corresponding power ratio of 25.2. The central peak amplitude is about 0.1777.

It is well known that amplitude and frequency modulation can generate a side-peak spectrum similar to the one observed (see e.g., Szeidl & Jurcsik 2009; Benko et al. 2011; Szeidl et al. 2012), however a quasiperiodic dynamical system will also generate such a spectrum. A quasiperiodic system is one whose oscillations contain two base frequencies that are not rationally related. These would correspond to the frequencies of PO1 and PO2 and the apparent modulation frequency is the difference between them. The complete spectrum of such a system consists of the set of frequencies given by Equation (1). It can be obtained by a two-dimensional Fourier series expansion as is shown in Appendix B. Nonlinear processes can change the amplitudes of the (side-peak) frequencies including those which are initially zero. We take as a starting point a spectrum that includes PO1 (the fundamental mode frequency with its harmonics) and PO2 (a sine wave). There are two nonlinear mechanisms for generating the complete side-peak array. One is the process that translates the pulsational motion of the star into a variation in the light output. There is no reason to assume this process is entirely linear. So, for example, if the light output depends slightly on the square of the sum of the radial velocities of the two oscillations, then there will be cross terms generated by the square which will introduce components at the triplet (or first order) side peaks for the fundamental and all of its harmonics. Higher order nonlinear terms will generate higher order side peaks. A second mechanism is generated by nonlinear coupling between the active modes. Although such coupling is (by assumption) insufficient to produce phase-locking between the modes, it can introduce side peaks that were previously missing, change the amplitudes of peaks that were already present and shift the base frequencies. In the first mechanism, the side-peak frequencies are only present in the light curve, while in the second they are also physically present in the stellar vibration. Both mechanisms may be active at the same time.

For purposes of matching the Kepler data, PO1 and PO2 are approximated as constant in amplitude and frequency throughout the Blazhko cycle, which is equivalent to saying that the nonlinear interactions between them are neglected as relatively small. PO1 is a highly non-sinusoidal waveform whose frequency ω₁ is the pulsation frequency. Based on success at fitting the data with the model, it appears that this oscillation is made up of two sub-components which we will call PO1a and PO1b both of which are sawtooth-like waveforms. PO1a is larger in amplitude, while PO1b is shifted in phase and generates the bump feature. The most likely possibility is that PO1 is generated by a single mode with complicated dynamics. It is generally thought that the bump is an echo, either a reflection from the stellar core or the atmosphere (Guggenberger & Kolenberg 2006) and thus PO1b is the echo of PO1a. An alternate, but less likely, possibility is that PO1b is actually a separate mode, probably nonradial, that is phase-locked to the fundamental. This would make it similar to Cepheids which are thought to exhibit a bump due to a resonance with another mode (Simon 1977; Buchler et al. 1990; Bono et al. 2000). But regardless of the reason for the bump, the important point is that the model generates a reasonably good approximation to the observed waveform, enabling a test of the beating-modes hypothesis. The model output for PO1 is shown in Figure 2. Details of how it is generated are given below. PO2, on the other hand, is approximated as a pure sinusoid as this seems to do a good job of reproducing the observed Blazhko behavior. The frequency of PO2 is ω₂ = ω₁ + ω_B, where ω_B is the Blazhko frequency. In Section 5 it is discussed how the current work relates to the "hybrid mode model" developed previously by the author (Bryant 2014).

Zoom In Zoom Out Reset image size

For PO1a and PO1b we will use the same one-zone model developed in Bryant (2014), although here we are only using the model to generate the sawtooth-like waveform needed for PO1. The model details are given in Appendix A. This one-zone stellar model is similar to one developed by Baker (1966) and used extensively to study Cepheid and RR Lyrae variable stars (see, e.g., Stellingwerf 1972, 1986; Stellingwerf et al. 1987; Munteanu et al. 2005). When excited to high amplitude, this model produces a sawtooth-like velocity function, a waveform that is commonly seen in the light curves of RR Lyrae and Cepheid variable stars. For simplicity, we will use identical model parameters to generate PO1a and PO1b, the only differences are that PO1b has its amplitude reduced by a scale factor and it is shifted in phase relative to PO1a. Both sub-components are made equally non-sinusoidal, which seems to produce a very good fit to the observational data.

Let v_c be the combined velocity of PO1 and PO2. For PO1 we use v₁(t) from Equation (A4) in Appendix A and for PO2 we use a sine wave, with phase and scale factor chosen to fit the Kepler data:

$$\begin{eqnarray}&&{{v}_{c}}(t)={{v}_{1}}(t)+0.05{\rm sin} ({{\omega }_{2}}t-0.6).\end{eqnarray} \tag{ 2 }$$

The light flux output L(t), as plotted in Figure 1(b), is obtained from:

$$\begin{eqnarray}&&L(t)/{{10}^{7}}=1.25+1.3{{v}_{c}}(t).\end{eqnarray} \tag{ 3 }$$

The coefficients used were selected to obtain good visual agreement with the Kepler data in Figure 1(a).

We observe that for the Kepler data, the upper and lower Blazhko envelopes have differing amplitudes with the lower one being considerably smaller than the upper. The model on the other hand has upper and lower amplitudes of the envelope functions that are equal to each other and to the amplitude of PO2. We take a phenomenological approach to this problem and assume that the observed stellar flux is approximately equal to some nonlinear function of the combined radial velocity. In particular we add a correction term proportional to $v_{c}^{2}(t)$ to the expression for the light flux. The result is shown in Figure 1(c) and obtained from

$$\begin{eqnarray}&&{{L}_{{\rm corr}}}(t)/{{10}^{7}}=1.18+1.3{{v}_{c}}(t)+1.56v_{c}^{2}(t).\end{eqnarray} \tag{ 4 }$$

The effect of the correction term is to stretch the upper part of the graph increasing the amplitude of the upper envelope and compress the lower part decreasing the amplitude of the lower envelope. Note that there may be other ways to achieve the same effect, e.g., it could conceivably result from a nonlinear interaction between the modes that generate PO1 and PO2.

3.1. Multiple Frequency Blazhko

3.2. Different Triplet Symmetry

3.3. Long Period Blazhko

4.1. The Ninth Overtone Model

4.2. The Shock Model

4.3. The Stothers Model

In a previous paper (Bryant 2014) it was shown that an accurate fit to second quarter Kepler data for RR Lyr could be achieved by a hybrid mode consisting of two component modes of the same frequency. The first of these components is a highly non-sinusoidal sawtooth-like waveform, while the second is approximated as a sinusoid. Although locked in frequency, the two components are allowed to change independently as a function of the Blazhko phase. The first component varies little in amplitude or frequency throughout the Blazhko cycle, while the second changes quite strongly both in amplitude and in phase relative to the first component. By optimizing the fit to the Kepler data, the amplitude and phase behavior was obtained and displayed in Figure 4 of Bryant (2014).

The current work is related to that earlier work in the following way: It is found that the phase and amplitude dynamics of the second component of the earlier description can be approximately generated by adding together two constant amplitude sine waves, one with frequency ω₁ and a second of lower amplitude with frequency ω₁ + ω_B. This was noted previously by Guggenberger et al. (2008) who found this type of phase and amplitude behavior in several RR Lyrae stars and attributed it to this type of beating process. The first component corresponds to PO1a in the current model. The second sine wave corresponds to PO2. The first sine wave corresponds to PO1b, although in the current model this is made non-sinusoidal. Making this non-sinusoidal helps to correct one of the problems with that earlier model, in that the bump feature generated by that model was not as pronounced as in the Kepler data.

It has been shown that a simple beating-modes model of the Blazhko effect can provide a remarkably good fit to Kepler data for RR Lyr. The model consists of two periodic oscillations of slightly different frequency that are simply added together. This is in contrast to the much more complex picture of resonant mode models that depend on nonlinear coupling between modes leading to an oscillatory exchange of energy between them. The current model reproduces for the first time the apparent motion, disappearance and reappearance of the bump feature, the phase difference between the upper and lower modulation envelope functions and their relationship to the motion of the bump. In addition, it has been shown that the amplitude disparity between the upper and lower envelope functions can be reproduced if the light curve is assumed to be a nonlinear function of the radial velocity. Previous work by Cox (1993, 2009) suggests that the modes involved could be the fundamental and the g₄, l = 1 nonradial mode. The model is easily extended to cases exhibiting multiple modulation frequencies by simply associating additional large near-resonant peaks in the spectrum with active modes of the star, i.e., one additional mode is required for each additional modulation frequency. Application of this extended model was discussed for the case of V445 Lyr. The mysterious double-maxima waveform observed for that case is explained providing very strong evidence that the hypothesis of this paper is correct.

In addition, problems with other recent models of the Blazhko effect have been discussed in some detail, including the ninth overtone model, the shock model, and the Stothers model.

The author thanks Katrien Kolenberg for a very helpful discussion and comments about the model, the manuscript, and the spectral analysis. The author thanks Elisabeth Guggenberger for helpful comments on the manuscript.

The equations of motion are derived from a simplified analysis of pressure and gravity acting on the portion of a star above a certain radius R₀. We may take R₀ to be the effective bottom of the stellar envelope if we are considering the fundamental mode (as we are in this paper) or the outermost radial node for other modes such as the first radial overtone. For large amplitude oscillations, the motion is ballistic or gravity dominant for most of the cycle, with pressure becoming dominant only near the minimum position of the cycle. This leads to a "bouncing ball" type of oscillation whose derivative is the sawtooth-like waveform mentioned previously. Regardless of the accuracy of this simplified picture of the dynamics, the ability to reproduce the sawtooth-like waveform shows that the one-zone model seems to capture an essential aspect of the true dynamics.

The model takes the form of an ordinary differential equation for a harmonic oscillator, but with the usual linear forcing term replaced with a nonlinear function f(x) that, in a very simplified way, captures the essence of the gravity and pressure forces, i.e., $\dot{x}=v$ and $\dot{v}=f(x)$ where x and v are position and velocity. This model is conservative, excitation and damping of the oscillations are omitted, but could be added to the model if desired. Also currently assumed to be small and ignored are nonlinear interactions between modes that could lead to energy transfer or phase shift. In Bryant (2014) the following expression is derived for the forcing function:

$$\begin{eqnarray}\begin{array}{rcl} f(x) & = & A{{({{R}_{1}}+x)}^{2}}/{{({{({{R}_{1}}+x)}^{3}}-1)}^{\gamma }} \\ {} & {} & -B/{{({{R}_{1}}+x)}^{2}}, \\ \end{array}\end{eqnarray} \tag{ A1 }$$

where:

$$\begin{eqnarray}&&A={{\omega }^{2}}V_{1}^{\gamma }{{(3\gamma R_{1}^{4}V_{1}^{-1}-4{{R}_{1}})}^{-1}},\end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray}&&B={{\omega }^{2}}R_{1}^{4}{{(3\gamma R_{1}^{4}V_{1}^{-1}-4{{R}_{1}})}^{-1}},\end{eqnarray} \tag{ A3 }$$

γ = 1 + 1/n_p, ${{V}_{1}}=R_{1}^{3}-1$, R₁ is a characteristic radius for the stellar material that lies above r = R₀ and n_p is the polytropic index. The problem has been rescaled so that R₀ = 1. To keep the values near unity, time is in units of 0.1 days. The values used in the model are ω = 1.5 rad/0.1 day, n_p = 1.5, R₁ = 1.1, ω₁ = 1.108328, and ω_B = 0.015393. The Blazhko period is 72 pulsation periods. Initial condition for PO1a: ${{x}_{1a}}(0)=-0.0691$, v_1a(0) = 0. Initial condition for PO1b is chosen so that it produces the identical oscillation except that it is leading in phase by 2.1058232 radians or 0.19 days: x_1b(0) = 0.185258, v_1b(0) = 0.078770. Note that ω is the frequency for small amplitude oscillations; the actual oscillation frequency drops as the amplitude increases and is equal to ω₁ for the specified initial conditions. The velocity for PO1, v₁(t), is the combined velocity of PO1a and PO1b, with scale factors chosen to fit the Kepler data:

$$\begin{eqnarray}&&{{v}_{1}}(t)={{v}_{1a}}(t)+0.15{{v}_{1b}}(t).\end{eqnarray} \tag{ A4 }$$

This function is plotted in Figure 2.

When a nonlinear dynamical system is quasiperiodic, it will have two base frequencies which are incommensurate, i.e., they do not have a rational ratio. In the phase space the dynamics will lie on a two-torus or doughnut surface. The dynamics will wind its way around this surface, taking an infinite amount of time to completely cover it. Motion on a two-torus can be described by two angle-like variables, θ₁ and θ₂, corresponding to the two ways in which one can go around the torus. We can define the angles so that they satisfy ${{\theta }_{1}}=({{\omega }_{1}}t){\rm mod} (2\pi )$ and ${{\theta }_{2}}=({{\omega }_{2}}t){\rm mod} (2\pi )$. The state of the system can be expressed as a function of these two angles:

$$\begin{eqnarray}&&{\boldsymbol{X}} ={\boldsymbol{X}} ({{\theta }_{1}},{{\theta }_{2}})\;=\;{\rm state\ vector\ of\ the\ system}.\end{eqnarray} \tag{ B1 }$$

The state is $2\pi$ periodic in both angles (since a $2\pi$ shift in either angle corresponds to one pass around the torus). Thus it can be expanded as a two-dimensional Fourier series:

$$\begin{eqnarray}&&{\boldsymbol{X}} =\mathop{\sum }\limits_{m=-\infty }^{\infty }\ \mathop{\sum }\limits_{n=-\infty }^{\infty }{{{\boldsymbol{A}} }_{mn}}{\rm exp} (im{{\theta }_{1}}+in{{\theta }_{2}})\end{eqnarray} \tag{ B2 }$$

where

$$\begin{eqnarray}\begin{array}{rcl} {{{\boldsymbol{A}} }_{mn}} & = & \frac{1}{{{(2\pi )}^{2}}}\int _{0}^{2\pi }d{{\theta }_{1}}\int _{0}^{2\pi }d{{\theta }_{2}}{\boldsymbol{X}} ({{\theta }_{1}},{{\theta }_{2}}) \\ {} & {} & \times {\rm exp} (-im{{\theta }_{1}}-in{{\theta }_{2}}). \\ \end{array}\end{eqnarray} \tag{ B3 }$$

In terms of time this becomes:

$$\begin{eqnarray}&&{\boldsymbol{X}} =\mathop{\sum }\limits_{m=-\infty }^{\infty }\ \mathop{\sum }\limits_{n=-\infty }^{\infty }{{{\boldsymbol{A}} }_{mn}}{\rm exp} (i{{\omega }_{nm}}t)\end{eqnarray} \tag{ B4 }$$

where the frequencies in the spectrum are

$$\begin{eqnarray}&&{{\omega }_{mn}}=m{{\omega }_{1}}+n{{\omega }_{2}}.\end{eqnarray} \tag{ B5 }$$

The amplitudes ${{{\boldsymbol{A}} }_{mn}}$ can also be obtained as an integral over time:

$$\begin{eqnarray}&&{{{\boldsymbol{A}} }_{mn}}=\mathop{{\rm lim} }\limits_{T\to \infty }\frac{1}{T}\int _{0}^{T}dt{\boldsymbol{X}} (t){\rm exp} (-i{{\omega }_{mn}}t).\end{eqnarray} \tag{ B6 }$$

The amplitude of these components must fall off as m and n become large.

Cases of three frequency quasiperiodicity (or higher) require a simple extension of the above analysis to a higher dimensional Fourier series.