format: text

text:
# Dynamics of magnetic flux tubes in close binary stars

## 1 Introduction

Owing to their rapid rotation and deep convection zone, evolved (sub-)giants and cool main sequence components in synchronised close binaries like RS CVn and BY Dra systems typically show high levels of magnetic activity (Strassmeier et al. 1993), e.g., in the form of large starspots and enhanced chromospheric and coronal emission in the UV and X-ray spectral ranges. Analyses of light curves with inversion techniques allow to derive the characteristic properties of spots like their number, position, size, and temperature (Eaton & Hall 1979; Rodonò et al. 1986; Budding & Zeilik 1987). Utilising rotation-modulated, asymmetric spectral line profiles of spectro- and spectropolarimetric observations, inversion techniques like the Doppler Imaging (Vogt & Penrod 1983; Vogt et al. 1987) and Zeeman-Doppler Imaging (Semel 1989; Donati et al. 1992, 1999) reveal information about temperature and magnetic field inhomogeneities on the stellar surface, which are interpreted as cool spots caused by emerging magnetic fields. Recent observations are complemented by photometric data which cover up to several decades and thus enable the determination of cyclic variations (Rodonò et al. 1995, 2000).

Surface reconstructions of close, fast-rotating binaries show extended starspots which cover large fractions of the stellar surface (O’Neal et al. 1998; Lanza et al. 2001) and frequently occur, in contrast to the case of the Sun, also at intermediate, high and even polar latitudes. Non-uniform spot distributions in longitude indicate the existence of *preferred longitudes* (PL), where spots occur more frequent or last longer than at other longitudes (Henry et al. 1995; Jetsu 1996). Preferred longitudes are often found to be separated by about $`180\mathrm{°}`$, i.e., in opposite active quadrants of the hemisphere, but their orientation with respect to the companion star is not unique. In short-period systems with rotation periods $`1\mathrm{d}`$ active regions tend to prefer the quadrature longitudes (Zeilik et al. 1989, 1990a, 1990b, 1994; Oláh et al. 1994; Heckert et al. 1998), i.e., the longitudes perpendicular to the line connecting both stellar centres. In systems with longer periods the situation is less clear: several systems show PL at the substellar point, its antipode or other fixed longitudes with respect to the companion star (e.g., Lanza et al. 1998, 2001; Oláh et al. 2002), whereas other stars reveal a longitudinal migration of the spotted regions with respect to the companion star, which is usually ascribed to differential rotation (Strassmeier 1994; Rodonò et al. 1995, 2000). Considering the dependence of PL on the orbital period, Heckert & Ordway (1995) find a transition from fixed PL for orbital periods less that one day to migrating and eventually no PL at all for more separated systems with periods exceeding about 10 days, leading to the suggestion that tidal effects might influence the spot distribution. However, observations of close binaries with giant components show that spot clusters at migrating PL may also occur in systems with longer orbital periods (Berdyugina & Tuominen 1998; Berdyugina et al. 1999), giving rise to the assumption that the formation and orientation of clusters are not exclusively dependent on the binary separation alone.

In this and a subsequent paper we study the influence of tidal effects on the dynamics and evolution of magnetic flux tubes in active binary components, assuming – in analogy to the case of the Sun – that starspots are formed by erupting flux tubes which originate from the bottom of the stellar convection zone. The present work delves into the equilibrium and linear stability properties of flux tubes stored in the convective overshoot region below the convection zone proper in order to examine whether the influence of the companion star is able to trigger rising loops at preferred longitudes which may penetrate to the convection zone above. The decomposition of small perturbations of an equilibrium flux tube in terms of eigenmodes and the determination of the corresponding eigenfrequencies are accomplished by a linear stability analysis which extends earlier treatments by Ferriz-Mas & Schüssler (1993, 1995) and Holzwarth & Schüssler (2000). In contrast to rotating single stars, in binary stars the existence of the companion star breaks the rotational symmetry. Consequently, the dynamics of flux tubes is not only subject to the tidal force but is additionally affected by the deviation of the stellar structure from spherical symmetry. The linear stability analysis deals with the response of equilibrium flux tubes to small perturbations inside the overshoot region, whereas the growing displacements of unstable loops have to be followed by non-linear numerical simulations. The evolution of flux loops rising through the convection zone and the resulting surface distribution upon their eruption is described in a following paper.

In Sect. 2 we briefly summarise the underlying (solar) activity model based on erupting magnetic flux tubes and outline the binary model. Sect. 3 describes the equilibrium properties and Sect. 4 gives the results of the linear stability analysis of toroidal flux tubes. Sect. 5 contains a discussion of the results and Sect. 6 gives our conclusions.

## 2 Model assumptions

### 2.1 The ‘flux tube paradigm’

It is widely accepted that sunspots and bipolar spot groups are caused by erupting magnetic flux tubes which originate from the base of the solar convection zone. The magnetic field is believed to be amplified in the rotational shear layer (tachocline) near the base of the convection zone and stored in the form of toroidal flux tubes in the convective overshoot region (Moreno-Insertis et al. 1992; Schüssler et al. 1994). This stably stratified subadiabatic layer at the interface to the radiative core is generated by gas motions which penetrate from the convection zone above. While Hale’s polarity rules and the East-West-orientation of bipolar sunspot groups suggest equilibrium flux rings parallel to the equatorial plane, overshooting gas motions cause perturbations which can trigger the onset of an undulatory (Parker-type) instability (Spruit & van Ballegooijen 1982; Ferriz-Mas & Schüssler 1993, 1995). If a tube segment is lifted upward, a net mass downflow from the crest increases the density contrast with respect to the environment. Once a critical magnetic field strength is exceeded the resulting buoyancy surpasses the opposing magnetic curvature force, the loop continues to rise through the convection zone, and eventually erupts at the stellar surface (Parker 1955; Moreno-Insertis 1986). Upon emergence it causes the various signatures of magnetic activity like active regions, dark spot groups, bright faculae, chromospheric activity, X-ray bright coronal loops, and flares. The buoyancy-driven mechanism implies a crucial dependence on the stellar stratification, viz. the superadiabaticity, $`\delta =_{\mathrm{ad}}`$. In the case of the Sun, magnetic flux tubes with field strengths of about $`10^5\mathrm{G}`$ are required to obtain theoretical results which are in accordance with observational constrains like flux emergence latitudes, tilt angles, and asymmetries of bipolar spot groups (Moreno-Insertis 1992; D’Silva & Choudhuri 1993; Fan et al. 1994; Schüssler et al. 1994; Caligari et al. 1995, 1998). The quantitative elaboration of the applied flux tube model is also consistent with the emergence of polar spots on rapidly rotating young stars (Granzer et al. 2000) and with the conspicuous decline of coronal X-ray emission in giant stars across the ‘coronal dividing line’ (Holzwarth & Schüssler 2001).

We consider the magnetic flux in the convection zone in the form of isolated magnetic flux tubes, which are embedded in a field-free environment. The stellar convection zone is treated as an ideal plasma with vanishing viscosity and infinite conductivity, implying the concept of frozen-in magnetic flux. All calculations are carried out in the framework of the ‘thin-flux-tube approximation’ (Roberts & Webb 1978; Spruit 1981): the diameter of a circular flux tube is assumed to be small with respect to the other relevant length scales of the problem like the pressure scale height, radius of curvature, or the wave length of a perturbation along the tube. As a consequence of very short signal travel times it remains in lateral (total) pressure equilibrium with its environment.

### 2.2 Binary model

Observations indicate that the orbital motion of binary systems with periods $`T10\mathrm{d}`$ is effectively synchronised with stellar rotation and typically exhibit very small eccentricity (Duquennoy & Mayor 1991; Strassmeier et al. 1993). Hence, we consider a close, detached, synchronised and circularised binary system, for which the orbital and stellar rotation periods are equal and the spin axes of the components are perpendicular to the orbital plane. The active star with mass $`M_{}`$ and the companion star with mass $`M_{\mathrm{co}}=qM_{}`$ move on circular orbits with a period of a few days. According to Kepler’s third law, they are separated by the constant distance

$$\left(\frac{a}{R_{\mathrm{}}}\right)4.21\left(1+q\right)^{1/3}\left(\frac{M_{}}{M_{\mathrm{}}}\right)^{1/3}\left(\frac{T}{\text{d}}\right)^{2/3}.$$
(1)
Our *reference binary system* consists of two solar-type stars ($`M_{}=M_{\mathrm{co}}=1M_{\mathrm{}}`$, i.e. stellar mass ratio $`q=1`$) with a period of $`T=2\mathrm{d}`$. The active star is described by a perturbed single star (here, standard solar) model, whereas the companion star is considered as an idealised point mass. In the case of the Sun, the differential rotation profile affects the stability properties only slightly (e.g., Ferriz-Mas & Schüssler 1993). Owing to tidal interactions and synchronisation, the differential rotation in close binary stars is much weaker than in the solar case (Donahue 1993; Rüdiger et al. 1998), so that we assume solid-body rotation.

For systems with orbital periods of a few days, tidal effects are sufficiently small to be treated in lowest-order perturbation theory. Using spherical coordinates in a co-rotating frame of reference (radius $`r`$, latitude $`\lambda `$, azimuth $`\varphi `$, with the companion star in the direction $`\varphi =0`$), the effective potential, $`\mathrm{\Psi }_{\mathrm{eff}}=\mathrm{\Psi }_{}+\mathrm{\Psi }_{\mathrm{tide}}`$, is the sum of the spherical gravitational potential of the star, $`\mathrm{\Psi }_{}(r)`$, and the potential of the tidal perturbation111Equations (2) and (3) also contain the centrifugal potential and acceleration, respectively, although these contributions do not strictly belong to tidal interactions but occur in rotating single stars as well.,

$`\mathrm{\Psi }_{\mathrm{tide}}(r,\varphi ,\lambda )`$ $`=`$ $`ϵ^3{\displaystyle \frac{GM_{}}{2r}}[(1+{\displaystyle \frac{5}{2}}q)\mathrm{cos}^2\lambda q`$ (2)
$`+{\displaystyle \frac{3}{2}}q\mathrm{cos}^2\lambda \mathrm{cos}2\varphi ]+𝒪(ϵ^4),`$
where $`ϵ=r/a`$ is the expansion parameter and $`G`$ the gravitational constant. The potential $`\mathrm{\Psi }_{\mathrm{tide}}`$ gives rise to the tidal acceleration

$`g_{\mathrm{tide}}(r,\varphi ,\lambda )`$ $`=`$ $`ϵ^3g_{}(r)[e_r+3q(e_re_a)e_a`$ (3)
$`(1+q)(e_re_\mathrm{\Omega })e_\mathrm{\Omega }]+𝒪(ϵ^4),`$
where $`e_r,e_\mathrm{\Omega }`$, and $`e_a`$ are unit vectors according to the directions indicated in Fig. 1.

Following Eq. (3), the component of the tidal acceleration in azimuthal direction $`e_\varphi `$, for example, is

$$\left(g_{\mathrm{tide}}e_\varphi \right)=\frac{3}{2}g_{}(r)ϵ^3q\mathrm{cos}\lambda \mathrm{sin}2\varphi +𝒪\left(ϵ^4\right).$$
(4)
Our model includes the deviation of the stellar structure from spherical symmetry arising from the tidal acceleration $`g_{\mathrm{tide}}`$. For a not to strongly deformed star, we suppose that the stellar structure is still described appropriately by functions of the effective potential, $`\mathrm{\Psi }_{\mathrm{eff}}`$, only, i.e. $`f(r,\varphi ,\lambda )=f(\mathrm{\Psi }_{\mathrm{eff}})`$, where $`f`$ stands for any stellar quantity (like gas pressure, density,$`\mathrm{}`$). This permits an analytical expression for the tidal deformation, namely

$`f(r,\varphi ,\lambda )`$ $`=`$ $`f(\mathrm{\Psi }_{}+\mathrm{\Psi }_{\mathrm{tide}})f(\mathrm{\Psi }_{})+{\displaystyle \frac{\mathrm{d}f}{\mathrm{d}\mathrm{\Psi }}}|_r\mathrm{\Psi }_{\mathrm{tide}}`$ (5)
$`=`$ $`f_0\left(1+{\displaystyle \frac{r}{H_f}}{\displaystyle \frac{\mathrm{\Psi }_{\mathrm{tide}}r}{GM_{}}}\right)`$
$`=`$ $`f_0\left(1+\overline{r}{\displaystyle \frac{r}{H_f}}+\widehat{r}{\displaystyle \frac{r}{H_f}}\mathrm{cos}2\varphi \right),`$
with the dimensionless coefficients

$`\overline{r}`$ $`=`$ $`ϵ^3{\displaystyle \frac{1}{2}}\left[\left(1+{\displaystyle \frac{5}{2}}q\right)\mathrm{cos}^2\lambda q\right],`$ (6)
$`\widehat{r}`$ $`=`$ $`ϵ^3{\displaystyle \frac{3}{4}}q\mathrm{cos}^2\lambda ,`$ (7)
and the scale height $`H_f=(\mathrm{d}\mathrm{ln}f/\mathrm{d}r)^1`$. Equation (5) includes the rotational flattening and the tidal (ellipsoidal) deformation, with tidal bulges centred on the line connecting the two stars. The deviation of a quantity, $`\mathrm{\Delta }f=f(r,\varphi ,\lambda )f_0(r)`$, from the single star stratification, $`f_0`$, depends on its local scale height, $`H_f`$. Both values, $`f_0`$ and $`H_f`$, are taken from the spherical model at radius $`r`$. The quality of the approximation can be estimated by considering the ratio between the geometrical deformation,

$$\frac{\delta r}{r}\frac{\mathrm{\Psi }_{\mathrm{tide}}}{\mathrm{\Psi }_{}}=\overline{r}+\widehat{r}\mathrm{cos}2\varphi +𝒪(ϵ^4),$$
(8)
and the *smallest* scale height $`H_f`$ of all quantities $`f`$. Since in the overshoot region, the supposed storage location of magnetic flux, the scale height of the superadiabaticity, $`H_\delta `$, is considerably smaller than of any other quantity (see Tab. 1), the tidal influence on the stability properties of flux tubes is mainly controlled by the variation of $`\delta `$ (Holzwarth & Schüssler 2000).

For our reference system ($`T=2\mathrm{d},M_{}=M_{\mathrm{}}`$ and $`q=1`$), Eq. (1) yields a separation $`a8R_{\mathrm{}}`$ and $`ϵ^310^3`$. The tidal effects are of small magnitude only, but they imply a significant *qualitative* difference, since both Eqs. (2) and (3) exhibit a periodic azimuthal dependence, which breaks the axial symmetry present in the single star problem.

## 3 Stationary equilibrium

We consider magnetic flux rings in mechanical stationary equilibrium, which are embedded in the overshoot region and parallel to the equatorial plane. The tubes are non-buoyant and in lateral pressure equilibrium with the environment. The inward directed magnetic curvature force is compensated by the outward directed Coriolis force, which arises from a prograde flow inside the tube with a velocity $`v`$ in excess of the stellar solid-body rotation. The equilibrium tube forms a torus with constant radius of curvature, $`R_0=r_0\mathrm{cos}\lambda _0`$, where $`r_0`$ and $`\lambda _0`$ are the radial and latitudinal equilibrium positions, respectively. Since the flux ring is situated in a tidally deformed star (Fig. 2),

the internal density, $`\rho `$ ($`=\rho _\mathrm{e}`$, where $`\rho _\mathrm{e}`$ is the density of the environment), depends, according to the approximation in Eq. (5), on the azimuth as

$$\frac{\rho (\varphi )}{\rho _0}=1+\overline{r}_0\frac{r_0}{H_{\rho 0}}+\widehat{r}_0\frac{r_0}{H_{\rho 0}}\mathrm{cos}2\varphi +𝒪(ϵ^4).$$
(9)
This entails a periodic variation of the magnetic field strength,

$$\frac{B(\varphi )}{B_0}=\frac{1M_{\alpha 0}^2}{\left(\frac{\rho (\varphi )}{\rho _0}\right)M_{\alpha 0}^2}\left(\frac{\rho (\varphi )}{\rho _0}\right),$$
(10)
and internal flow velocity,

$$\frac{v(\varphi )}{v_0}=\frac{1M_{\alpha 0}^2}{\left(\frac{\rho (\varphi )}{\rho _0}\right)M_{\alpha 0}^2},$$
(11)
where $`M_{\alpha 0}=v_0/c_{\mathrm{a0}}`$ is the *Alfvénic Mach number*, i.e., the flow velocity in units of the Alfvén velocity, $`c_{\mathrm{a0}}=B_0/\sqrt{4\pi \rho _0}`$. Quantities with index ‘0’ refer to a comparable axisymmetric equilibrium flux tube in a slowly rotating single star, where the tidal and centrifugal acceleration as well as the deformation of the stellar structure are negligible (e.g., Moreno-Insertis et al. 1992). For flux tubes in mechanical equilibrium, these quantities obey the relation

$$v_0=R_0\mathrm{\Omega }\left[\sqrt{1+\left(\frac{c_{\mathrm{a0}}}{R_0\mathrm{\Omega }}\right)^2}1\right]>0.$$
(12)
In the following, we consider flux rings in the middle of the overshoot region. Table 1 gives values of the characteristic parameters of this reference configuration for a model of a $`1M_{\mathrm{}}`$ star.

## 4 Linear stability analysis

### 4.1 Determination and interpretation of eigenmodes

Once all quantities of the external stratification are determined by the equilibrium depth $`r_0`$, the stability properties of the flux ring depend only on its latitude, $`\lambda _0`$, and magnetic field strength, $`B_0`$. Lagrangian displacements are decomposed in eigenmodes,

$$\mathrm{\Xi }(\varphi ,t)=\xi (\varphi )\mathrm{exp}\left(i\omega t\right),$$
(13)
characterised by the eigenfrequency $`\omega `$ and the eigenfunction $`\xi (\varphi )=(\xi _\mathrm{t},\xi _\mathrm{n},\xi _\mathrm{b})^T`$ (in the co-moving trihedron with the local tangent, normal, and binormal unit vectors, $`[t,n,b]`$, respectively). For unstable eigenmodes with increasing amplitudes $`\tau =1/|\mathrm{}(\omega )|`$ provides the characteristic growth time. The displacement vector $`\xi `$ is determined by the linearised equations of motion,

$$\xi ^{\prime \prime }+\left(_\varphi +i\omega _{\varphi t}\right)\xi ^{}+\left(_\xi +i\omega _t\omega ^2_{tt}\right)\xi =0,$$
(14)
where primes denote derivatives with respect to the azimuthal coordinate, $`\varphi `$ (see Appendix A). Taking the tidal force, Eq. (3), and the azimuthal variation of the equilibrium quantities, Eqs. (9) – (11), into account, the coefficient matrices depend periodically on $`\varphi `$:

$`(\varphi )`$ $`=`$ $`(\varphi +\pi )`$ (15)
$`=`$ $`_0+ϵ^3_\mathrm{c}\mathrm{cos}2\varphi +ϵ^3_\mathrm{s}\mathrm{sin}2\varphi +𝒪(ϵ^4).`$
Closed flux tubes require $`\xi (\varphi +2\pi )\xi (\varphi )`$ and $`\xi ^{}(\varphi +2\pi )\xi ^{}(\varphi )`$, suggesting the ansatz

$`\xi (\varphi )`$ $`=`$ $`\widehat{\xi }(\varphi )\mathrm{exp}\left(im\varphi \right)`$ (16)
$`=`$ $`\left[{\displaystyle \underset{k=\mathrm{}}{\overset{\mathrm{}}{}}}\widehat{\xi }_k\mathrm{exp}\left(ik\varphi \right)\right]\mathrm{exp}\left(im\varphi \right).`$ (17)
In analogy to the single-star problem, a phase factor with the azimuthal wave number $`m`$ has been separated. Due to the $`\varphi `$-dependent contributions in Eq. (15), the substitution of Eq. (17) into Eq. (14) yields the 3-term recursion formula

$$_{m+k}\widehat{\xi }_{k2}+𝒞_{m+k}\widehat{\xi }_k+_{m+k}\widehat{\xi }_{k+2}=0,k.$$
(18)
In contrast to the single-star problem, where an eigenmode corresponds to a single azimuthal wave number $`m`$, i.e., $`\widehat{\xi }_0=\mathrm{const}.`$ and $`\widehat{\xi }_k=0`$ for $`k0`$, an eigenmode here consists of a spectrum of coupled wave modes $`n=m+k`$ with $`k=0,\pm 2,\pm 4,\mathrm{}`$, whose amplitudes $`\widehat{\xi }_k`$ are recursively related to each other through Eq. (18). This relation is illustrated in Fig. 3.

Since the coupling is weak \[$`𝒪(ϵ^3)`$\], the amplitudes decrease quite fast with increasing $`|k|`$, so that more distant wave modes (here, with $`|k|4`$) do practically not contribute to the eigenmode. In most cases the amplitude $`\widehat{\xi }_0`$ is much larger than $`\widehat{\xi }_{\pm 2}`$, so that the resulting eigenmode can still be characterised by the wave number $`m`$.

The linear eigenfunctions characterise the structure of small displacements of flux tubes in the overshoot region *prior* to their rapid, non-linear rise through the convection zone. Owing to the stellar stratification, the most important aspect of an eigenfunction is its component in the radial direction, $`\widehat{\xi }_r=(\widehat{\xi }e_r)`$: the larger the radial displacements, the larger are the changes in the environment of the corresponding tube segments, which become located in either more stable or more unstable layers. Assuming that larger radial displacements toward more unstable layers favour the penetration of a growing flux loop to the convection zone, we use the envelope, $`|\widehat{\xi }_r(\varphi )|`$, of the radial component of the (normalised) eigenfunction as a measure for the onset of rapidly rising flux loops at azimuth $`\varphi `$. Since for eigenmodes with $`\mathrm{}(\omega )<0`$ the magnitude of the radial envelope increases exponentially in time,

$$|\mathrm{\Xi }_r(\varphi ,t)|=|\widehat{\xi }_r(\varphi )|\mathrm{exp}\left[\mathrm{}(\omega )t\right],$$
(19)
an unstable flux tube in the overshoot region will eventually enter the superadiabatical region. Its actual displacement depends, however, also on the azimuthal wave number $`m`$, the frequency $`\mathrm{}(\omega )`$, and the phase relation

$$\mathrm{arg}\mathrm{\Xi }_r(\varphi ,t)=m\varphi +\mathrm{arg}\widehat{\xi }_r(\varphi )+\mathrm{}(\omega )t,$$
(20)
corresponding to an azimuthal propagation of the tube’s crest along the envelope, Eq. (19), with the phase velocity

$$v_\mathrm{p}(\varphi )=\frac{\mathrm{}(\omega )R_0}{m}\left(1+\frac{1}{m}\frac{\mathrm{d}}{\mathrm{d}\varphi }\mathrm{arg}\widehat{\xi }_r(\varphi )\right)^1.$$
(21)
The radial envelope thus represents a kind of ‘statistical weight’, which is applicable either to a large ensemble of comparable flux tubes with similar equilibrium properties and phase shifts evenly distributed in the range $`0\mathrm{arg}\mathrm{\Xi }_r<2\pi `$, or to an individual flux tube with a growth time much longer than the wave frequency, $`\tau 2\pi /\mathrm{}(\omega )`$. If an eigenmode is dominated by the wave mode $`m`$ and only slightly modified by the adjacent wave modes $`m\pm 2`$, i.e., $`\widehat{\xi }_{r,\pm 2}\widehat{\xi }_{r,0}:=1`$, the radial envelope is approximately given by

$`\left|\widehat{\xi }_r(\varphi )\right|`$ $`=`$ $`\left|\widehat{\xi }_{r,0}+\widehat{\xi }_{r,2}e^{i2\varphi }+\widehat{\xi }_{r,2}e^{i2\varphi }+\mathrm{}\right|`$ (22)
$``$ $`1+\left|\widehat{\xi }_{r,2}\right|\mathrm{cos}\left(2\varphi +\alpha _2\right)+\left|\widehat{\xi }_{r,2}\right|\mathrm{cos}\left(2\varphi \alpha _2\right)`$
$`+𝒪(\widehat{\xi }_{r,\pm 2}^2),`$
with $`\alpha _{\pm 2}=\mathrm{arg}\widehat{\xi }_{r,\pm 2}`$. The eigenfunctions thus inherit the properties of the underlying problem and exhibit a $`\pi `$-periodicity in azimuthal direction. Assuming that the shape of $`|\widehat{\xi }_r(\varphi )|`$ represents a measure for the onset of rising loops, Eq. (22) implies the existence of *preferred longitudes* at the maxima of $`|\widehat{\xi }_r|`$.

### 4.2 Growth times and eigenfrequencies

The influence of the companion star modifies the stability properties of flux tubes by breaking the axial symmetry. The resulting *quantitative* difference in growth times $`\tau `$ with respect to a rotating single star is very small; as illustrated in Fig. 4, the relative deviations in all cases are considerably smaller than $`1\%`$.

*Qualitative* differences occur, however, for equilibrium configurations which had been stable in the single-star problem. In previously stable parameter domains new types of instabilities with very long growth times lead to an ‘instability background’ consisting of extended ‘plateaus’ which are pervaded by narrow ‘ridges’ (see Fig. 5).

The ‘plateaus’ are due to the tidal force and the azimuthal variation of the density along the equilibrium flux ring. Using the lateral (total) pressure equilibrium and the variation of the magnetic field strength, Eq. (10), it can be shown that the internal specific entropy, $`S`$, is also a function of $`\varphi `$ and varies approximately as

$`{\displaystyle \frac{1}{c_p}}{\displaystyle \frac{S}{\varphi }}`$ $``$ $`{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\mathrm{d}\mathrm{ln}p_\mathrm{e}}{\mathrm{d}\varphi }}{\displaystyle \frac{ϵ^3}{\beta }}{\displaystyle \frac{r}{H_{p,\mathrm{e}}}},`$ (23)
where $`\beta =8\pi p/B^2`$ is the ratio between the gas pressure and the magnetic pressure, and $`H_{p,\text{e}}`$ is the local pressure scale height. Although we have $`r/H_{p,\text{e}}9`$, the azimuthal variation is nevertheless small compared to the radial variation of the stellar stratification, since the relevant field strengths considered here lead to values in the regime $`\beta 10^4\mathrm{}10^61`$. Due to its velocity excess, $`v`$, a gas element inside the tube would have to change its entropy to preserve the equilibrium conditions. But since there is no mechanism available to accomplish this entropy change, the azimuthal variation leads to a loss of the mechanical equilibrium and to very slow monotonic displacements of the flux tube out of its toroidal configuration. For fast-rotating stars, the internal flow velocity $`v`$ is small and thus the growth times of plateau instabilities are typically several $`10^4`$ days (up to $`\tau 10^7\mathrm{d}`$). The effects of the azimuthal density variation and tidal force along the tube are negligible because the stability properties are dominated by the azimuthal variation of the superadiabaticity,

$$\delta (\varphi )=\delta (r_0)\left(1+\frac{\delta r}{H_\delta }\right),$$
(24)
owing to its relative small scale height in the overshoot region (with $`H_\delta /H_p0.08`$, see Tab. 1). If all azimuthal variations other than that of $`\delta `$ are omitted in the stability analysis, the plateau instabilities vanish and only the ridge instabilities remain.

The ‘ridge’ instabilities with typical growth times of several thousand days owe their existence to a resonant coupling of wave modes with different wave numbers. As an example, the top panel of Fig. 6 shows the growth times of both Parker-type and background instabilities along a cut at $`\lambda =1.64\mathrm{°}`$ in Fig. 5.

The panels below show the frequencies $`\mathrm{}(\omega )`$ of the wave modes involved. In the single-star problem (dotted lines, dashed-dotted if unstable) the wave modes corresponding to different values of the azimuthal wave number, $`m`$, are independent. But the eigenmodes in the binary problem consist of a spectrum of *coupled* wave modes (see Fig. 3) with wave numbers $`n=m,m\pm 2,m\pm 4,\mathrm{}`$: modes of wave number $`m`$ with frequency $`\mathrm{}(\omega )`$ stimulate excitations of the wave modes $`m\pm 2`$, which, supposed $`\mathrm{}(\omega )`$ is close to the eigenfrequency of the adjacent wave mode, results in a resonant amplification and eventually in the onset of an instability (solid lines in Fig. 6, dashed if unstable).

In the following, we dismiss instabilities with very long growth times $`\tau >10^4\mathrm{d}(27\mathrm{years})`$, since they are not expected to play any role for the eruption of magnetic flux tubes in active stars.

### 4.3 Eigenfunctions

The eigenfunctions, $`\widehat{\xi }`$, inherit the property of the underlying problem and exhibit a $`\pi `$-periodicity in azimuthal direction, implying a $`\varphi `$-dependent probability that a loop penetrates to the superadiabatic part of the convection zone and rises rapidly toward the surface. The preferred longitudes, $`\varphi _{\mathrm{max}}`$ (and $`\varphi _{\mathrm{max}}+\pi `$), of loop penetration are determined by the maxima of the radial envelope, $`|\widehat{\xi }_r|_{\mathrm{max}}=|\widehat{\xi }_r(\varphi _{\mathrm{max}})|`$. They are shown in Fig. 7 as a function of the equilibrium latitude $`\lambda _0`$ and field strength $`B_0`$.

Background instabilities with growth times exceeding $`10^4\mathrm{d}`$ have been neglected, leaving behind only a fraction of a ridge at low latitudes and a plateau at $`\lambda 60\mathrm{°}`$ with $`\tau 310^3\mathrm{d}`$. This high-latitude feature is caused by the slow loss of equilibrium due to the internal flow described above, which exhibits only small displacements in radial direction (here, $`|\widehat{\xi }_r|/|\widehat{\xi }_\varphi |𝒪(10^2)`$). In contrast to the case of stable eigenfunctions, for which the envelope maxima are located exactly at the ‘symmetry points’, $`Q_{1,2}`$ or $`V_{1,2}`$, the orientation of $`|\widehat{\xi }_r|_{\mathrm{max}}`$ of unstable eigenfunctions show a phase shift with respect to these points. For instabilities at equatorial and low latitudes dominated by the wave number $`m=1`$, the maxima are located in the vicinity of $`V_{1,2}`$, i.e., around the line connecting both stellar centres, whereas for modes with $`m=2`$ they are distributed in a broad interval around the quadrature points $`Q_{1,2}`$; preferred longitudes of loop penetration thus depend on the equilibrium parameters $`(B_0,\lambda _0)`$ as well as on the dominating wave mode.

The radial envelope, $`|\widehat{\xi }_r|`$, is used as a measure for the probability of loop penetration into the convection zone, subject to the normalisation $`|\widehat{\xi }_r|_{\mathrm{max}}:=1`$. The larger the relative peak-to-peak variation,

$$\mathrm{\Delta }|\widehat{\xi }_r|:=\frac{|\widehat{\xi }_r|_{\mathrm{max}}|\widehat{\xi }_r|_{\mathrm{min}}}{|\widehat{\xi }_r|_{\mathrm{max}}},$$
(25)
the larger is the possibility for the transition of a loop from the overshoot region to the convection zone around $`\varphi _{\mathrm{max}}`$. Figure 8 shows $`\mathrm{\Delta }|\widehat{\xi }_r|`$ as a function of $`\lambda _0`$ and $`B_0`$.

For $`B10^5\mathrm{G}`$ the peak-to-peak variation is $`5\%`$ at intermediate latitudes and up to about $`20\%`$ near the equator. Typically $`\mathrm{\Delta }|\widehat{\xi }_r|`$ decreases with increasing field strengths since the stability properties strongly depend on $`\beta \delta \delta /B^2`$. High field strengths reduce the effect of the azimuthal $`\delta `$-variation, so that magnetic flux tubes with smaller field strengths are, in general, more susceptible to tidal effects. Particularly for ridge instabilities and instability islands at low field strengths, $`\mathrm{\Delta }|\widehat{\xi }_r|`$ is significant and can even reach values close to unity, so that radial displacements are almost suppressed near the envelope minima. The approximation used in Eq. (22) is not applicable here, since the dominance of the wave mode $`\widehat{\xi }_{r,0}`$ is not always given. Ridge instabilities are due to resonant interactions of adjacent wave modes (see Sec. 4.2) with comparable amplitudes while the amplitudes of all other modes are small. The case $`|\widehat{\xi }_{r,0}|,|\widehat{\xi }_{r,2}|\widehat{\xi }_{r,k},k0,2`$, for example, yields the radial envelope

$$|\widehat{\xi }_r|^21+2\frac{|\widehat{\xi }_{r,0}||\widehat{\xi }_{r,2}|}{|\widehat{\xi }_{r,0}|^2+|\widehat{\xi }_{r,2}|^2}\mathrm{cos}\left(2\varphi +\alpha _2\alpha _0\right),$$
(26)
with $`\alpha _{0/2}=\mathrm{arg}\widehat{\xi }_{r,0/2}`$. For $`|\widehat{\xi }_{r,0}||\widehat{\xi }_{r,2}|`$, the radial envelope drops to very small values at longitudes which are determined by the phase factors $`\alpha _0`$ and $`\alpha _2`$. The low magnetic field strengths corresponding to these particular instabilities imply, however, long growth times which are probably not relevant for the evolution of rising flux loops in stars.

### 4.4 Parameter dependence

We examine the dependence of the stability properties on tidal effects in various binary systems which are characterised by the orbital period, $`T`$, the mass ratio, $`q`$, and the mass of the primary star, $`M_{}`$. Equations (3) and (5) show that the magnitude of the azimuthal variation caused by tidal effects is governed by the term

$`ϵ^3q`$ $`=`$ $`{\displaystyle \frac{4\pi ^2}{G}}{\displaystyle \frac{q}{1+q}}{\displaystyle \frac{r^3}{M_{}T^2}}`$ (27)
$``$ $`10^2\left({\displaystyle \frac{q}{1+q}}\right)\left({\displaystyle \frac{r}{R_{\mathrm{}}}}\right)^3\left({\displaystyle \frac{M_{}}{M_{\mathrm{}}}}\right)^1\left({\displaystyle \frac{T}{\mathrm{d}}}\right)^2.`$
Here, the solar-type model is retained, which fixes $`M_{}`$ and $`r_0`$. The variation of Eq. (27) due to a different value of $`r_0`$ is negligible since the equilibrium radius is restricted to the rather thin overshoot layer, which is prescribed by the stellar model.

First we consider binary systems with orbital periods in the range $`T=0.75\mathrm{}10\mathrm{d}`$, mass ratio $`q=1`$, and, following Eq. (1), separations $`a4.4\mathrm{}25R_{\mathrm{}}`$. Since the expansion parameter covers the range $`ϵ^3=4.510^3\mathrm{}2.510^5`$, the analytical approximation for the stellar stratification, Eq. (5), is well applicable even for the shortest period. A comparison of growth times for different $`T`$ is shown in Fig. 9.

A smaller orbital period has a stabilising effect, which is indicated by the shifting of corresponding isolines of growth times to higher field strengths. This effect is reinforced by the rotational flattening of the star since flux rings at a given, constant equilibrium depth, $`r_0`$, are embedded in slightly deeper and more subadiabatic layers of the overshoot region. The tidal effects on $`\tau `$ are nevertheless small and primarily affect the background instabilities, leading to shifts of ridges and plateaus. With decreasing system period, the high-latitude plateau grows toward smaller field strengths since the azimuthal variation of the density and tidal force along the flux ring become larger; however, the growth times remain very large. The orientations of the radial envelope maxima do not show a significant dependence on orbital period and are well represented by the reference case with $`T=2\mathrm{d}`$ (Fig. 7). The relative peak-to-peak variations, $`\mathrm{\Delta }|\widehat{\xi }_r|`$, are shown in Fig. 10.

$`\mathrm{\Delta }|\widehat{\xi }_r|`$ increases significantly with decreasing orbital period, which results in strongly preferred longitudes of loop penetrations to the convection zone in short-period systems. The tidal effects are very effective for flux rings with small field strengths, particularly for those which are liable to ridge instabilities with $`\mathrm{\Delta }|\widehat{\xi }_r|`$ close to unity. For Parker-type instabilities with short growth times, the relative peak-to-peak variation remains much smaller. Figure 10 also shows that the peak-to-peak variation strongly decreases with increasing system period for $`T3\mathrm{d}`$ and eventually becomes insignificant.

Considering a system with period $`T=2\mathrm{d}`$, we find that a variation of the mass ratio in the range $`q=0.1\mathrm{}10`$, has only marginal influence on the stability properties. In the domain of Parker-type instabilities, the relative changes of the growth time and of the peak-to-peak variation with respect to the reference case are well below $`0.1\%`$ and $`10\%`$, respectively. Furthermore, the orientations of the envelope maxima exhibit hardly any dependence on the mass ratio.

In summary, the growth times show only a weak dependence on the binary parameters, except for the instability background and the rotational stabilisation of flux tubes at short system periods. The latter effect is mainly due to the conservation of angular momentum and not innate to the tidal effects. The relative peak-to-peak variation of the radial envelope exhibits a considerable dependence on the orbital period $`T`$, while the dependence on the mass ratio $`q`$ is almost negligible. Finally, the orientations of the envelope maxima are largely unaffected by variations of the binary parameters.

## 5 Discussion

Tidal effects and high rotation rates in close binaries cause characteristic changes of both the equilibrium and the stability properties of magnetic flux tubes stored in the convective overshoot region. Owing to tidal effects, a flux ring in stationary equilibrium experiences periodic variations of its environment in azimuthal direction, which modify its characteristic eigenmodes. To first order, the azimuthal variation of a given quantity depends on its local scale height. Since the stability properties are dominated by the superadiabaticity, which is the most strongly depth-dependent quantity in the overshoot region, tidal effects are mediated mainly by the azimuthal variation of $`\delta `$.

The growth times of Parker-type instabilities with $`\tau 1000\mathrm{d}`$ are only marginally affected by the presence of the companion star, whereas for instabilities with $`\tau 1000\mathrm{d}`$ the influence of tidal effects is considerably stronger. This applies particularly to the instability background, since ridge and plateau instabilities owe their existence to the binary character and are most susceptible to changes of the tidal effects. However, since the eigenfrequency is a global parameter describing the whole flux tube and thus representing an average over the properties of all tube segments, the dependence of $`\omega `$ on the binary parameters is dominated by the axially symmetric influence of stellar rotation and not so much by the non-axially symmetric tidal interactions: for axially symmetric changes of the equilibrium configuration, e.g., those due to the rotational flattening, the changes are identical for each tube segment and thus may add up to a considerable deviation of $`\omega `$, whereas for non-axial, periodic changes along the flux ring, like those due to tidal effects, the changes of individual tube segments mutually cancel and eventually average to only minor deviations of $`\omega `$. The growth times of background instabilities are typically several thousand days and more and therefore very long compared to the system period as well as to the growth times which are relevant for rising flux tubes in the Sun. Whether these instabilities are important for loop formation at preferred longitudes seems questionable, although they possess the most conspicuously asymmetric eigenfunctions.

In contrast to the modest effect on the eigenfrequencies and growth times, tidal effects can significantly alter the azimuthal structure of the eigenfunctions since the periodic variations along the tube affect the displacement of each segment individually. The resulting eigenfunction, $`\xi (\varphi )`$, characterised by its radial envelope $`|\widehat{\xi }_r|`$, consists of a superposition of coupled wave modes and inherits the $`\pi `$-periodicity of the azimuthal variation of the equilibrium properties. In contrast to stable eigenmodes, whose envelope maxima are located either at the line of centres or at the quadrature points, the maxima of unstable eigenmodes typically exhibit a phase offset with respect to these points; this result is schematically illustrated in Fig. 11.

The orientations of the maxima are not unique but depend on the equilibrium configuration, predominantly on its latitude and on the wave number of the dominating wave mode. Considering Fig. 2, one might expect that the quadrature points would be favoured for rapid loop growth, since they are located in less stable layers of the overshoot region and nearest to the superadiabatic convection zone. However, the Parker-type instability is driven by a downflow from the loop summit, which increases the density contrast with respect to the environment. Following the equation of continuity, a net downflow implies an asymmetry of mass flux at the loop summit, i.e., the amount of plasma flowing down in the preceding slope of the loop outbalances the mass flux coming up in the following loop slope. The azimuthal variation of the internal equilibrium density and velocity, Eqs. (9) and (11), enhance the asymmetry in both slopes and thereby favour maximal displacements at intermediate longitudes. In contrast, the extrema of stable eigenmodes prefer either the quadrature points or the directions along the line connecting both stellar centres since at these ’symmetry points’ (concerning the internal equilibrium density and velocity structure) the downflow and thus buoyancy effects can be minimised.

The peak-to-peak variation, $`\mathrm{\Delta }|\widehat{\xi }_r|`$, serves as a measure for the probability of loop penetration into the convection zone, based on the assumption that the loop segments lifted furthest toward the superadiabatic region are the most likely to leave the overshoot region. The radial envelope determines the displacement of the perturbation wave and thus the formation of rapidly growing loops inside the convection zone is favoured around its maxima at $`\varphi _{\mathrm{max}}`$ and $`\varphi _{\mathrm{max}}+\pi `$. In the case of close binaries with orbital periods of a few days, the peak-to-peak variation of fast growing Parker-type instabilities is of the order of several percent, up to about $`20\%`$. In the case of ridge instabilities, radial displacements are strongly suppressed at some longitudes, giving rise to distinct preferred longitudes in the vicinity of $`\varphi _{\mathrm{max}}`$ and $`\varphi _{\mathrm{max}}+\pi `$. Although it can not be excluded that ridge instabilities contribute to the non-uniform spot distribution, their importance for the formation of preferred longitudes at the stellar surface is questionable since the relevant equilibrium configurations are restricted to rather small parameter domains $`(B_0,\lambda _0)`$, where instabilities have large growth times. Based on the parameter study in Sect. 4.4, we would predict that preferred longitudes appear (in binaries with solar-type components) for orbital periods $`T5\mathrm{d}`$, which is in agreement with the results of Heckert & Ordway (1995). However, in some binaries with longer periods (consisting of giant components) preferred longitudes have also been observed (Berdyugina & Tuominen 1998).

The linear stability analysis establishes the existence of preferred longitudes at the bottom of the convection zone. However, the determination of the actual longitude where an *individual* loop enters the convection zone is not possible since our procedure only allows *statistical* predictions. The tendency to constitute preferred longitudes is instead discernible for a large ensemble of similar flux tubes, supposed that the environment conditions in the overshoot region, e.g., the amount and distribution of overshooting convective motions, are uniform and unaltered for a sufficiently long time. Preferred longitudes at the bottom of the convection zone are not necessarily simply related to the observed preferred longitudes of starspots at stellar surfaces, because the evolution of rising loops through the convection zone is still subject to tidal effects222If the tube evolution would be linear throughout the convection zone, the maxima of the radial envelope would, in fact, also determine preferred longitudes of flux eruption at the surface.. The preferred longitudes of flux eruption at the stellar surface based on the flux tube model described above have to be determined by non-linear numerical calculations. Results of such calculations are presented in the subsequent paper II.

## 6 Conclusion

The equilibrium and stability properties of toroidal magnetic flux tubes in the overshoot region of an active component of a close binary system are considerably affected by tidal effects. The breaking of the axial symmetry due to the presence of the companion star leads to azimuthal variations along the tube of the equilibrium configuration and of the probability of the penetration of unstable loops to the convection zone. This demonstrates the existence of preferred longitudes *at the bottom of the convection zone*, which are well pronounced for orbital periods of a few days but strongly fade with increasing period. The orientations of the preferred longitudes depend on the equilibrium latitude, the field strength, and the dominating wave mode of the unstable eigenmode. Despite their rather small magnitude, tidal effects are nevertheless capable of influencing the dynamics of magnetic flux tubes considerably because they are enhanced by the sensitive dependence of the stability properties on the superadiabaticity of the local environment and by the resonant coupling between interacting wave modes.

###### Acknowledgements.

Volkmar Holzwarth thanks Prof. S. K. Solanki and Prof. F. Moreno-Insertis for valuable discussions and the Max-Planck-Institut für Aeronomie in Katlenburg-Lindau and the Kiepenheuer–Institut für Sonnenphysik in Freiburg/Brsg. for financial support during the accomplishment of this work.

## Appendix A Numerical treatment

The perturbations of all physical and geometrical quantities (magnetic field, flow velocity, density, position, curvature,…) arising from the small perturbations of an equilibrium flux ring are expressed in terms of the displacement vector $`\xi (\varphi )`$ (see, e.g., Ferriz-Mas & Schüssler 1993). The substitution of the corresponding expressions in the equation of motion followed by a linearision with respect to $`\xi `$ yields Eq. (14), with the coefficient matrices333In the case $`\beta 1`$.

$`_{\varphi t}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Omega }}}`$ (28)
$`_{tt}`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \frac{1}{\mathrm{\Omega }^2}}{\displaystyle \frac{1M_\alpha ^2}{M_\alpha ^2}}`$ (29)
$`_\varphi `$ $`=`$ $`\left(_1{\displaystyle \frac{g_\mathrm{n}}{g_{}}}1\right)𝒴+_1{\displaystyle \frac{g_\mathrm{b}}{g_{}}}𝒵`$ (30)
$`_t`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{\mathrm{\Omega }}}{\displaystyle \frac{1+M_\alpha ^2}{M_\alpha ^2}}𝒴`$ (31)
$`_\xi `$ $`=`$ $`_1{\displaystyle \frac{1}{g_{}}}\left(\begin{array}{ccc}g_\mathrm{t}^{}& g_\mathrm{t}& 0\\ g_\mathrm{t}& 2g_\mathrm{n}& g_\mathrm{b}\\ 0& g_\mathrm{b}& 0\end{array}\right)+_3{\displaystyle \frac{1}{g_{}^2}}\left(\begin{array}{ccc}g_\mathrm{t}g_\mathrm{t}& g_\mathrm{t}g_\mathrm{n}& g_\mathrm{t}g_\mathrm{b}\\ g_\mathrm{n}g_\mathrm{t}& g_\mathrm{n}g_\mathrm{n}& g_\mathrm{n}g_\mathrm{b}\\ g_\mathrm{b}g_\mathrm{t}& g_\mathrm{b}g_\mathrm{n}& g_\mathrm{b}g_\mathrm{b}\end{array}\right)`$ (42)
$`+_2{\displaystyle \frac{1}{g_{}}}\left(\begin{array}{ccc}g_\mathrm{t}& g_\mathrm{n}& g_\mathrm{b}\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),`$
where the constant matrices

$$=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),𝒴=\left(\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right),𝒵=\left(\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ 1& 0& 0\end{array}\right),$$
(43)
and the abbreviations

$`_1`$ $`=`$ $`{\displaystyle \frac{1}{\left(1M_\alpha ^2\right)}}{\displaystyle \frac{\rho _\mathrm{e}}{p_\mathrm{e}}}{\displaystyle \frac{1}{\gamma }}{\displaystyle \frac{g_{}}{\kappa }}`$ (44)
$`_2`$ $`=`$ $`_1^2\left[\left(1\right)M_\alpha ^2+\left(1M_\alpha ^2\right)\right]\gamma {\displaystyle \frac{g_\mathrm{t}}{g_{}}}`$ (45)
$`_3`$ $`=`$ $`_1^2{\displaystyle \frac{\gamma ^2}{2}}\mathrm{\Delta }\left(1M_\alpha ^2\right)`$ (46)
$`\mathrm{\Delta }`$ $`=`$ $`\beta \delta {\displaystyle \frac{2}{\gamma }}\left({\displaystyle \frac{1}{\gamma }}{\displaystyle \frac{1}{2}}\right)`$ (47)
have been used. The components $`g_\mathrm{x}=\left(g_{\mathrm{eff}}x\right),x[t,n,b]`$ of the effective gravitational acceleration $`g_{\mathrm{eff}}`$ are expressed in the co-moving trihedron.

Owing to the azimuthal variation of the equilibrium quantities induced by the presence of the companion star, the coefficient matrices (29) – (42) can be expressed in the form of Eq. (15) and comprise small contributions which are, in lowest order, $`\pi `$-periodic in the longitude $`\varphi `$. Following the ansatz in Eq. (17), these contributions cause the coupling of wave modes $`\widehat{\xi }_\mathrm{n}`$ with wave numbers $`n=m,m\pm 2,m\pm 4,\mathrm{}`$, which is expressed either by the 3-term recursion formula in Eq. (18) or by the Hill-type algebraic system

$$\underset{(\omega )}{\underset{}{\left(\begin{array}{ccccc}_{m2K}& 𝒞_{m2K}& _{m2K}& & \\ & \mathrm{}& \mathrm{}& \mathrm{}& O\\ & _m& 𝒞_m& _m& \\ O& & \mathrm{}& \mathrm{}& \mathrm{}\\ & & _{m+2K}& 𝒞_{m+2K}& _{m+2K}\end{array}\right)}}\left(\begin{array}{c}\widehat{\xi }_{2K}\\ \mathrm{}\\ \widehat{\xi }_0\\ \mathrm{}\\ \widehat{\xi }_{2K}\end{array}\right)=0,$$
(48)
with $`K\mathrm{}`$ and the block matrices

$`_n`$ $`=`$ $`{\displaystyle \frac{1}{2}}[(n2)(i_{\varphi ,c}+_{\varphi ,s}\omega (_{\varphi t,c}i_{\varphi t,s}))`$ (49)
$`+_{\xi ,c}i_{\xi ,s}+\omega \left(i_{t,c}+_{t,s}\right)`$
$`\omega ^2(_{tt,c}i_{tt,s})],`$
$`_n`$ $`=`$ $`{\displaystyle \frac{1}{2}}[(n+2)(i_{\varphi ,c}_{\varphi ,s}\omega (_{\varphi t,c}+i_{\varphi t,s}))`$ (50)
$`+_{\xi ,c}+i_{\xi ,s}+\omega \left(i_{t,c}_{t,s}\right)`$
$`\omega ^2(_{tt,c}+i_{tt,s})],`$
and

$`𝒞_n`$ $`=`$ $`n^2n\omega _{\varphi t,0}\omega ^2_{tt,0}+_{\xi ,0}`$ (51)
$`+i\left(n_{\varphi ,0}+\omega _{t,0}\right).`$
The eigenfrequencies $`\omega `$ are determined by the roots of the dispersion relation, $`det\left[(\omega )\right]=0`$. A treatment of the case $`K\mathrm{}`$ is not practical, but since the coupling between adjacent wave modes is small \[$`,𝒪(ϵ^3)`$\], the wave mode spectrum can be truncated at a finite number of constituents, here $`K=3`$. A verification of the eigenfrequency and the determination of the eigenfunctions $`\xi (\varphi )`$ is based on the following approach. Let the independent solutions $`x_{(i)}(\varphi ),i=1,\mathrm{},n`$ of a system

$$x^{}=𝒜(\varphi )x\text{with}𝒜(\varphi +\pi )=𝒜(\varphi )$$
(52)
represent the fundamental matrix $`𝒳(\varphi )=(x_{(1)},\mathrm{},x_{(n)})`$ with $`det[𝒳(\varphi )]0`$. Furthermore, we have $`𝒳(\varphi +2\pi )=𝒳(\varphi )`$, where $``$ is the constant *transition matrix* with $`det()0`$. Since an arbitrary solution of Eq. (52) can be represented by a constant vector $`c`$ in the form $`x(\varphi )=𝒳(\varphi )c`$, it follows that $`x(\varphi +2\pi )=𝒳(\varphi +2\pi )c=𝒳(\varphi )c`$. Periodic solutions with $`x(\varphi +2\pi )x(\varphi )`$ thus require $`c=c`$ and an eigenvalue of the transition matrix $``$ equal to unity. If the fundamental matrix fulfils the initial condition $`𝒳(0)=`$, the transition matrix is calculated by integrating over a period of $`2\pi `$, i.e., $`𝒳(2\pi )=`$. In our case, the transition matrix $``$ and its eigenvalues depend on the eigenfrequency, $`\omega `$. If an eigenfrequency determined by the dispersion relation above does not entail $`c=c`$, $`\omega `$ is iteratively improved and the eigenfunction $`\xi (\varphi )`$ determined.