© ESO 2018

1. Introduction

Precise timing measurements in close binaries routinely reveal variations in the eclipse timings. While for binaries with very short periods of less than three hours gravitational wave emission is the dominating means of angular momentum loss, magnetic braking can account for such effects in binaries with longer periods (see, e.g., Kraft et al. 1962; Faulkner 1971; Verbunt & Zwaan 1981; Parsons et al. 2013). A subclass of close binaries, mainly RS Canum Venaticorum (RS CVn) and post-common envelope binary (PCEB) systems, feature cyclic or nearly-periodic orbital period variations on timescales of a few years to decades that are incompatible with classical gravitational wave or magnetic braking models (see, e.g., Lanza et al. 1998; Brinkworth et al. 2006; Zorotovic & Schreiber 2013). For some systems, planetary companions may explain these variations (see, e.g., Qian et al. 2011; Beuermann et al. 2013; Nasiroglu et al. 2017; Han et al. 2017), while some of the planetary solutions were found to be dynamically unstable (Horner et al. 2013) or have been proven incorrect observationally (Hardy et al. 2015). On the other hand, the period variations in other systems such as QS Vir are still not well understood (Parsons et al. 2010).

Based on and inspired by earlier work by Matese & Whitmire (1983) and Applegate & Patterson (1987), Applegate (1992) proposed a new mechanism to explain cyclic orbital period variations in close binaries. The author assumed a time-dependent gravitational quadrupole moment modulated by the stellar activity cycle. Given a constant orbital angular momentum, an increasing quadrupole moment results in a stronger gravitational field and finally in a decreasing orbital radius and increasing orbital velocity, which can be observed as a reduction of the binary period.

To calculate the required energy as well as the expected quadrupole changes, Applegate (1992) considered a thin shell rotating in a point mass potential representing the rest of the star. The author assumed that the torque necessary to perform the angular momentum exchange is provided by subsurface magnetic fields of a few kG. Given such a field, luminosity variations on the 10% level and angular velocity variations on the 1% level are expected, leading to binary period variations compatible with the observed values in RS CVn systems (ΔP/P  ∼  10⁻⁵). Applegate (1992) made a number of testable predictions, which included significant period luminosity changes due to temperature variations and a 1:1 relation between the modulation period and the stellar activity cycle. Both processes have been observed in systems such as CG Cygni, accompanied by orbital period variations (Hall 1991), while luminosity variations can be observed by studying photospheric temperatures (see, e.g., Gray & Baliunas 1994).

Lanza et al. (1998) generalized the magnetohydrodynamic aspects of the pioneering work of Applegate (1992), included the effects of internal magnetic fields, and elaborated on the connection between the Applegate mechanism and different types of dynamo models, emphasizing that a careful study of the Applegate mechanism may allow distinguishing between them. The improved Lanza et al. (1998) model is based on the momentum balance equation and assumes a uniform ratio of plasma pressure to magnetic pressure over the star. Expansion of the gravitational potential in spherical harmonics allowed Lanza et al. (1998) to calculate the quadrupole potential at the surface of the star, finding that the angular velocity changes required to produce a certain amount of period variation are reduced by a factor of 2 compared to the original Applegate (1992) approach.

Furthermore, Lanza et al. (1998) found that the driver of the period variations is the magnetic dynamo of the active component, which can effectively transform kinetic energy from nuclear reactions into magnetic energy and is more likely to be an α²Ω dynamo, introducing a further observational aspect. Lanza & Rodonò (1999) stressed that the original Applegate (1992) ansatz may systematically overestimate the energies required to power the quadrupole moment variations as the exchanges between kinetic and magnetic energy are cyclic and partially reversible processes. Therefore, the authors concluded that the original model only provides upper limits for the required energy.

Lanza (2005) revisited the original work by Applegate (1992) and introduced an entirely new approach. Assuming that the angular velocity of the active component is only a function of the distance from its rotation axis, Lanza (2005) considered angular momentum redistributions not only between two layers, but rather general radial transport and redistribution modes within the convection zone. Imposing a strictly adiabatic convection zone and neglecting density perturbations, the author considered the equations of mass continuity and angular momentum conservation within a mean-field framework. Solving the angular momentum equation, Lanza (2005) calculated individual angular momentum redistribution modes from which he derived both the amplitude of the quadrupole moment change and the kinetic energy dissipated during the redistribution processes. Contrary to prior works, Lanza (2005) found that the power required to drive the observed levels of binary period variation in HR 1099 exceeds the luminosity provided by the active component by one or two orders of magnitude. Lanza (2006) extended the Lanza (2005) approach to angular velocity distributions that are functions of both radius and latitude, finding that the Applegate effect is still not a viable option in the case of HR 1099 and RS CVn systems in general.

Another route has been taken by Brinkworth et al. (2006), who generalized the original Applegate (1992) thin-shell ansatz to angular momentum exchanges between a finite core and a finite shell, including the core back-reaction to oblateness changes of the surrounding shell. Völschow et al. (2016) extended the Brinkworth et al. (2006) model to realistic stellar density profiles, derived analytic expressions to estimate the energy required to produce a certain level of period variation in a given system, and applied their full model to a set of 16 close binary systems, including 11 PCEBs. In line with previous works, the authors found that the Applegate effect cannot uniquely explain period variations in close binaries, and only 4 of the 16 systems were identified as potential Applegate candidates. Most recently, Navarrete et al. (2018) employed the analytic two-zone model by Völschow et al. (2016) to investigate how rotation affects the energetic feasibility of the Applegate effect and on which scale the activity cycle matches the observed modulation period, noting that the systems with the highest rotation rates are the most likely Applegate candidates. However, one critical simplification of both the Brinkworth et al. (2006) and Völschow et al. (2016) models is that the authors set the core-shell transition where it minimizes the energy necessary to power the Applegate effect.

In this paper, we extend the Lanza (2006) model by assuming a time-dependent magnetic field, velocity field fluctuations, and magnetic field fluctuations in the convection zone from which we can explicitly calculate the temporal evolution of the dissipated energy as well as the binary period variation. In addition, we consider superpositions of angular momentum redistribution modes instead of calculating the quadrupole moment changes caused by individual modes, as previously done by Lanza (2005, 2006). For the sake of consistency, we first apply our model to HR 1099 and use the system as an illustrative example for the basic predictions we make, before we extend our analysis to PCEB systems and identify the most likely Applegate mechanism candidates.

2. Model description

The formalism we employed is based on Lanza (2006). Sections 2.1 and 2.2.1 closely follow the author’s description, but our presentation puts an emphasis on a more algorithmic description. Starting from Sect. 2.3.1, we extend on the crucial aspects of the calculation by making an ansatz for the stellar parameter fluctuations, explicitly solving the temporal part of the underlying differential equation (see Sect. 2.2.1), elaborating on the evolution of the mechanism, and presenting an approach for describing superpositions of angular momentum redistribution modes.

2.1. Setup

Lanza (2006) assumed an inertial reference frame originating at the barycenter of the active component star, with the z-axis pointing in the direction of the stellar rotation axis. In line with Lanza (2006), we employ a spherical coordinate system where r is the distance from the origin, the colatitude θ is measured from the north pole and the azimuthal angle is ϕ, and we impose that all variables are independent of ϕ. The hydrodynamics of the turbulent convection zone are described employing a mean-field approach $\mathit{V}=\overline{\mathit{\upsilon }}+\mathit{\upsilon }\prime$ with mean velocity $\overline{\mathit{\upsilon }}$ and mean value fluctuation v′. Furthermore, we neglect the effect of meridional circulations inside the star by assuming that the mean velocity field arises purely from stellar rotation (Lanza 2006). The equation for the angular velocity ω = v_ϕ/(r sin θ) reads (see, e.g., Lanza 2005, 2006)

$$\begin{array}{cc}& \frac{\partial \omega }{\partial t}-\frac{1}{\rho r^4}\frac{\partial }{\partial r}\left(r^4{\eta }_{\mathrm{t}}\frac{\partial \omega }{\partial r}\right)-\frac{{\eta }_{\mathrm{t}}}{\rho r^2}\frac{1}{(1-{\mu }^2)}\frac{\partial }{\partial \mu }\left({(1-{\mu }^2)}^2\frac{\partial \omega }{\partial \mu }\right)\hfill \\ \hfill & =S(r,\mu ,t).\hfill \end{array}$$(1)

Here, η_t = η_t(r) is the turbulent dynamical viscosity, μ = cosθ. S(r, μ, t) is a source term that controls the temporal evolution of the angular momentum redistribution and is specified in Sect. 2.3.1. Equation (1) is solved assuming a stress-free boundary condition

$$\begin{array}{c}\hfill {\left(\frac{\partial \omega }{\partial r}\right)}_{r_{\mathrm{b}},R}=0,\end{array}$$(2)

r_b denotes the base of the stellar convection zone, and R is the radius of the star (Lanza 2006).

2.2. Angular momentum equation

2.2.1. Model framework

Introducing the state of rigid rotation Ω₀, the angular velocity of the star can be written as

$$\begin{array}{c}\hfill \mathrm{\Omega }(r,\mu ,t)={\mathrm{\Omega }}_0+\omega (r,\mu ,t),\end{array}$$(3)

where ω(r, μ, t) describes the deviation from rigid rotation (Lanza 2006). According to Lanza (2006), solutions of the angular velocity Eq. (1) have the form

$$\omega (r,\mu ,t)={\displaystyle \sum _{n=0}^{\infty }{\alpha }_n}(t){\zeta }_n(r)P_n^{(1,1)}(\mu )$$(4)

where $P_n^{(1,1)}$ are Jacobian polynomials and the index n is referred to as angular order. Furthermore, following Lanza (2006), the source term S (r, μ, t) can be developed into

$$S(r,\mu ,t)={\displaystyle \sum _{n=0}^{\infty }{\beta }_n}(t){\zeta }_n(r)P_n^{(1,1)}(\mu ).$$(5)

Based on the angular momentum conservation equation and using the properties of the Jacobian polynomials, Lanza (2006) derived a partial differential equation for the α_n, β_n, and ζ_n functions that can be separated into two equations. They can be solved independently from each other, namely

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\alpha }_n(t)}{\mathrm{d}t}-{\beta }_n(t)+{\alpha }_n(t){\lambda }_n=0,\end{array}$$(6)

which we refer to as temporal equation, as well as

$$\begin{array}{c}\hfill \frac{1}{\rho (r)r^4}\frac{\mathrm{d}}{\mathrm{d}r}(r^4{\eta }_{\mathrm{t}}(r){\zeta }_n^{\prime }(r))-n(n+3)\frac{{\eta }_{\mathrm{t}}(r)}{\rho (r)r^2}{\zeta }_n(r)+{\lambda }_n{\zeta }_n(r)=0,\end{array}$$(7)

which we refer to as the radial equation, and λ_n are eigenvalues of Eq. (7). The radial Eq. (7) together with the boundary conditions (2) define a regular Sturm–Liouville problem in the interval [r_b, R] for r_b >  0 (Lanza 2006).

2.2.2. Properties of the radial equation

For any given n, an infinite number of eigenvalues exist. Therefore, Lanza (2006) introduced a new index k (radial order) and denominated as λ_nk the kth eigenvalue of the nth radial equation. The first eigenvalue λ₀₀ is zero, the associated eigenfunction ζ₀₀ vanishes, and all eigenvalues are positive (see Lanza 2006).

Accounting for the radial order k, Eq. (4) can be recast into

$$\omega (r,\mu ,t)={\displaystyle \sum _k{\displaystyle \sum _{n=0}^{\infty }{\alpha }_{nk}}}(t){\zeta }_{nk}(r)P_n^{(1,1)}(\mu ),$$(8)

and the radial Eq. (7) now reads

$$\begin{array}{c}\hfill \frac{1}{\rho r^4}\frac{\mathrm{d}}{\mathrm{d}r}(r^4{\eta }_{\mathrm{t}}\zeta {{}^{\prime }}_{\mathit{nk}})-n(n+3)\frac{{\eta }_{\mathrm{t}}}{\rho r^2}{\zeta }_{\mathit{nk}}+{\lambda }_{\mathit{nk}}{\zeta }_{\mathit{nk}}=0.\end{array}$$(9)

Imposing |ω| ≪ Ω₀ (cf. Eq. (3)) and employing a linear approximation for the variation of the gravitational quadrupole moment (see Sect. 2.5), we only have to consider the cases n = 0 and n = 2 (Lanza 2006). Strictly speaking, in the case of a convective zone that spans the entire star (r_b = 0), the regular Sturm–Liouville problem turns into a singular Sturm–Liouville problem. The boundary condition at the lower end of the domain must be replaced by a regularity condition, and a different treatment of Eq. (9) is required to deal with singularities at r = 0. However, the stellar structure models we employ (see Sect. 3) all start at some radial coordinate r >  0, which ensures r_b >  0 even for fully convective stars.

2.2.3. Solving the radial equation

The most essential input to solve the radial equations is a stellar structure model that provides the mass distribution within the star M(r), the luminosity L(r), the density ρ(r), the temperature T(r), and the mean molecular weight m_μ(r). Based on this input, we can derive the gravitational acceleration

$$\begin{array}{c}\hfill g(r)=\frac{GM(r)}{r^2},\end{array}$$(10)

and the convective velocity employing mixing length theory u_c

$$\begin{array}{c}\hfill u_{\mathrm{c}}(r)={\left(\frac{{\alpha }_{\mathrm{ml}}L(r)}{40\pi r^2\rho (r)}\right)}^{1/3},\end{array}$$(11)

where we use a mixing length parameter α_ml = 1.5 (Lanza 2005). The pressure scale height h_p is given by

$$\begin{array}{c}\hfill h_{\mathrm{p}}(r)=\frac{k_{\mathrm{B}}T(r)}{m_{\mathrm{H}}{m}_{\mu }(r)g(r)},\end{array}$$(12)

with Boltzmann’s constant k_B and the hydrogen mass m_H (Lanza 2005). Finally, we can calculate the turbulent viscosity from

$$\begin{array}{c}\hfill {\eta }_{\mathrm{t}}(r)=\frac{1}{3}\rho (r){\alpha }_{\mathrm{ml}}{u}_{\mathrm{c}}(r)h_{\mathrm{p}}(r).\end{array}$$(13)

The stress-free boundary condition specified by Eq. (2) implies ${\zeta }_{nk}^{\prime }(r)=0$ at r = r_b and r = R as boundary conditions for ζ_nk. We assume that ζ_nk vanishes outside [r_b, R], that is, angular velocity variations only occur in the convection zone. Intermediate stellar structure values are calculated by a linear interpolation scheme. Solutions for ζ_0k and ζ_2k are computed with a shooting method (see, e.g., Fehlberg 1987; Press et al. 1992). While the boundary conditions imply ${\zeta }_{nk}^{\prime }(r_{\text{b}})=0$, the eigenfunction normalization at the bottom of the convective zone ζ_nk(r_b) can be regarded as a free parameter of the model only restricted by ζ_nk(r_b) ≪ Ω₀. The effect of this parameter is addressed in Sect. 3.

Lanza (2006) adopted ζ_nk(r_b) = 0.01Ω₀ throughout his calculations, that is to say, he normalized individual angular momentum redistribution modes to some level. In contrast, we solve the full temporal equation to find solutions for the α_nk(t) and β_nk(t) functions, which allows us to calculate the evolution of the redistribution processes and consider superpositions of the elemental redistribution functions ζ_nk(r).

2.3. Temporal equation

2.3.1. Framework

Following Lanza (2006), solutions for the temporal Eq. (6) are

$${\alpha }_{nk}(t)=\exp (-{\lambda }_{nk}t){\displaystyle \underset{0}{\overset{t}{\int }}{\beta }_{nk}}(t\prime )\exp ({\lambda }_{nk}t\prime )\text{d}t\prime +{\alpha }_{nk}(0),$$(14)

where β_nk is calculated from

$${\beta }_{nk}=E_{nk}{\displaystyle \underset{r_{\text{b}}}{\overset{R}{\int }}{\displaystyle \underset{-1}{\overset{1}{\int }}\rho }}{r}^{4}S(r,\mu ,t){\zeta }_{nk}{P}_n^{(1,1)}(1-{\mu }^2)\text{d}\mu \text{d}r.$$(15)

E_nk is a normalization constant¹ given by

$$E_{nk}^{-1}={\displaystyle \underset{r_{\text{b}}}{\overset{R}{\int }}\rho }{r}^4{\zeta }_{nk}^2\text{d}r{\displaystyle \underset{-1}{\overset{1}{\int }}{P}_n^{(1,1)}}{P}_n^{(1,1)}(1-{\mu }^2)\text{d}\mu .$$(16)

The integral involving the Jacobian polynomials simplifies to

$${\displaystyle \underset{-1}{\overset{1}{\int }}{P}_n^{(1,1)}}{P}_n^{(1,1)}(1-{\mu }^2)\text{d}\mu =\frac{8(n+1)}{(2n+3)(n+2)}.$$(17)

In a linear approximation for the quadrupole moment variation, we only have to consider the angular orders n = 0 and n = 2 (see Sect. 2.5) for which the term on the right-hand side can be evaluated as 4/3 for n = 0 and 6/7 for n = 2.

The source term S in Eq. (1) describes angular momentum transfer by Reynolds stresses and magnetic torques and takes the form

$$\begin{array}{c}\hfill S(r,\mu ,t)=-\frac{\mathrm{div}\mathit{\tau }}{\rho r^2(1-{\mu }^2)}.\end{array}$$(18)

The vector τ_i has the components

$$\begin{array}{c}\hfill {\mathit{\tau }}_i=rsin\theta ({\mathrm{\Lambda }}_{i\varphi }+\frac{1}{\tilde{\mu }}(B_i{B}_{\varphi }+M_{i\varphi })).\end{array}$$(19) ${\Lambda }_{i\varphi }=\rho \overline{{\upsilon }_i^{\prime }{\upsilon }_{\varphi }^{\prime }}$ is the Reynolds stress tensor, $\tilde{\mu }$ is the magnetic permeability, B is the mean magnetic field, and $M_{i\varphi }=\overline{B_i^{\prime }{B}_{\varphi }^{\prime }}$ is the Maxwell stress tensor (Lanza 2005). Finally, we arrive at the expression

$${\beta }_{nk}=-E_{nk}{\displaystyle \underset{-1}{\overset{1}{\int }}{\displaystyle \underset{r_b}{\overset{R}{\int }}{r}^2}}\text{div}\tau {\zeta }_{nk}{P}_n^{(1,1)}\text{d}\mu \text{d}r$$(20)

to calculate the β_nk functions from which we can then calculate the α_nk functions necessary to describe the temporal evolution of the star’s angular momentum distribution.

2.3.2. Stellar parameter fluctuations

As the Applegate mechanism is linked to and triggered by magnetic activity, we impose a phase factor of the form

$$\begin{array}{c}\hfill f(t)=sin\left({\omega }_{\mathrm{act}}t\right)\end{array}$$(21)

with ω_act = 2π/P_act to describe the fluctuations of all involved quantities where P_act is the activity cycle period (see, e.g., Rüdiger et al. 2002). For the magnetic field, we assume a simple azimuthal structure with B_r = 0, B_θ = 0 and

$$\begin{array}{c}\hfill B_{\varphi }(t)=B_{\mathrm{surf}}sin\left({\omega }_{\mathrm{act}}t\right)\end{array}$$(22)

accompanied by magnetic field fluctuations of the form

$$\begin{array}{c}\hfill B{{}^{\prime }}_i(r,t)=A_{\mathrm{B}}{B}_{\mathrm{surf}}sin\left({\omega }_{\mathrm{act}}t\right)\end{array}$$(23)

and velocity field fluctuations in the convection zone

$$\begin{array}{c}\hfill {\upsilon }_i^{\prime }(r,t)=A_{\mathrm{\upsilon }}{u}_{\mathrm{c}}(r)sin\left({\omega }_{\mathrm{act}}t\right),\end{array}$$(24)

where A_B and A_v are amplification coefficients and B_surf is the surface magnetic field amplitude. Assuming that these equations describe proper long-term ensemble means and imposing that both $B_i^{{}^{\prime }}{B}_j^{{}^{\prime }}$ and ${\upsilon }_i^{\prime }{\upsilon }_j^{\prime }$ are independent statistics for i ≠ j, we can make an ansatz for the Reynolds tensor

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{i\varphi }(r,t)=\rho (r)A_{\mathrm{\upsilon }}^2{u}_{\mathrm{c}}^2(r){sin}^2({\omega }_{\mathrm{act}}t)\end{array}$$(25)

as well as the Maxwell stress tensor:

$$\begin{array}{c}\hfill M_{i\varphi }(r,t)=A_{\mathrm{B}}^2{B}_{\mathrm{surf}}^2{sin}^2({\omega }_{\mathrm{act}}t).\end{array}$$(26)

The effect of a constant phase lag Δϕ between the phase factor of different coordinates is addressed in Sect. 3.5. We note that the stellar parameter fluctuation description we adopt for our model does not come from dynamo theory or hydrodynamic calculation, but is rather a simplified ad hoc approach to generate cyclic Maxwell and Reynolds stresses.

2.3.3. Solution of the temporal equation

We can now evaluate the source term S(r, μ, t). For a vector field F(r, θ, ϕ) in spherical coordinates, the divergence takes the form

$$\begin{array}{c}\hfill \mathrm{div}\mathit{F}={(\frac{1}{r^2}\frac{\partial }{\partial r}(r^2{\mathit{F}}_r),\frac{1}{rsin\theta }\frac{\partial }{\partial \theta }(sin\theta {\mathit{F}}_{\theta }),\frac{1}{rsin\theta }\frac{\partial }{\partial \varphi }({\mathit{F}}_{\varphi }))}^{\mathrm{T}}.\end{array}$$(27)

For our given magnetic field configuration and fluctuation setup, the divergence of the vector τ can be calculated as

$$\begin{array}{c}\hfill \mathrm{div}\mathit{\tau }=\sqrt{1-{\mu }^2}\frac{1}{r^2}\frac{\partial }{\partial r}\left[r^3{\tilde{B}}_{r\varphi }(r,t)\right]+2\mu {\tilde{B}}_{\theta \varphi },\end{array}$$(28)

where we defined

$$\begin{array}{c}\hfill {\tilde{B}}_{i\varphi }(r,t)={\mathrm{\Lambda }}_{i\varphi }(r,t)+\frac{B_i(r,t)B_{\varphi }(r,t)}{\tilde{\mu }}+\frac{M_{i\varphi }(r,t)}{\tilde{\mu }}.\end{array}$$(29)

Using this, the function β_nk(t) can be evaluated as

$$\begin{array}{c}\hfill {\beta }_{\mathit{nk}}=-E_{\mathit{nk}}{\displaystyle \underset{r_{\mathrm{b}}}{\overset{R}{\int }}{\zeta }_{\mathit{nk}}(r)\frac{\partial }{\partial r}\left[r^3{\tilde{B}}_{r\varphi }(r,t)\right]\mathrm{d}r}{\displaystyle \underset{-1}{\overset{1}{\int }}{P}_n^{(1,1)}\sqrt{1-{\mu }^2}}\mathrm{d}\mu .\end{array}$$(30)

The second term from div τ vanishes because of the symmetry of the μ integral. The remaining Jacobian polynomial integral can be evaluated analytically. First, we substitute μ = sinu and use the identity sin²u + cos²u = 1, which yields

$$\begin{array}{cc}& {\displaystyle \underset{-1}{\overset{1}{\int }}{P}_n^{(1,1)}\sqrt{1-{\mu }^2}}\mathrm{d}\mu \hfill \\ \hfill & ={\displaystyle \underset{-\pi /2}{\overset{\pi /2}{\int }}{P}_n^{(1,1)}(sinu)}\mathrm{d}u-{\displaystyle \underset{-\pi /2}{\overset{\pi /2}{\int }}{P}_n^{(1,1)}(sinu){sin}^2(u)}\mathrm{d}u.\hfill \end{array}$$(31)

Next, $P_n^{(1,1)}(\mu )$ can be written as

$$\begin{array}{c}\hfill P_n^{(1,1)}(\mu )=\frac{1}{n+2}{\displaystyle \sum _{m=0}^n\frac{(n+m+2)!}{m!(n-m)!(m+1)!2^m}{(\mu -1)}^m},\end{array}$$(32)

and for (μ − 1)^m, we have

$$\begin{array}{c}\hfill {(\mu -1)}^m={\displaystyle \sum _{k=0}^m{(-1)}^k\frac{m!}{k!(m-k)!}{\mu }^{m-k}}.\end{array}$$(33)

Using

$$\begin{array}{c}\hfill {\displaystyle \underset{-\pi /2}{\overset{\pi /2}{\int }}{sin}^q(u)\mathrm{d}u=\pi \frac{2^{-q-1}q!({(-1)}^q+1)}{(q/2){!}^2}},\end{array}$$(34)

the integrals in Eq. (31) are given by

$$\begin{array}{cc}& {\displaystyle \underset{-\pi /2}{\overset{\pi /2}{\int }}{P}_n^{(1,1)}(sinu)}\mathrm{d}u\hfill \\ \hfill & =\frac{\pi }{2(n+2)}{\displaystyle \sum _{m=0}^n\frac{(n+m+2)!}{(n-m)!(m+1)!2^{2m}}}{\displaystyle \sum _{k=0}^m\frac{{(-2)}^k}{k!}\frac{({(-1)}^{m-k}+1)}{((m-k)/2){!}^2}}\hfill \end{array}$$(35)

as well as

$$\begin{array}{cc}& {\displaystyle \underset{-\pi /2}{\overset{\pi /2}{\int }}{P}_n^{(1,1)}(sinu){sin}^2}u\mathrm{d}u\hfill \\ \hfill & =\frac{\pi }{8(n+2)}{\displaystyle \sum _{m=0}^n\frac{(n+m+2)!}{(n-m)!(m+1)!2^{2m}}}\hfill \\ \hfill & \times {\displaystyle \sum _{k=0}^m\frac{{(-2)}^k}{k!}\frac{({(-1)}^{m-k+2}+1)(m-k+1)(m-k+2)}{((m-k+2)/2){!}^2}}.\hfill \end{array}$$(36)

Because of the symmetry of the integrand, the integral vanishes for odd values of n, implying β_nk = 0 and α_nk = 0. Using Eq. (35) and Eq. (36) allows us to explicitly calculate the Jacobian integral Eq. (31) (see Table 1). In order to normalize and calculate α_nk, we have to make use of the initial conditions, which correspond to the state of rigid rotation. Formally, we have

$$\begin{array}{c}\hfill \mathrm{\Omega }(r,0,0)={\mathrm{\Omega }}_0,\end{array}$$(37)

which implies

$$\begin{array}{c}\hfill \omega (r,0,0)=0\end{array}$$(38)

and can be satisfied by demanding α_nk(0) = 0.

2.4. Energy dissipation

According to Lanza (2006), the variation of the rotational kinetic energy can be calculated from

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathcal{T}}_{\mathit{nk}}=\frac{8\pi (n+1)}{(2n+3)(n+2)}{\alpha }_{\mathit{nk}}^2{\displaystyle \underset{r_{\mathrm{b}}}{\overset{R}{\int }}\rho r^4{\zeta }_{\mathit{nk}}^2}\mathrm{d}r,\end{array}$$(39)

with a kinetic energy dissipation rate P_diss given by

$$\begin{array}{c}\hfill P_{\mathrm{diss}}=-2{\displaystyle \sum _k}{\displaystyle \sum _n{\lambda }_{\mathit{nk}}\mathrm{\Delta }{\mathcal{T}}_{\mathit{nk}}}.\end{array}$$(40)

2.5. Gravitational quadrupole moment variation

According to Ulrich & Hawkins (1981) and Lanza (2006), the variation of the quadrupole moment potential ΔΦ₁₂ can be calculated by solving

$$\begin{array}{ccc}& {\displaystyle \frac{{\partial }^2(\mathrm{\Delta }{\mathrm{\Phi }}_{12})}{\partial r^2}=}\hfill & \hfill -\frac{2}{r}\frac{\partial (\mathrm{\Delta }{\mathrm{\Phi }}_{12})}{\partial r}+\frac{6}{r^2}(\mathrm{\Delta }{\mathrm{\Phi }}_{12})+\frac{4\pi r^2}{M(r)}\\ \hfill & \times [\left(\frac{\mathrm{d}\rho }{\mathrm{d}r}\right)(\mathrm{\Delta }{\mathrm{\Phi }}_{12})-\frac{\partial }{\partial r}\left(r^2\rho b_2\right)-r\rho a_2],\hfill \end{array}$$(41)

where M(r) is the integrated star mass up to a given radius r. The Lebovitz (1970) coefficients a₂ and b₂ are calculated from the radial eigenfunctions ζ_nk

$$\begin{array}{c}\hfill a_2={\mathrm{\Omega }}_0{\displaystyle \sum _k\frac{12}{7}}{\alpha }_{2k}{\zeta }_{2k}-\frac{4}{3}{\alpha }_{0k}{\zeta }_{0k}\end{array}$$(42)

and

$$\begin{array}{c}\hfill b_2={\mathrm{\Omega }}_0{\displaystyle \sum _k\frac{2}{3}{\alpha }_{0k}{\zeta }_{0k}+\frac{4}{7}{\alpha }_{2k}{\zeta }_{2k}}.\end{array}$$(43)

Solutions for Eq. (41) must verify the condition

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Phi }}_{12}^{\prime }+3\mathrm{\Delta }{\mathrm{\Phi }}_{12}/R=0\end{array}$$(44)

to match the outer gravitational potential (Lanza 2006).

In order to solve Eq. (41) with the boundary condition Eq. (44), we perform a shooting-method search and convert this boundary value problem into an initial value problem. The initial conditions are given by ΔΦ₁₂(r) = C r² and $\mathrm{\Delta }{\mathrm{\Phi }}_{12}^{\prime }(r)=2Cr$ for r → 0. By performing subsequent integrations of the quadrupole moment Eq. (41) for varying trial constants C, we find a trial constant C for which condition (44) holds.

From ΔΦ₁₂(R), we can calculate the quadrupole moment variation ΔQ via (Lanza 2006)

$$\begin{array}{c}\hfill \mathrm{\Delta }Q=-\frac{R^3\mathrm{\Delta }{\mathrm{\Phi }}_{12}(R)}{3G},\end{array}$$(45)

giving a relative orbital period variation of Applegate (1992)

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }P}{P}=-9\frac{\mathrm{\Delta }Q}{M{a}^2},\end{array}$$(46)

where a is the semi-major axis of the binary and M is the mass of the active component.

3. Results

We applied the formalism described in Sect. 2 to an EZAMS² model of HR 1099³ (M = 1.3 M_⊙, solar metallicity, time step 426, t = 4.48 Gyr, R ∼ 4.05 R_⊙, L = 8.62 L_⊙, Ω₀ = 2.569 ⋅ 10⁻⁵ s⁻¹, and depth of the convection zone r_b/R = 0.1968), which provides the necessary stellar data: radial mass profile M(r), luminosity profile L(r), density profile ρ(r), temperature profile T(r), and mean molecular weight m_μ(r), as well as the base of the convection zone r_b. Using this, we can calculate derived quantities such as the radial gravitational acceleration profile g(r), convective velocity u_c(r), pressure scale height h_p(r), and turbulent viscosity η_t(r), for which we employ mixing length theory with a mixing length parameter α_ml = 1.5. We plot the radial profiles of u_c(r), h_p(r) and η_t(r) in Fig. 1.

For the surface magnetic flux density, we assume B_surf = 1 kG. Details on the calculation of the Maxwell stress tensor and Reynolds stress tensor are described in Sect. 2. We adopt A_v = A_B = 0.1 as relative oscillation amplitudes for the velocity field and magnetic field inside the star.

3.1. Analytical approach

In order to understand the temporal evolution of the Applegate effect in the model presented here, we start with an analytical approach toward the most fundamental quantities involved in the mechanism. For a given set of functions that describe the oscillation of stellar parameters, the mechanism is controlled by the functions β_nk and α_nk, which can be calculated from the source term S(r, μ, t), which is in turn controlled by the function ${\tilde{B}}_{i\varphi }(r,t)$ as given by Eq. (29). ${\tilde{B}}_{i\varphi }(r,t)$ consists of three terms:

• a magnetic field term B_iB_ϕ(r, t) that controls the mean magnetic field,

• the Reynolds tensor Λ_iϕ(r, t) that describes the velocity field fluctuations, and

• the Maxwell tensor M_iϕ(r, t) that describes the magnetic field fluctuations.

Let f(t) be a periodic function with period P_act (implying f(t + P_act) = f(t)) that describes the temporal variation of the stellar parameters involved in the Applegate mechanism, that is, $B_i^{\prime }(r,t)=B_i(r)f(t)$), and ${\upsilon }_i^{\prime }(r,t)=u_c(r)f(t)$. Furthermore, we employ a cyclic purely azimuthal magnetic field B_ϕ(t) = B_surf f(t) with constant amplitude B_surf. We impose that all these quantities are tightly coupled and oscillate in phase. Under such conditions, all terms in ${\tilde{B}}_{i\varphi }(r,t)$ have the same time dependence, implying $\tilde{B}(r,t)\propto f{(t)}^2$ for fixed r. Substituting into Eq. (15) yields β_nk ∝ f(t)². Given this result, we can evaluate Eq. (14): Following Lanza (2006) and employing that the first eigenvalues are typically associated with angular momentum redistribution modes on timescales shorter than the activity cycle of the star, we have α_nk = β_nk/λ_nk, which yields α_nk ∝ f(t)².

Finally, we arrive at

$$\begin{array}{c}\hfill P_{\mathrm{diss}}\propto f{(t)}^4.\end{array}$$(47)

The main insight of this exercise is the change in periodicity: In case of sinusoidal variations with period P_act in the magnetic field and the velocity field fluctuations, α_nk, β_nk, and P_diss oscillate with half of the activity period to the leading order⁴ which is the observed binary modulation period P_mod as the quadrupole moment variation closely follows the dissipated energy (see following subsections). Another interesting result can be found by calculating the dimensionless quantity

$$\begin{array}{c}\hfill \frac{\tilde{\mu }\mathrm{\Lambda }}{M}=\tilde{\mu }\frac{\underset{r_{\mathrm{b}}}{\overset{R}{\int }}{A}_{\upsilon }^2{u}_c^2\rho \mathrm{d}r}{\underset{r_{\mathrm{b}}}{\overset{R}{\int }}{A}_B^2{B}_{\mathrm{surf}}^2\mathrm{d}r},\end{array}$$(48)

which compares the Reynolds tensor integral with the Maxwell tensor integral, quantifying the relative strength of the velocity field fluctuation amplitude versus the magnetic field fluctuation amplitude. Given A_v = A_B = 0.1, a surface magnetic field of B_surf = 1 kG, and our EZAMS model of HR 1099, we have $\tilde{\mu }\mathrm{\Lambda }/M\gg 1$, implying that the velocity fluctuations dominate the magnetic field fluctuations. The same holds for the zero-age main-sequence stars investigated in the parameter study in Sect. 3.6. As a result, the choice of B_surf and A_B does not significantly alter our results.

3.2. Eigenvalues and radial eigenfunctions

Figures 2 and 3 show the first four eigenfunctions ζ_0k(r) and ζ_2k(r) associated with the first four eigenvalues of λ_0k and λ_2k (see Table 2). A radial eigenfunction of order k has k nodes in [r_b, R] (see, e.g., Lanza 2006).

3.3. Eigenfunction normalization and k cutoff

For a fixed r_b and assuming that all radial eigenfunction ζ_nk share one normalization parameter ζ_nk(r_b), we find that the scaling relations ω ∝ ζ_nk(r_b) and P_diss ∝ ζ_nk(r_b)² hold. Consider ζ_nk(r) = ζ_nk(r_b) ξ_nk(r) with dimensionless functions ξ_nk(r). Plugging into Eq. (9) and making use of the linearity of differentials yields

$$\begin{array}{c}\hfill \frac{1}{\rho r^4}\frac{\mathrm{d}}{\mathrm{d}r}(r^4{\eta }_{\mathrm{t}}{\xi }_{\mathit{nk}}^{\prime })-n(n+3)\frac{{\eta }_{\mathrm{t}}}{\rho r^2}{\xi }_{\mathit{nk}}+{\lambda }_{\mathit{nk}}{\xi }_{\mathit{nk}}=0,\end{array}$$(49)

which is equivalent to Eq. (9). This implies that for a fixed value of r_b, n, and k, solutions for ζ_nk(r) with varying normalizations ζ_nk(r_b) are multiples of one another. In particular, the eigenvalues λ_nk are constant. With this, we obtain for the dissipated energy

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathcal{T}}_{\mathit{nk}}={\zeta }_{\mathit{nk}}^2(r_{\mathrm{b}})\frac{8\pi (n+1)}{(2n+3)(n+2)}{\alpha }_{\mathit{nk}}^2{\displaystyle \underset{r_{\mathrm{b}}}{\overset{R}{\int }}\rho r^4{\xi }_{\mathit{nk}}^2}\mathrm{d}r,\end{array}$$(50)

and the kinetic energy dissipation rate becomes

$$\begin{array}{c}\hfill P_{\mathrm{diss}}=-2{\displaystyle \sum _k}{\displaystyle \sum _n{\lambda }_{\mathit{nk}}\mathrm{\Delta }{\mathcal{T}}_{\mathit{nk}}\propto {\zeta }_{\mathit{nk}}^2(r_{\mathrm{b}})}\end{array}$$(51)

as the eigenvalues λ_nk are independent of the normalization chosen. In an analogous way, we can write

$$\begin{array}{c}\hfill \omega (r,\mu ,t)={\zeta }_{\mathit{nk}}(r_{\mathrm{b}}){\displaystyle \sum _{n=0}^{\infty }}{\displaystyle \sum _k{\alpha }_{\mathit{nk}}(t){\xi }_{\mathit{nk}}(r)P_n^{(1,1)}(\mu )}\end{array}$$(52)

for the angular velocity variation inside the star, which indeed implies ω(r, μ, t) ∝ ζ_nk(r_b) because the functions α_nk are invariant under variations of ζ_nk(r_b).

As we investigate the question whether the Applegate effect – if operative – can explain the observed period variations in close binaries, we choose the following normalization strategy: First, we start with ζ_nk(r_b) = 1 and calculate the eigenvalues λ_nk and the functions ζ_nk, β_nk and α_nk for n = 0 and n = 2 up to some maximum order k_max. The maximum order is determined by the eigenvalues, which have the dimension of inverse time and are associated with angular momentum redistributions on a timescale of 1/λ_nk (see Lanza 2005; Lanza 2006). For increasing orders of n and k, these timescales gradually shrink until they fall below the typical travel time of a soundwave when crosses the convection zone, which we calculate from

$$\begin{array}{c}\hfill t_{\mathrm{c}}=\frac{{(R-r_{\mathrm{b}})}^2}{\underset{r_{\mathrm{b}}}{\overset{R}{\int }}{c}_{\mathrm{s}}\mathrm{d}r},\end{array}$$(53)

where the speed of sound is

$$\begin{array}{c}\hfill c_{\mathrm{s}}=\sqrt{5P/3\rho }\end{array}$$(54)

assuming monoatomic gas. From this point, we therefore assume that higher orders are not able to contribute to the global angular momentum redistribution process. Depending on the details of the stellar interior, typicalk cutoffs range between 50 for evolved subgiants and >100 for red dwarfs. Given a k cutoff and the eigenfunctions ζ_nk (starting with ζ_nk(r_b) = 1), β_nk as well as α_nk, we calculate the total energy dissipation P_diss(t) and locate its maximum P_diss, max. Given this, we adjust the eigenfunction normalization factor ζ_nk(r_b) in order to have a maximum dissipated power equal to A_P L_HR1099. Using this normalization, the kinetic energy dissipation rate is limited to some fraction A_P of the stellar power output L_HR1099, which is typically on the order of 10%, and we can finally calculate the corresponding angular velocity variations inside the star.

3.4. Angular velocity variation and energy dissipation

The dissipated energy is calculated from Eq. (51) up to some radial order k. For k = 48, the timescale associated with the respective angular momentum redistribution mode is shorter than the calculated soundwave travel time. We fix the kinetic energy fluctuation limit to A_P = 0.1 as described in the previous subsection. When α_nk and ζ_nk are determined and our normalization to limit the kinetic energy dissipation is applied, we can calculate the effective angular velocity variation inside the star as functions of both time and stellar radius. In Fig. 4 we present the angular velocity changes inside the star for three different times between t = 0.05P_mod and t = 0.25P_mod. Starting from the state of rigid rotation at t = 0, the angular velocity variation builds up continuously with peak fluctuations of a few percent in the outermost layers of the star.

3.5. Period modulation

In Fig. 5 we plot the calculated period variation in HR 1099 for B_surf = 1 kG, a cycle length of P_act = 2 P_mod, a magnetic field fluctuation amplitude of A_B = 0.1, and a velocity field fluctuation amplitude A_V = 0.1 for three different cases of phase lag Δϕ between the radial and the azimuthal velocity and magnetic field fluctuations. The dissipated energy is limited to A_P = 0.1 of the star’s luminosity. For Δϕ = 0, the oscillation is purely negative, that is, the quadrupole moment changes have a positive sign. On the other hand, for Δϕ = π/2, the period modulation oscillates between the positive and the negative regime, which is indicative of a cyclic transition between an oblate and a prolate state. Figure 6 shows the expected level of period variation in case of a surface magnetic field of 1 kG and varying fluctuation parameters A_V, A_B, and A_P (see Sect. 2.2.1). Even with fluctuation amplitudes as high as A_V = A_B = A_P = 0.5, which are not observed in HR 1099 (see Frasca & Lanza 2005), the resulting binary period variation amplitudes are still off by two orders of magnitude; no physically reasonable parameter range yields period variations as high as those observed in the system (9.0 × 10⁻⁵, Frasca & Lanza 2005). Given that the redistribution processes are almost entirely governed by convection (see Sect. 3.1), the only degree of freedom left to produce larger period variations is the convective velocity inside the star, for instance, through the mixing-length parameter α_ml. Reproduction of the observed modulation period of 35 yr imposes an activity cycle length of twice that value as the model predicts a 2:1 relation between the observed binary modulation period and the activity cycle length. Therefore, if the observed period variations are energetically and mechanically feasible, we would expect an activity cycle in HR 1099 of roughly 70 yr, while several authors claim to have found cycle lengths between 14.1 ± 0.3 yr and 19.5 ± 2 yr (see Lanza et al. 2006; Muneer et al. 2010; Perdelwitz et al. 2018, and references therein).

3.6. Active component mass parameter study

We showed in the previous sections that RS CVns such as HR 1099 are not expected to produce significant levels of period modulation. We now expand our analysis to pre-cataclysmic PCEB systems (see, e.g., Zorotovic & Schreiber 2013), which have been proposed as potential Applegate systems by Völschow et al. (2016).

Consider a tidally locked binary system with a = 1 R_⊙, consisting of a typical white dwarf primary with M_WD = 0.5 M_⊙ and a zero-age main-sequence secondary with varying mass. Just as in the previous sections, we assume a constant magnetic field of B_surf = 1 kG and fix the fluctuation parameters to A_P = A_V = A_B = 0.1 with a phase lag Δϕ = 0. For the modulation period, we assume P_mod = 15 yr. We varied the mass of the active component between 0.1 M_⊙ and 0.6 M_⊙ and calculated the amplitude of the orbital period variation using an EZAMS model of age t = 0. Figure 7 plots the expected period variation amplitude as a function of the active component mass. Between 0.1 M_⊙ and 0.36 M_⊙, we see a clear increasing trend toward higher binary period variation levels. For fluctuation parameters slightly above the fiducial 10% level, stars between 0.30 M_⊙ and 0.36 M_⊙ are expected to produce amplitudes ≃10⁻⁷, which is a typical order of magnitude observed in PCEB systems (Parsons et al. 2010; Zorotovic & Schreiber 2013; Völschow et al. 2016), while stars outside of this range only support low levels of period variations ≲10⁻⁷. For all simulated systems, the peak angular velocity variations are ≲1%.

The physical explanation for the peak around M ≃ 0.35 M_⊙ is the emergence of a radiative core, which shifts the boundary of the convective zone outward and reduces the total mass and angular momentum involved in the redistribution process. In Fig. 8 we plot the bottom of the convective zone as function of mass. Stars with M ≳ 0.37 M_⊙ develop a radiative core, which reduces the domain available for angular momentum redistribution processes. While the calculated binary period modulation amplitudes again increase toward higher masses, the upper end of the investigated mass range may be affected by Roche-lobe overflow, which introduces additional effects that are not considered in our model. We also recall that fully convective stars with r_b/R = 0 are strictly speaking not covered by the regular Sturm–Liouville problem, but the structure files adopted in our calculations for stars up to M = 0.36 M_⊙ have relative convection zone depths of r_b/R ≃ 0.01.

4. Conclusions and discussion

We presented a model that connects periodic fluctuations in the magnetic field and velocity field of the convection zone to periodic modulations of the quadrupole moment of the star, resulting in cyclic modulations of the binary period. Using the example of HR 1099, which has been discussed as one of the most promising Applegate candidates, we calculated the angular momentum redistribution modes ζ_nk(r) and considered the angular equation to calculate the α_nk(t) and β_nk(t) functions that are required to calculate the temporal evolution of the dissipated energy and the quadrupole moment, assuming that the star dedicates a fiducial fraction of A_P = 0.1 of its luminosity to drive the Applegate mechanism. The order up to which we calculate the eigenfunctions is determined by a simple soundwave travel time argument. Our constant azimuthal field ansatz leads to a simple source function only controlled by the Λ_rϕ and M_rϕ components of the Reynolds and Maxwell tensors, where the former clearly dominates in the case of simple harmonic velocity and magnetic field fluctuations, which we control through the A_B and A_V parameters with fiducial values of 0.1.

In this model, harmonic oscillations of the magnetic field fluctuations and velocity field fluctuations on a timescale of the activity cycle P_act result in orbital period variations on a timescale of P_mod = 0.5 P_act. In other words, we expect that the activity cycle of any Applegate candidate satisfying our assumptions is double the observed binary modulation period, whereas alternative models predict a different ration (see, e.g., Applegate 1992; Lanza et al. 1998; Lanza & Rodonò 2004).

Furthermore, depending on the phase shift between the radial fluctuations and the azimuthal fluctuations, purely negative, positive or mixed quadrupole moment changes may be observed similar to findings by Rüdiger et al. (2002).

In line with previous work, we confirm that HR 1099 is not expected to produce period variations on the level observed today through the Applegate mechanism, which would require unphysically high field fluctuation amplitudes, angular velocity changes, and dissipated power. This is consistent with previous findings by Lanza (2005, 2006), who did not explicitly consider the temporal evolution of the mechanism.

While HR 1099 in particular and RS CVn systems in general seem highly unlikely Applegate candidates, our extensive parameter study identified short-period PCEBs with active component masses between 0.30 M_⊙ and 0.36 M_⊙ as the strongest Applegate candidates with expected period variation amplitudes similar to those typically observed in such systems, which are between 10⁻⁷ and a few 10⁻⁶ (see, e.g., Zorotovic & Schreiber 2013; Völschow et al. 2016). Consequently, we propose a careful reanalysis of PCEB systems with cyclic period variations to exclude the planetary hypothesis and find direct observational evidence that intrinsic effects such as quadrupole moment changes are at work.

To distinguish between Applegate modulations and period variations caused by planets through the light travel time effect (see, e.g., Pribulla et al. 2012), long-term measurements are required to determine the correlations between activity indicators, luminosity variations that are due to photospheric temperature changes, and the binary period modulation (Applegate 1992,). Changes in the rotational profile on the 1% level may lead to tiny deformations on a similar order of magnitude as the active component oscillates between an oblate and a prolate state, which is also expected to correlate with the activity cycle (Applegate 1992; Lanza et al. 1998; Rüdiger et al. 2002,). However, matters are complicated by spots and tidal deformations.

The model presented here is a significant step forward to understand the physics of the Applegate scenario. Nevertheless, it is still based on simplifying assumptions, such as a purely azimuthal magnetic field and the adoption of a mean-field framework. The Reynolds and Maxwell stress tensors do not self-consistently follow from a dynamo model either, but we adopt a (physically motivated) prescription for the kinetic and magnetic fluctuations that is externally imposed. Similarly, this approach still neglects effects such as viscosity quenching as a result of fast rotation, which may reduce the dissipated power and allow for higher levels of period variation in rapidly rotating systems (see Lanza 2006). Furthermore, our model does not guarantee the hydrodynamical stability of the differential rotation profile that emerges during the oscillation, although the maximum relative deviations from the state of rigid rotation are on the 2% level in HR 1099 and well below the 1% level in our PCEB parameter study (see Knobloch & Spruit 1982). Another limitation ist that we do not account for meridional circulation or for the effect of energy advection inside the star. Finally, additional quadrupole moment changes may occur as results of anisotropic Lorentz forces inside the star (Lanza et al. 1998).

While we cannot draw a final conclusion at this point, our framework nevertheless provides the most complete assessment of the Applegate mechanism so far, and it shows that a consideration of physical processes allows determining under which conditions the Applegate mechanism may be feasible. Our model yields testable predictions, suggesting that the Applegate mechanism may be favored for certain PCEB systems and disfavored in HR 1099. These predictions should be tested both on the observational side and through further refinement of the model to provide a clear picture where the Applegate mechanism is able to work.

---

Acknowledgments

We would like to express our gratitude to the anonymous referee for carefully reading our manuscript and providing us with valuable advice that improved the quality of our work. DRGS thanks for funding via Fondecyt regular (project code 1161247), the “Concurso Proyectos Internacionales de Investigación, Convocatoria 2015” (project code PII20150171) and the BASAL Centro de Astrofísica y Tecnologías Afines (CATA) PFB-06/2007.

References

All Tables

All Figures