Magnetar Oscillations I: strongly coupled dynamics of the crust and the core

Quasi-Periodic Oscillations (QPOs) observed during Soft Gamma Repeaters giant flares are commonly interpreted as the torsional oscillations of magnetars. The oscillatory motion is influenced by the strong interaction between the shear modes of the crust and Alfven-like modes in the core. We study the dynamics which arises through this interaction, and present several new results: (1) We show that global {\it edge modes} frequently reside near the edges of the core Alfven continuum. (2) We compute the magnetar's oscillatory motion for realistic axisymmetric magnetic field configurations and core density profiles, but with a simplified model of the elastic crust. We show that one may generically get multiple gaps in the Alfven continuum. One obtains discrete global {\it gap modes} if the crustal frequencies belong to the gaps. (3) We show that field tangling in the core enhances the role of the core discrete Alfven modes and reduces the role of the core Alfven continuum in the overall oscillatory dynamics of the magnetar. (4) We demonstrate that the system displays transient and/or drifting QPOs when parts of the spectrum of the core Alfven modes contain discrete modes which are densely and regularly spaced in frequency. (5) We show that if the neutrons are coupled into the core Alfven motion, then the post-flare crustal motion is strongly damped and has a very weak amplitude. Thus magnetar QPOs give evidence that the proton and neutron components in the core are dynamically decoupled and that at least one of them is a quantum fluid. (6) We show that it is difficult to identify the high-frequency 625 Hz QPO as being due to the physical oscillatory mode of the magnetar, if the latter's fluid core consists of the standard proton-neutron-electron mixture and is magnetised to the same extent as the crust. (Abstract abridged)


Introduction
Since the discovery of quasi periodic oscillations (QPOs) in the lightcurves of giant flares from soft gamma repeaters (SGR) (Israel et al., 2005;Strohmayer & Watts, 2005;Watts & Strohmayer, 2006;Barat et al., 1983) there has been considerable interest in their physical origin. One of the appealing explanations is that the QPOs are driven by torsional oscillations 1 of the neutron stars whose mag-presence of an ultra strong magnetic field (B ∼ 10 14 − 10 15 G) that is frozen into the neutron star and penetrates both the crust and the core. The field provides a channel for an intense hydro-magnetic interaction between the motion of the crust and the core, which becomes effective on the timescale of 1 second (Levin 2006). Since the QPOs are observed for hundreds of seconds after the flare, it is clear that the coupled motion of the crust and the core must be considered. In recent years, significant theoretical effort has gone into the study of this problem (e.g., Glampedakis et al. 2006, Levin 2007, Gruzinov 2008b, Lee 2008. This paper's analysis is based, in part, on an approach of Levin (2007, L07).
To make progress in computing the coupled crust-core motion, L07 has studied the time evolution of an axisymmetric toroidal displacement of a star with axisymmetric poloidal magnetic field. In that case the Alfven-type motions on different flux surfaces decouple from each other, a well-known fact from previous MHD studies (for a review see Goedbloed & Poedts 2004, hereafter GP). One can then formulate the full dynamics of the system in terms of discrete modes of the crust which are coupled to a continuum of Alfven modes in the core. L07 demonstrated that (i) the global modes with frequencies inside the continuum are strongly damped via a phenomenon known in MHD as resonant absorption (see GP), and (ii) in many cases, asymptotically the system tends to oscillate with the frequencies close to the continuum edges. This result was later confirmed by Gruzinov 2008b, who has used a powerful analytical technique to solve the L07's normal-mode problem (Gruzinov noted that the resonant absorption is mathematically equivalent to Landau damping). Oscillations near the continuum edge frequencies were also observed in a number of numerical general-relativistic MHD simulations of purely fluid stars (Sotani et al. 2008, Colaiuda et al. 2009, Cerda-Duran et al. 2009).
Apart from finding QPOs near the continuum edges, L07's dynamical simulations identified transient QPOs with drifting frequencies; these were transiently amplified near the crustal frequencies. No explanation for the origin of the drifts was given.
In this paper, we extend the previous analyses of the hydro-magnetic crust-core coupling in an essential way. In section 2, we re-analyse L07's toy model of a single oscillator coupled to a continuum, and we show that this system generically contains the edge normal modes with frequencies near the continuum edges. We show that these modes dominate the late-time dynamics of the system, and develop a formalism which allows one to predict analytically the edge mode's amplitude from the initial data. We then explore the effect of viscosity on the system (introduced as a friction between the neighbouring continuum oscillators), and show that the edge mode is longer lived than all other motions of the system. We also provide a non-trivial analytical formula for the time dependence of the overall energy dissipation.
In section 3, we describe how transient QPOs, not associated with the normal modes of the system, are obtained when parts of the core spectrum consist of densely and regularly spaced discrete modes (and in section 5 we show that such an array of discrete modes is expected when the magnetic field in the core is not perfectly axisymmetric but has some degree of tangling). As a by-product of our analysis, we explain the origin of the QPO frequency drifts seen in L07 simulations. We provide simple analytical fits to the drifts, and show that when the regularity of the continuum sampling is removed (e.g, when the frequencies are sampled as random numbers picked from the continuum range), the drifts disappear.
In section 4, we set up models with a more realistic hydro-magnetic structure of the neutron-star core. We show how to find the continuum modes and their coupling to the crust for an arbitrary axisymmetric poloidal field, with an arbitrary density profile on the core (L07s calculations, for simplicity and concreteness, were restricted to constant-density core and homogeneous magnetic field). We treat a more complicated case of a mixed axisymmetric toroidal-poloidal field, with radial stratification, in the Appendix B. We demonstrate that for realistic field configurations, the Alfven continuum of modes coupled to the crust may show a number of gaps. If a crustal mode frequency belongs to one of these gaps, a strong global discrete mode arises which dominates the late-time dynamics and whose frequency also belongs to the gap. The frequency of the gap global mode does not generally coincide with, but is often close to that of the crust. We suggest that it was these gap modes that appeared in Lee's (2008) calculations as well-defined discrete global modes.
So far, only axisymmetric magnetic fields have been considered in the magnetar-QPO literature, with the Alfven continuum modes occupying the flux surfaces of the field. In section 5 we argue that if the field is not axisymmetric but instead is highly tangled, then the Alfven continuum modes become localized within small regions of individual field lines, and therefore become dynamically unimportant. Instead, a set of discrete Alfven modes appears, with the spacing between the modes strongly dependent on the degree of field tangling. We devise a phenomenological prescription which allows us to parametrize the field tangling for computing the dynamically important modes, and in-troduce an easily solvable "square box" model suitable for exploring the parameter range.
Finally, in section 6, we use the suite of models built in the previous sections to explore their connection to the QPO phenomenology. We find that (a) within the standard magnetar model, it is possible to produce strong long-lived or transient QPOs with frequencies in the range of around 20-150Hz, but only if the neutrons are decoupled from the Alfven-like motion of the core; this implies that at least one of the baryonic components of the core is a quantum fluid. (b) Our models could not produce the high-frequency 625Hz QPO within the standard paradigm of a magnetar core composition.
2. An oscillator coupled to a continuum: edge modes In this section, we study the motion of a harmonic oscillator (which we hereafter call the large oscillator) which is coupled to a continuum of modes. 2 . This model was introduced in L07 and it provides a qualitative insight into the behaviour of crustal modes (represented by the large oscillator) coupled to a continuum of Alfven modes in the core of a magnetar. L07 found that if the large oscillator's proper frequency was within the range of the continuum frequencies, then the late-time behaviour of the system was dominated by oscillatory motion near the edges of the continuum interval. Here, we give an explanation of this phenomenon in terms of the edge modes. Our analysis allows us to use initial data and predict the displacement amplitudes and frequencies of the system at late times.
The model consists of the large mechanical oscillator with mass M and proper frequency ω 0 , representing a crustal elastic shear mode. Attached to the large oscillator is a set of N smaller oscillators of mass m n and proper frequency ω n constituting a quasi-continuum of frequencies ω n (where n = 1, 2, ..., N ). The continuum is achieved when N → ∞ while the total small-oscillator mass Σm n remains finite. The convenient pictorial representation is through suspended pendulae, as shown in Fig. 1 (see also Fig. 2 of L07).
The equations of motion are obtained as follows. Each small oscillator is driven by the the motion of the large 2 In many areas of physics similar models have been studied, notably in quantum optics and plasma physics. By contrast with the case studied here, in these models the range of the continuum frequencies is not limited. oscillator:ẍ where x n is the displacement of the n th small oscillator in the frame of reference of the large oscillator, x 0 is the displacement of the large oscillator in the inertial frame of reference, and the right-hand side represents the non-inertial force acting on the small oscillator due to the acceleration of the large one. The large oscillator experiences the combined pull of the small ones: Hereω 0 is the frequency of the big pendulum corrected for the mass loading by the small pendulae, i.e.ω 2 0 = ω 2 0 (M + i m i ) /M .

Time-dependent behavior.
In this subsection we explore the behavior of this system by direct numerical simulations. We found this to be helpful in the building of our intuition. We defer a semianalytical normal-mode analysis to the next subsection. We follow L07 and for concreteness concentrate on a specific example; it will be clear that the conclusions we reach are general. We choose ω 0 = 1rad/second and mass M = 1. We choose a total number of 1000 small pendulae with frequencies ω n = (0.5 + n/1000)rad/second and masses m n = m = 10 −4 , to mimic the continuum frequency range between 0.5rad/second and 1.5rad/second. The simulation is initiated by displacing the large oscillator while keeping the small pendulae relaxed (this mimics the stresses in the crust), and then releasing. The subsequent motion of the system is then followed numerically by using a second order leapfrog integration scheme which conserves the energy with high precision. The resulting motion of the large pendulum can be decomposed into three stages (see (1) During the first 50-60 seconds, there is a rapid exponential decay of the large oscillator's motion, during which most of the energy is transferred to the multitude (i.e., the 'continuum') of small oscillators. This is the so-called phenomenon of "resonant absorption", which has been studied for decades in the MHD and plasma physics community (e.g., Ionson 1978, Hollweg 1987, Goedbloed & Poedts 2004, L07, Gruzinov 2008b. In this first stage, the amplitude of the big pendulum motions drops by a factor of ∼ 100. (2) After ∼ 60 seconds, the exponential decay stops abruptly as the large oscillator now reacts to the collective pull of the small ones. This second stage is characterized by a slow algebraic decay of the amplitude of the big pendulum displacement. Gruzinov (2008b) explains this as being due to the branch cut in the oscillator's response function.
(3) The motion of the large oscillator stabilizes at a constant level (L07 missed this stage in his simulations, which he stopped too early). Fourier transform reveals 2 QPOs at the frequencies close to the continuum edges, ω = 0.5 and ω = 1.5; the same QPO frequencies can be observed in the previous stage (2) as well.
What is the origin of the QPOs, and how is this eventual stability established? In Fig. 4 and 5, we show how the amplitude of the small oscillators evolves with time.
After the initial resonant absorption phase, the amplitude is distributed as a Lorentzian centered on the frequency around ω = 1; this is because the small oscillators in resonance with the large one are the ones which gain the most energy. However, in subsequent times we see that the energy exchange occurs between the small oscillators 3 , and that the net result of this exchange is the energy flow towards the oscillators whose frequencies are near the edges. By the time the third stage begins, the amplitudes of the oscillators near the edge stabilize and their phases become locked. They are pulling and pushing the large oscillator in unison. In the next subsection, we show that this behavior is due to the presence of the edge normal modes, and we shall derive their frequencies and amplitudes.

Finding eigenmodes
In this section we deal with the system of coupled harmonic oscillators, and one should be able to find its normal modes using the standard techniques (Landau and Lifshitz mechanics,§23). However, the fact that all small oscillators are attached to the large one, and there is no direct coupling between the small oscillators, allows us a significant shortcut (in Appendix A, we treat a more general problem of several large oscillators coupled to a multitude of the core modes). We proceed as follows: Suppose that we impose on the large oscillator a periodic motion with angular frequency Ω, by driving it externally with the force F ext = F 0 (Ω) exp(iΩt). This motion in turn drives the small oscillators according to Eq. (1): which has the steady state solution: where we have omitted the time dependent factor exp(iΩt) on both sides. The combined force f cont of the small oscillators acting back on the large one (see Eq. (2)) is given by According to Newton's second law, If Ω corresponds to the normal-mode frequency, then F 0 (Ω) = 0. Hence by substituting Eq. (5) into Eq. (6) we get the following eigenvalue equation for Ω: In the continuum limit N → ∞, the above equation becomes where ρ(ω) = dm/dω is the mass per unit frequency of the continuum modes.
In the discrete case, the solutions of Eq. (7) are N − 1 frequencies Ω i that are within the quasi-continuum (ω i < Ω i+1 < ω i+1 , for i = 1, 2, ...N − 1; 'quasi-continuum modes'), and 2 modes with frequencies Ω low and Ω high (the 'edge-modes') that are near the edges, but outside, of the continuum (i.e. Ω low < ω 1 and Ω high > ω N ). It can be shown from Eq. (7) that in the limit N 1 and m n M , the contribution of the small oscillator to the ith quasi-continuum mode is completely dominated by the pendulae that are in close resonance with the mode. More precisely, one can show that as the number of oscillators N increases and m n decreases, the number of small oscillators contributing to the mode energy remains constant. However, for the two edge modes there is no such singular behavior in the limit of large N , and consequently they play a special role in the dynamics of the system. This last point is clearly seen in the continuum case represented by Eq. (8). The eigenvalue equation has no real solutions in the range of small-oscillator continuum ω min < Ω < ω max , since the response function G(Ω) is ill-defined in this interval 4 . However, the edge modes on both sides of the continuum interval remain, and their frequencies can be found by numerically evaluating the zero points of G(Ω) in Eq. (8). For the numerical of the previous subsection, one finds Ω low = 0.5 − 8.2· 10 −6 and Ω high = 1.5 + 8.6· 10 −4 . Analytically, one can find the following scaling for the distance δω edge between the mode frequency and the nearest edge ω edge of the continuum range: where C is a constant of order unity. The larger is the density of continuum modes at the edge ρ(ω edge ), the further is the edge mode pushed away from the continuum range. It is particularly interesting to consider the case when the continuum interval is limited by a turning point (L07) with the divergent density of states near the edge, ρ(ω) = A/ |ω − ω edge |, where A is a constant. In this case the distance from the edge-mode frequency to the nearest continuum edge is given by The quadratic dependence in Eq. (10) vs. the exponential dependence in Eq. (9) implies that the continua with 4 There is a complex solution if the integration in the expression for G(Ω) is performed along the contour chosen according to the Landau rule. One then obtains a "resonantly absorbed" or "Landau-damped" mode (Gruzinov 2008b, L07), which exactly represents the exponential decay of stage (1) in our numerical experiment of the previous subsection.
turning points typically feature much more pronounced edge modes and stronger QPOs than the ones with linear edges. In the next section, we show how to calculate the edge-mode amplitudes and QPO strengths from the initial data.

Late time behavior of the system
In the continuum limit, the only modes with real oscillatory frequency are the edge modes. Thus, as we demonstrate explicitly below, they dominate the late-time dynamics of the system when the number N of small oscillators becomes large. Our analysis proceeds as follows: Lets define a new set of variables, expressed as a vector X with components for n = 1, ..., N . With these new variables, the kinetic energy of the system is a trivial quadratic expression where the inner product of two vectors A and B is de- The potential energy is a positive-definite quadratic form, whose exact form is unimportant here. The mutually orthogonal eigenmodes X i can be found via a procedure outlined in the previous section 5 . Their eigenfreguencies Ω i are identified by finding zeros of G(Ω) in Eq. (7), and the corresponding eigenvector components are given by Lets assume that we initiate our simulation by displacing the large oscillator by an amount x 0 (0) while keeping the small oscillators relaxed x n (0) = 0 and all initial velocities at zero. In the new variables, the initial state of the system is given by the vector X(0), where X 0 = √ M x 0 (0) and X n = √ m n x 0 (0). The time evolution of the system is given by: Substituting the initial conditions, and the expression in Eq. (12) for the eigenvector components, we get The coordinate of the large oscillator is simply given by x 0 (t) = X 0 (t)/ √ M . For the continuum of small modes, the above expansion breaks down, since the eigenvalue equation has no real solutions inside the continuum range. However, the edge modes are well defined, and they determine the dynamics at late times. Therefore, for the continuum case we can still write down the analogous expression which is valid only at late times: The sums of Eq. (14) are replaced with the corresponding integrals, and we have the following expression for the displacement of the large oscillator at late times: (16) This expression is in excellent agreement with the numerical simulations. In the numerical example of subsection 2.1, the upper edge mode dominates the late-time behavior of the system, and its calculated contribution is plotted in Fig. 3, together with the numerically simulated motion.

The effect of viscosity
We now add an extra degree of realism by introducing viscous friction into the system. In MHD, continuum modes are spatially localized, and the effect of viscosity is to frictionally couple the neighboring modes (see, e.g., Hollweg 1987). In our simple model we introduce viscosity by adding frictional forces between the neighboring oscillators, where f n,n+1 is the force from the n'th oscillator acting on the (n + 1)'th. We now calculate how the system dissipates energy as a function of time. We will show that it occurs in two stages (see Fig. 6): (1) Initially, the small oscillators are strongly and simultaneously excited by the "Landau-damped" large oscillator, then they become dephased, with the average relative motion between the neighboring oscillators growing linearly in time. This leads to a very rapid dissipation of the bulk of the initial energy. (2) The edge modes persist, since the participating small oscillators move in phase and the energy dissipation is small. The energy of the modes is damped exponentially on a timescale much longer than that of the first stage. The dissipated energy is given by In the continuum limit, the small oscillators are labeled not by a discrete index n, but by a continuous variable λ. The expression for the dissipated energy is then whereγ is the viscous coefficient. After the initial exponential damping of the large oscillator and the excitation of the small oscillators, the latter initially move independently, with wherex(λ) is the amplitude of the λ'th oscillator. From the above equation, we then obtain where the ... stands for time-averaging over many oscillation periods. For times t d log x λ /dω λ the second term on the right-hand side of Eq. (21) dominates. For a simple model with dω λ /dλ = const and ρ(ω) = const, where E is the total energy of the system and A = (γ/ρ)(dω λ /dλ). The analytical solution for the energy and dissipated power, agrees very well with numerical simulations, see Fig. 6. While the equations above were derived for restrictive assumptions (dω λ /dλ = const and ρ(ω) = const), we found that the analytical formulae in Eq. (23) provide a good fit for a large variety of simulations. This is because it is the small oscillators with the frequencies near that of the large oscillator which carry most of the energy, and in that narrow band our approximations hold. After the energy dissipation due to dephasing is over, only the edge modes remain. This is illustrated in Fig. 7, where we show how the energies of the small oscillators evolve with time. At late times, only the oscillators taking part in the edge modes move substantially; this is because they remain in phase and do not dissipate much. At this stage the energy is drained by slow exponential decay of the edge modes.

Transient and drifting QPOs
Finite-size MHD systems feature a mix of continuum and discrete modes (see Poedts et al., 1985 and GP). For axisymmetric field configurations the continuum modes occupy the whole flux surfaces and play an important role in the oscillatory dynamics; this was the motivation for L07 and our study of the previous section. We will argue in section 5 that if the core field is highly tangled, the continuum modes become localized in space and discrete core modes will play a more important role. Thus it is important to study the case when the the crustal modes are coupled to a set of discrete core modes. In this section we show that if the frequencies of the discrete modes are regularly spaced in some frequency intervals, then the system displays transient QPOs that are entirely missed by its normal-mode analysis. This is interesting from the observational point of view, since many of the observed magnetar QPO features are transient.
Suppose that a set of discrete modes are located in the interval ∆ω around frequency ω 0 and are separated by a regular frequency interval δω, and assume the following hierarchy: δω After the modes are excited, they are initially in phase but will de-phase rapidly on the timescale 1/∆ω. However, at times t n = 2πn/δω the modes come into phase again and pull coherently on the large oscillator. Therefore, a transient QPO feature should appear around these times at a frequency close to ω 0 . In Fig. 8 and Fig. 9 we show the dynamical spectrum from a simulation where the conditions for transient QPOs were engineered for two frequencies.
The transient QPOs agree well with the expectations. As is seen from the figures, the strongest transient QPOs are those whose frequencies are the closest to that of the large oscillator; this is because the response of the large oscillator is the strongest around its proper frequency.
One can now easily understand the frequency drifts in Fig. 10 of L07 ( Fig. 10 in this paper). In the simulations of that paper, the core continuum was sampled with a set of densely and regularly-placed Alfven modes by slicing the field into finite-width flux shells. The spacing δω between the modes was not constant but a function of the Alfven frequency ω. In that case, the QPO drifts with the QPO frequency ω(t) given by the inverse relation With this relation we are able to fit all of L07 drifting QPOs, as shown in Fig. 10 and 11. Note that multiple QPOs correspond to different branches of the Alfven continuum. As expected, the drifting QPOs amplified near the crustal frequencies, since there the response of the crust to the core modes' pull is the strongest.

More realistic magnetar models
In this section we extend the constant magnetic field and constant-density magnetar model from L07 to include more realistic pressure and density profiles and more general (but still axisymmetric) magnetic field configurations. Our aim is to use this model to: 1) calculate numerically Alfven eigenmodes and their eigenfrequencies on different flux surfaces inside the star, in order to determine the continuous spectrum of the fluid core, and 2) use these modes to simulate the dynamics of a realistic magnetar. In order to calculate the Alfven eigenmodes and eigenfrequencies for a realistic magnetar model, we employ the linearized equations of motion for an axisymmetric magnetized, self-gravitating plasma. The general equations, which are derived in detail in Poedts et al. (1985, hereafter P85) and given in their equations (53) and (54), constitute a fourth order system of coupled ordinary differential equations in the case of a mixed poloidal and toroidal magnetic field. The formalism of P85 is briefly summarized in the Appendix B. In the case of a purely poloidal magnetic field, the system simplifies to two uncoupled second order differential equations (P85, equations (70) and (71)).

The model
We assume our star is non-rotating and neglect its deformation due to the magnetic pressure, which is expected to be small. Therefore, we consider a spherically symmetric background model that is a solution of the Tolman-Oppenheimer-Volkoff equation (TOV equation). The hydrostatic equilibrium is calculated using a SLy equation of state (Douchin & Haensel 2001;Haensel & Potekhin, 2004;Haensel, Potekhin & Yakovlev 2007), which can be found in tabulated form on the website http://www.ioffe.ru/astro/NSG/NSEOS/. The integration of the TOV equation is performed using a 4th order Runge-Kutta scheme, integrating from the center of the star outward until we reach a mass density ρ = 1.3· 10 14 g/cm 3 , which is consistent with the crust-core interface in the equation of state from Douchin & Haensel (2001).
The resulting model has a central mass density ρ c = 10 15 g/cm 3 , a total mass of 1.40 M and a radius of R core = 1.07· 10 6 cm. To this spherical model we add a poloidal magnetic field, which we generate by placing a circular current loop of radius a and current I around the center of the star. The field is singular near the current loop, however all the field lines which connect to the crust (and thus are physically related to observable oscillations) carry finite field values. This particular field configuration is cho-sen as an example; there is an infinite number of ways to generate poloidal field configurations. In the appendix B we will add to this field a toroidal component and calculate the corresponding Alfven continuum of the core.

The continuum
In order to find the equations of motion for the magnetized material in the neutron star core, we would need to add self-gravity to the ideal magnetohydrodynamic equations. This problem has been solved by P85 in a tour the force mathematical approach. In that paper the authors assume a self-gravitating axisymmetric equilibrium with a field geometry consisting of mixed poloidal and toroidal field components, and they derive linearized equations of motion. For this field geometry it is convenient to work with so-called flux-coordinates (ψ, χ, φ). The basic concept behind this curvilinear coordinate system is the magnetic flux-surface, which is defined as the surface perpendicular to the Lorentz force F L ∝ j × B. From this definition it is clear that the magnetic field lines lie in flux surfaces. If one considers a closed loop on a flux surface which makes one revolution around the axis of symmetry, then the magnetic flux ψ through the loop depends on the flux surface only and is the same for all of the loops. Therefore ψ is chosen as the coordinate labeling the flux surfaces. In each flux-surface we can denote a position by its azimuthal angle φ and its 'poloidal' coordinate χ, which is defined as the length along φ = const line. In P85, it is shown that the equations of motion allow for a class of oscillatory solutions that are located on singular flux surfaces, constituting a continuum of eigenmodes and eigenfrequencies. In the case of a purely poloidal field (B = B χ ), the continuum solutions are degenerate and polarized either parallel (ξ χ ) or perpendicular (ξ φ ) to the magnetic field lines. In the latter case the displacement is φ-independent. It is clear that in contrast to the χ-polarized modes, the φ-polarized modes are purely horizontal and are therefore unaffected by gravity. This latter case is considered here.
The equation of motion is in this case simply the Alfven wave equation: where the operator F is given by Here x is the distance to the magnetic axis of symmetry. Although in the presence of a mixed poloidal and toroidal field the equations still give rise to a continuous set of solutions, the calculations are significantly complicated as the continuum modes are affected by the toroidal component of the field, by gravity, and by compressibility. For the sake of simplicity we will ignore toroidal fields in our dynamic simulations. We will however, calculate the continuum frequencies for a mixed poloidal and toroidal field in the Appendix B.
For determining the spectrum of the core continuum, the appropriate boundary conditions are ξ φ (χ = χ c ) = 0, where χ c (φ) marks the location of the crust-core interface. With this boundary condition, Equation (26) constitutes a Sturm-Liouville problem on each separate flux surface ψ. Using the stellar structure model and magnetic field configuration from section 4.1, we can calculate the eigenfunctions and eigenfrequencies for each flux surface ψ. The reflection symmetry of the stellar model and the magnetic field with respect to the equatorial plane assures that the eigenfunctions of Eq. (26) are either symmetric or anti-symmetric with respect to the equatorial plane. We can therefore determine the eigenfunctions by integrating Eq. (26) along the magnetic field lines from the equatorial plane χ = 0 to the crust-core interface χ = χ c (ψ). Let us consider the odd modes here for which ξ φ (0) = 0, and solve Eq. (26) with the boundary condition ξ φ (χ c ) = 0 at the crust-core interface; for even modes, the boundary condition is dξ φ (0) /dχ = 0. We find the eigenfunctions by means of a shooting method; using fourth order Runge-Kutta integration we integrate from χ = 0 to χ = χ c . The correct eigenvalues σ n and eigenfunctions ξ n (χ) are found by changing the value of σ until the boundary condition at ξ n is satisfied. In this way we gradually increase the value of σ until the desired number of harmonics is obtained. In figure 12 we show a typical resulting core-continuum.
According to Sturm-Liouville theory the normalized eigenfunctions ξ n of Eq. (26) form an orthonormal basis with respect to the following inner product: Where δ m,n is the Kronecker delta and r = 4πρ/B χ is the weight function. We have checked that the solutions we find satisfy the orthogonality relations.
We are now ready to compute the coupled crust-core motion. Here we follow L07 and assume that the crust is an infinitely thin elastic shell 6 . We label the lattitudinal 6 It is straightforward to relax this assumption, and carry out the analysis of this section for the finite crustal thickness. However, from Section 2 it is clear that the interesting dynamics is dominated by Fig. 12.-The red curves show the Alfven frequencies σn as a function of the angle θ(ψ), the polar angle at which the fluxsurface ψ intersects the crust. Since we are only considering odd crustal modes, the only Alfven modes that couple to the motion of the star are the ones with an odd harmonic number n. This particular continuum was calculated using a poloidal field with an average surface value B surface ∼ 6· 10 14 G, generated by a circular ring current of radius a = R * /2. location by the flux surface ψ intersecting the crust, and consider the crustal axisymmetric displacementsξ φ (ψ). In the MHD approximation, the magnetic stresses enforce a no-slip boundary condition at the crust-core interface, such that ξ φ (ψ, χ c ) =ξ φ (ψ, χ c ) instead of ξ φ (ψ, χ c ) = 0. It is useful to make the following substitution where we choose the function w (ψ, χ) so that (a) it corresponds to the static displacement in the core and hence satisfies F (w (ψ, χ)) = 0, and (2) w (ψ, χ c ) = 1. Therefore the new quantity satisfies the boundary condition ζ (ψ, χ c ) = 0 and can be expanded into the Alfven normal modes ξ n which satisfy the same boundary conditions. We now proceed by substituting Eq. (29) into Eq. (26) thus obtaining a simple equation of motion for ζ From the definition of the operator F it follows that for the odd modes Here the constant K (ψ) is chosen such that w (ψ, χ c ) = 1, in order that ζ = 0 on both boundaries. We expand ζ and w into a series of ξ n 's: ζ (ψ, χ, t) = n a n (ψ, t) ξ n (ψ, χ) w (ψ, χ) = n b n (ψ) ξ n (ψ, χ) Eq. (30) reduces to the following equations of motion for the eigenmode amplitudes a n ∂ 2 a n (ψ) ∂t 2 + σ 2 n (ψ) a n (ψ) = −b n (ψ) These equations show how the core Alfven modes are driven by the motion of the crust. To close the system, we must address the motion of the crust driven by the hydromagnetic pull from the core. The equation of motion for the crust is given by the spectral structure of the core Alfven waves; therefore in order to flesh out the physics we choose the simplified model of the crust.
Where the acceleration due to elastic stresses L el is where θ is the polar angle (cf. L07). The acceleration L B due to the magnetic stresses between the crust and the core can be expressed as where x is the distance to the axis of the star, Σ is column mass-density of the crust and α is the angle between the magnetic field line and the normal vector of the crust. It is convenient to express the crustal displacementξ φ as a Fourier series, being a sum normal modes of the free-crust problem. Using Eq. (36) is straightforward to show analytically that the eigenfunctions f l of the free-crust problem (Eq. (35) with eigenfrequencies Here Y l0 is the m = 0 spherical harmonic of degree l. The normalized functions f l form an orthonormal basis, so that where δ l,m is again the Kronecker delta. The crustal displacement can then be expressed in terms of f l Substituting Eq. (41) into Eq. (35) we obtain the equations of motion for the crustal mode amplitudes c l We can express L B as χ=χc Up to this point the derived equations of motion for the crust and the fluid core are exact. We are now ready to discretize the continuum by converting the integral of Eq. (42) into a sum over N points θ i . In order to avoid the effect of phase coherence (see section 3) which caused drifts in the results from L07, we sample the continuum randomly over the θ-interval [0, π/2]. In the following, functional dependence of the coordinate ψ or θ (ψ) is substituted by the discrete index i which denotes the i-th flux surface.
These are the equations that fully describe dynamics of our magnetar model. As with the toy model from section 2 we integrate them using a second order leap-frog scheme which conserves the total energy to high precision. As a test we keep track of the total energy of the system during the simulations. Further we have checked our results by integrating equations (45) and (45) with the fourth-order Runge-Kutta scheme and found good agreement with leapfrog integration.

Results
Based on our section-2 results, we have a good idea of what type of dynamical behavior should occur in our more realistic magnetar model. First, we expect that crustal modes with frequencies inside the Alfven continuum will be damped quickly by resonant absorption ("Landaudamping" in the terminology of Gruzinov 2008b). Second, as with our previous model we expect the late time behavior of the system to show QPO's near the edges of the continuum, or edge modes. The third important feature of our model is that the continuum may possibly contain gaps, as is shown in Fig. 13. In this case there is the possibility that crustal frequencies fall inside the gaps and remain undamped. In all of our simulations these expectations have come true. We will now show the results from a simulation which illustrate the above mentioned effects.
The basic freedom of choice that we have for our model is the strength and geometry of the equilibrium magnetic field. We choose here a purely poloidal magnetic field with an average strength at the surface of B surf = 10 15 G, induced by a circular current loop of radius a = 0.5R * . This field gives us a gap in the continuum at frequencies 53 < ω < 78 Hz.
We consider the lowest degree odd crustal modes with frequencies ω 2 = 40 Hz and ω 4 = 84.5 Hz, which we couple to 5000 continuum oscillators (the Alfven continuum). We sample the continuum at 1000 randomly chosen flux surfaces, and at each flux surface we consider 5 Alfven overtones.
Like with our toy model from section 2, we initiate the simulation by displacing the crust (c 2 = c 4 = 1) while keeping the continuum oscillators (the Alfven modes) relaxed (a in = 0).
In Figures 14 and 15 we show the resulting power spectrum for 2 different models. In the first one, the crustal frequencies are located inside the core continuum range, and the peaks due to the edge modes appear. By contrast, in the second case one of the crustal frequencies belongs to the gap, and a peak representing the global gap mode strongly stands our above the background. We note that the gap-mode's frequency lies close to but does not coincide with the crustal-mode frequency; we found this to be a generic feature of our models, with the difference of 10% for the typical model parameters. The gap modes are particularly interesting because they have relatively large amplitudes, and are not as strongly damped by viscosity as the edge modes.
It must be emphasized that for all persistent modes in the system, the position in the frequency space of the core Alfven continuum plays the key role in setting the globalmode frequency and in determining its longevity.
We note that Lee (2008) has used a different method to identify discrete modes in a magnetar with similar magnetic configuration to ours. These modes were not associated with crustal frequencies, and we strongly suspect that they were located in the gaps of the continuum spectrum and could be identified with the edge or gap modes presented in this work.

Tangled magnetic fields
Our preceding discussion of the continuum was predicated on the foliation of the axisymmetric magnetic field into the flux surfaces, with each of the singular continuum mode localized on the flux surfaces. These modes are "large"-they are coherent over the spacial extent compara-  netar with a single 'gap' in the Alfven continuum. The global mode within the gap is not damped, and its frequency is similar, but not identical, to that of the crustal mode in the same gap. ble to the size of the system, and thus they play an important role in the overall dynamics-they are responsible for the resonant absorption of the crust oscillations, and contribute to generating the edge and gap modes. But what happens if the field cannot be foliated into the flux surfaces, but is instead tangled in a complicated way? One can argue that the continuum part of the spectrum still persists, as follows: Consider an arbitrary field line anchored at the crustcore interface at both ends, and choose a tube of field lines of infinitesimal radius which is centered on the original field line (see Fig. 16). It is clear that a twisting Alfven mode exists in the tube: it is obtained by the circular rotation of the fluid around the central field line, propagating along the central field line with the local Alfven velocity. Since there is a continuum of the field-line lengths, there is also a continuum of Alfven modes. However, the modes we constructed are highly localized in space and and have a small leverage in the overall dynamics. We conjecture that the more tangled the fields are, the less role do the singular continuum modes play in the overall dynamics.
Whilst we cannot rigorously prove this conjecture, we can motivate it as follows: consider an area element δS of random orientation with the normaln inside the star, and consider a shearing motion along the element. This shearing motion will be resisted by the Bn component of the magnetic field, with the effective shear modulus of order where ... stands for averaging over the area element. For ordered field, it is possible to choose the orientation of the area element so that µ eff 0; the presence of such an orientation makes a fundamental difference between MHD and elasticity theory and is responsible for the presence of continuous spectrum in MHD. However, if the linear size of the δS is greater than the characteristic length on which the field is tangled, then µ eff is non-zero for all orientations ofn. Therefore, for highly-tangled fields there can be no large-scale singular continuum modes, and their existence is restricted to the small scales. Hence our assertion that for strongly tangled fields continuum modes play a secondary dynamical role.
One is then faced with the problem when crustal modes are coupled to a set of discrete core Alfven modes. In Appendix A we show how to find the eigen-solution of such a problem, provided that all of the coupling coefficients are known.
How does one quantify the degree to which the fields are tangled? Some insight comes from the numerical simulations of Braithwaite and colleagues, who have studied what type of fossil fields remain in a stratified star after an initial period of fast relaxation. Consider a stable fossil field field configuration, such as the one obtained in the simulations of Braithwaite and Spruit (2004) and Braithwaite and Nordlund (2006) [see also Gruzinov (2008a) for analytical considerations]. There, the final field is nearly, but not perfectly axisymmetric and has a small-scale random component. For a less-centrally concentrated initial field, Braithwaite (2008) shows that the final fossil field is in general non-axisymmetric and can have a complicated topology. 7 As a starting point, we shall consider the nearly axisymmetric field with a small random component. The latter acts like a small extra shear modulus µ eff and dynamically couples the flux surfaces of the axisymmetric component. We then quantify the degree of tangling by the relative value of µ eff and B 2 /(4π).

simple model: "square" neutron star
To study this idea further, we specify a very simple model of a neutron star, motivated by the one considered in Levin (2006, hereafter L06) see Fig. 17 that never-theless captures the essential physics.
Consider a perfectly conducting homogeneous fluid of density ρ contained in a box with width L x , length L y and depth L z . The magnetic field in this box is everywhere aligned with the y-axis and its strength is a function of x only. We assume that gravity is zero and consider a Lagrangian displacement ξ (x, y, t) of the fluid along the zdirection; we specify periodic boundary conditions in this direction (one should think of the z direction as azimuthal). We now add to this model a small effective shear modulus µ eff due to the field tangling. The fluid equation of motion is: Here c A (x) is the Alfven velocity and c s = µ eff /ρ is the µ eff -generated shear velocity. If we assume a small shear speed, i.e. c s << c A , Eq. (47) reduces to We now find the core Alfven eigenmodes. After adapting the no-slip boundary conditions the problem can be easily solved by separation of variables ξ (x, y, t) ∝ e iωt sin {πm[(y/L y ) + 1/2]} X (x), where m = 1, 2, .... Plugging this in Eq. (48) we find for the the x−dependent part of the solution: Here ω A,m (x) = πmc A (x) /L y can be interpreted as the frequency of the m-th Alfven overtone at x. From the above expression it is clear that in the limit of very small c s , the solution for X must be close to zero everywhere except in a very small neighborhood of ω A,m (x) = ω. It is in this limit that the solutions are located on singular flux surfaces. However, in the presence of the non-vanishing shear velocity c s , the eigenmodes spread out on neighboring field lines, effectively coupling the motion on different flux surfaces. The continuum of Alfven frequencies ω A,m (x) will in this case be no longer solutions of the system. Instead, the coupling term gives rise to a discrete set of solutions rather than a continuum. Eq. (50) is the mathematical equivalent of Schrödingers equation, which can in general cases be solved numerically. However, for many special case exact solutions exist. Let us consider, for the sake of simplicity, a field configuration in our box such that: We can rewrite Eq. (50) as follows: This differential equation is the mathematical equivalent of the quantum harmonic oscillator problem for which the exact solution is well known. The eigenfrequencies are given by ω 2 mn = π (1 + 2n) mc s √ a c A /L y + c 2 A,0 π 2 m 2 /L y . (53) Here n (= 0, 1, ...) is the 'quantum' number labeling the harmonic-oscillator wavefunctions. We see that instead of a continuum, we obtain a densely packed discrete set of frequencies with the frequency spacing ω m,n − ω m,n−1 ∼ πmc s √ a c A /L y ω m,n .
With the no-slip boundary conditions on the left and right sides x = ±L x /2, the eigenvalue equation must be solved numerically. An example of such calculation is shown in Fig. 18. There, the spacing between the discrete Alfven modes is shown to increase as one increases the effect of the field tangling characterized by the µ eff .
We now introduce the crustal modes into the problem by making the top and bottom of the box elastic and mobile. We allow their displacementξ t,b (x, t) in the z-direction, and impose the boundary conditions on the sides: Here the subscripts "t" and "b" stand for the top and bottom of the box, respectively. The top and bottom are assumed to be thin and have mass M cr and surface density σ = M cr /(L x L z ) . The crustal equation of motion is given by where v s is the shear velocity in the crust. The crustal angular frequencies are given by ω cr j = jπv s /L y with the corresponding crustal mode eigenfunctionsξ j = sin{jπ[(x/L x ) + 1/2]} , where j = 1, 2, ... is roughly equivalent to l in the spherical case. The symmetry of the problem allows one to consider either symmetricξ t =ξ b or antisymmetricξ t = −ξ b crustal modes. This will couple to the symmetric (m = 1, 3, 5, ...) or antisymmetric (m = 2, 4, 6, ...) Alfven modes of the core.
Just as in section 4, it is now convenient to define a new variable ζ(x, y, t) for the core displacement: where The new variable observes the regular boundary condition ζ = 0 on all the box edge, and satisfies the following inhomogeneous partial differential equation: where The advantage of the new variable is that it satisfies the regular boundary condition ζ = 0 on all the boundaries of the box. It can therefore be expanded as a series consisting of eigenfunctions ξ mn of the right-hand side of Equation (48): ζ(x, y, t) = Σ mn a mn (t)ξ mn (x, y).
The rest of the procedure is very similar to that in section 4. We expand the crustal displacement into a series consisting of the eigenmode wavefunctionsξ j : where p j (t) and q j (t) are real numbers. The magnetar deformation is now fully represented by a set of generalized coordinates [p j (t), q j (t), a mn (t)]. The coupled equations of motion are derived by following the procedure specified in section 4. We obtain the following system of equations: and β j(mn) = ∂ξmn(x,y) ∂y y=Ly/2ξ Thus we have obtained a system of linear second-order differential equations, which describes the time evolution of the square-box magnetar. These equations are solved by trancating all the series [i.e., rescricting the range of indices (m, n, j)] and then by either solving the eigenvalue problem in order to find the normal modes, or by integrating the equations numerically 8 . One then checks that the series trancation does not introduce errors in the magnetar's motion within the frequency range of our interest. So far we have worked in the approximation of the thin crust, i.e. we have effectively included the crustal modes which have no radial nodes in their wavefunction. However, several observed high-frequency QPOs, and in particular the strong QPO at 625Hz (Watts & Strohmayer 2006) necessitate introduction of higher radial-order modes into our model. In the square-box model we do this phenomenologically, as follows. We assume that higher radialorder crustal modes have amplitudes p sj (t) and q sj (t) and natural eigenfrequencies ω cr sj , with s being the number of radial nodes, and assume that they cause displacement at the crust-core interface given byξ j (x). This mirrors realistic spherically-symmetric case where the functional form of the crust-core displacement due to the torsional ∇×Y lm mode of the n'th radial order is a very weak function n. The amplitudes p sj (t) and q sj (t) are then introduced on into the equations of motions (62) and (63) in the same way as the other p j and q j amplitudes, with the same jdependent coupling coefficients but with ω cr sj instead of ω cr j on the left-hand side of Eq. (63). We now have the basic ingredients of building a phenomenological modes with tangled fields. To sum up, we (1) quantify tangling using the effective shear modulus, (2) Find discrete core eigenmodes and evaluate their coupling to the crustal model, and (3) Either find eigenfrequencies of the total star by diagonalizing the potential energy of the system, or simulate the time-dependent behavior directly.
An example of a resulting powerspectrum is shown in Fig 19 for the model described in this section. Fig. 19.-Power spectrum for the dynamics of a magnetized box as described in the text. In this particular model we have used the maximum possible shear modulus, corresponding to a maximally tangled field. The Alfven motion in the box is coupled to 9 of the lowest frequency 'crustal' modes, plus a high frequency crust mode at 630 Hz.

What do our models tell us about magnetar
QPOs?
In this paper we have developed a formalism which allows one to build a magnetar model with a variety of the spectral features of the core Alfven waves, including continua with gaps and edges, and the large-scale discrete modes generated by the field tangling. We have constructed a number of magnetar models and explored the resulting QPOS, both for the case of axisymmetric magnetar with core Alfven continuum, and for the "square" magnetar models with the tangled fields (see the previous section). The full range of model parameters, and detailed comparison with the data will be the subject of a separate study. For now, we have restricted ourselves to the standard magnetar model, in which the core is a perfect conductor, the field of ∼ 10 15 G penetrates both the core and the crust, and the proton fraction in the star is the one tabulated in Haensel, Potekin, and Yakovlev (2007). Our models give us the following robust conclusions, as compared against QPO observations: (1). A number of strong QPOs have been observed in the 1998 and 2004 giant flares, with frequencies ranging between 18 Hz and 150 Hz ( Israel et al. 2005, Strohmayer and Watts 2006, Watts and Strohmayer 2006. These QPOs can be qualitatively explained as gap and/or edge modes of sections 4 and 2, or even transient QPOs of section 3 9 . However, this was only possible if the neutrons were decoupled from the Alfven waves in the core. If the neutrons took part in the Alfven motion, then the effective mass of the Alfven modes shifted up by a factor of 20 − 40 and their frequencies shifted down by a factor 4 − 8 (Easson & Pethick 1979, Alpar et al. 1984, van Hoven & Levin 2008, Andersson et al. 2009). As a result, all modes at frequencies above ∼ 50Hz were strongly damped (see Fig.  20). Increasing the magnetic-field tension by a factor of 3 did not affect this conclusion (Fig 21). For the spherical magnetar models of section 4 we obtain similar results if couple the neutrons to the Alfven motion in the core. The key point that we would like the reader to appreciate is that Alfven modes in the core are key to determining the frequency and strength of the observable QPOs, and thus QPOs are very sensitive probe of the core interior.
(2) A number of the high-frequency QPOs have been measured in the 2004 giant flare by Watts and Strohmayer (2006), the strongest among them being the QPO at 625 Hz. This QPO is particularly strong and long-lived in the hard x-rays, reaching the amplitude of ∼ 25% over the time interval of ∼ 100 seconds (i.e., it persists for almost 10 5 oscillation periods!). Watts and Strohmayer (2006) argued that this frequency corresponds to the crustal shear mode with a single radial node (see also Piro 2005); this interpretation, if correct, would strongly constrain the thickness of the crust and rule out the fluid strange stars as magnetar candidates (Watts & Reddy, 2007). To investigate this suggestion, we have introduced several high-frequency lowj crustal modes into our square-box simulations. However, as is demonstrated in Figs. 14 and 15, the high-frequency modes are strongly damped and at no time during the simulations do we observe any significant power at those frequencies. This is to be expected. No natural axisymmetric model has gaps in the Alfven continuum at such high frequencies, so global modes are strongly absorbed. One could expect that if the Alfven modes are discrete in the core due to the field tangling, this problem would not arise. However, even in the discrete case the frequency spacing between the modes is around 20 Hz, which is much smaller than 600 Hz. Thus the grid of Alfven waves is so dense that it is effectively seen as the absorbing continuum by the modes around 600 Hz. Our detailed simulations, of the type shown in Figs. 14 and 15, fully confirm this qualitative picture.
The concern about the viability of high-frequency QPOs as being due to the physical oscillations of standard-model magnetars has been raised in the original L06 paper on the basis of rather simplistic calculations. As our work here shows, more detailed calculations partially alleviate the L06 concern for the frequencies below ∼ 150Hz, but only if the neutrons are decoupled from the Alfven motion in the core, i.e. if at least one baryonic superfluid (protons or neutrons) are present in the neutron-star core. Our analysis sustains L06 concern for the high-frequency QPOs, in particular for the strong long-lived QPO at 625Hz. Its explanation seems to require either QPO production in the magnetosphere, or a somewhat radical revision of the magnetar model. Just how radical this revision has to be will be explored in a separate study.
Our work presented here has several shortcomings. We have limited ourselves to the linear approximation, and a non-linear regime may bring surprises. Direct non-linear simulations of axisymmetric oscillations of a magnetised fluid star has been carried out recently by Cerda-Duran, Stergioulas, & Font (2009). At this stage it is difficult to say whether non-linearities introduce significantly new QPO features to their model; their results have largely been in agreement with the linear simulations of Colaiuda, Beyer, & Kokkotas (2009). However, the computational techniques seem promising and we do not exclude that large-amplitude simulations of stars with the crust will show qualitatively new features. Another limitation of our work is that we have assumed that once the flare sets the magnetar into motion, the magnetar's oscillations are not driven externally. This may not be the case in real flares: some energy stored in the pre-flare magnetar may be released gradually, and this release could be extended in time into the flare's tail 10 . The latter consideration is straightforward to incorporate phenomenologically into our model, and we plan to address it in our future work.

A. Multimodal crust-core system
In this Appendix we generalize the normal-mode treatment of Section 2.2, and write down the general prescription of how to find the eigenmodes when several "large" crustal shear modes are coupled to a multitude of small core Alfven modes, provided the coupling coefficients are known. In this paper, the coupling coefficients are worked out in simple models of sections 4 and 5; we shall postpone the discussion of how the coefficients are computed in more general case to the future paper.
Let us denote the displacement of the crustal and core modes by X n and x i respectively. Since both the crustal and the core modes are not directly coupled to themselves (i.e., X's are only coupled to x's), most general equations of motion take the formẌ n + Ω 2 X n = Σ i α ni x i (A1) where Ω n and ω j are the proper frequencies of the crustal and core modes, and α's and β's are the coupling coefficients. We look for an oscillatory solutions of the above equations with angular frequency Ω. One can trivially re-write these equation as a matrix eigen-equation with Ω 2 as an eigenvalue, and solve it using standard methods. However, if the number of crustal modes is not too large, it is convenient to make a shortcut. Using the second of Eq. (A1) to express x i 's through X n 's, and substituting into the first one, we get the following equation: where the elements of the matrix G are given by One obtains the eigenfrequencies by finding numerically the zeros of det G mn .
B. Core continua with a mixed axisymmetric toroidal-poloidal magnetic field In this appendix we will calculate the continuum of Alfven frequencies in a magnetar core in the case of a axisymmetric magnetic field with mixed toroidal and poloidal components. The general MHD equations of motion for spherically symmetric, self-gravitating equilibrium with an axisymmetric field, are derived in detail in P85. In contrast to the special case of a purely poloidal field (see section 4.2) which leads to two uncoupled differential equations, the continuum for a mixed toroidal-poloidal field is described by a system of fourth order coupled ODEs. Due to this coupling, the solutions are complicated as they are no longer polarized in the directions parallel (so-called "cusp solutions") and perpendicular (Alfven solutions) to the magnetic field lines, but rather have a mixed character. Strictly speaking, one can only speak of an "Alfven continuum" in the limit that the variations in ρ, P and B 2 are small in the magnetic flux-surfaces. The general equations of motion are given in Eqs. (53) and (54) of P85. We note however, that in magnetars the speed of sound c >> c A , and therefore we consider Poedts et al.'s equations (53) and (54) in the incompressible limit (P85, Eqs. (73) and (74)). For completeness we give the equations here, The variables Y ≡ i B 2 φ ξ χ − B φ B χ ξ φ /B χ B 2 and Z ≡ i (B χ ξ χ + B φ ξ φ ) /B 2 are components of the fluid displacement perpendicular and parallel to the magnetic field lines, the operator F ≡ i∂/∂χ is a differential operator along ψ intersects the crust. In the presence of a toroidal field, the degeneracy between the cusp-solutions and the Alfven solutions is broken and we find two separate solutions for each wave number n; waves with primarily Alfven character (red curves) and waves with primarily cusp character (blue curves). This particular continuum was calculated using a poloidal field with an average surface value B p,surface ∼ 6· 10 14 G (again generated by a circular ring current of radius a = R * /2) and a toroidal field strength at the equator and the crust-core interface of Bt,0 = 3· 10 14 G (see Eq. (B3)). the field lines, N χ ≡ − (1/B χ ρ) (∂ρ/∂χ) (∂P/∂χ) can be thought of as a Brunt-Väisälä frequency for displacements along the field lines. According to Gauss' law for magnetism, the toroidal component of the magnetic field is of the form B φ = f (ψ)/ , where is the distance to the polar axis, and f (ψ) is an arbitrary function of ψ. In the following calculation we adopt a toroidal field component of the form Here θ(ψ) is the polar angle at which the flux surface ψ intersects the stellar crust. Clearly this choice for B φ is completely arbitrary and one could in principle try many different toroidal geometries.
As with our calculation of the Alfven continuum in the case of a purely poloidal field (section 4.2), we adopt the zerodisplacement boundary conditions at the crust, and use the fact that our equilibrium model is (point-) symmetric with respect to the equatorial plane. This enforces the existence of classes of symmetric and anti-symmetric eigenfunctions, Y n (χ) and Z n (χ). We consider only the odd modes, and calculate the eigenfunctions by means of the shooting method; we use a fourth order Runge-Kutta scheme to integrate Eqs. (B1) and (B2). Starting with Y (0) = 0 and Z(0) = 0 at the equator, we integrate outward until we reach the crust at χ = χ c . We find the eigenfrequencies by changing the value of σ until we match the boundary conditions at the crust. A resulting continuum is plotted in Figure 22.