text:
# Resonant origin for density fluctuations deep within the Sun: helioseismology and magneto-gravity waves

## 1 Introduction

Helioseismology has become a precision tool for studying the inner workings of the Sun, providing one of the only direct probes of physical properties as a function of depth, below the convective layer and into the radiative zone. The frequencies of thousands of normal modes have been measured, and current levels of agreement between these measurements and calculations preclude deviations of density profiles and sound speeds (as functions of depth) in excess of one percent from solar-model predictions. Careful analysis of small discrepancies between models and measurements which have arisen in the past has helped refine the models, by winnowing out small errors in opacities and other quantities used as inputs.

The vast majority of helioseismic analyses to date are performed in the approximation which neglects the sun’s magnetic field. This approximation is generally very reasonable since the energy density, $`𝐁^2/8\pi `$, of the expected fields is much smaller than other energies in the problem, such as gas pressures, deep in the solar interior. This is expected to be particularly true in the solar radiative zone, below the turbulent convection which is believed to be responsible for the solar dynamo.

Very little is directly known about magnetic field strengths within the radiative zone – see  for early studies. A generally-applicable bound is due to Chandrasekhar, and states that the magnetic field energy must be less than the gravitational binding energy: $`B^2/8\pi GM_{}^2/R_{}^4`$, or $`B10^8`$ G. Stronger bounds are possible if one makes assumptions about the nature and origins of the solar magnetic field. For instance, if it is a relic of the primordial field of the collapsing gas cloud from which the sun formed , then it has been argued that central fields cannot exceed around 30 G . Similarly, the dynamo mechanism can only generate a global field in the radiative zone of a newly-born Sun with amplitude below 1 G . Even stronger limits, $`B10^6`$ G apply  if the solar core is rapidly rotating, as is sometimes proposed. On the other hand, it has recently been argued  that fields up to 7 MGauss could persist in the radiative zone for billions of years and are consistent with current observational bounds. Other authors  have recently entertained radiative-zone fields as large as 30 MGauss. Since the initial origin and current nature of the central magnetic field is unclear, we consider as admissible any magnetic field smaller than of order 10 MGauss.

The purpose of the paper you are now reading is to examine the influence of magnetic fields in more detail, deep within the solar radiative zone. Our motivation for so doing is threefold. First, although magnetic fields in the radiative zone are believed to be small, they cannot be directly measured, and so quantitative bounds on their size require the explicit calculation of their effects. Ours is a first step along these lines, where we assume a particularly simple geometry which we argue to be appropriate when the background magnetic field is orthogonal to the local density gradients. Such configuration would occur near the equatorial plane in a dipole magnetic field. Second, plasmas in magnetic fields are notoriously unstable, with small perturbations often giving rise to strong effects, perhaps amplifying the implications of a small magnetic field. Indeed our analysis uncovers such instabilities, which give rise to surprisingly large effects. Third, the recent advent of neutrino telescopes opens a new window on the solar interior, motivating the exploration of the implications of solar magnetic fields for solar-neutrino properties.

In a nutshell, our main findings are these:

* We find, contrary to naive expectations, that profiles of helioseismic $`g`$-modes as a function of solar depth can be appreciably altered in the radiative zone even for magnetic fields which are not unreasonably large . In particular, density profiles due to these waves tend to form comparatively narrow spikes at specific radii within the sun, corresponding to radii where the frequencies of magnetic Alfvén modes cross those of buoyancy-driven ($`g`$-type) helioseismic modes. Due to this resonance, energy initially in $`g`$-modes (presumably due to turbulence at the bottom of the convective zone) is directly pumped into the Alfvén waves, causing an amplification of the density profiles in the vicinity of the resonant radius. This amplification continues until it is balanced by dissipation, resulting in an unexpectedly large density excursion at the resonant radii.
* Since the only such level crossing which occurs within the radiative zone occurs with $`g`$-modes, and not $`p`$-modes, the latter are less substantially affected by radiative-zone magnetic fields. This makes it very unlikely that radiative-zone magnetic fields can significantly alter standard analyses of helioseismic data, which ignore solar magnetic fields.
* We identify the implications of magnetic fields for the amplitudes and frequencies of solar $`g`$-modes, in principle permitting inferences about radiative-zone fields to be drawn once these modes are experimentally observed. We find that the bound implied by the consistency of our analysis with existing discussions of $`p`$-mode helioseismic waves is at present uninterestingly weak.

The main message in these findings is that the detailed shape of helioseismic $`g`$-waves can be significantly changed by reasonable radiative-zone magnetic fields, and this sensitivity to magnetic fields is driven by the occurrence of level crossing between helioseismic $`g`$-modes and Alfvén waves. Of course, although this cartoon describes the Alfvén and helioseismic waves as if they were separate, a real calculation should simply diagonalize the normal modes of the linearized hydrodynamic problem including the magnetic fields. We perform such a diagonalization explicitly, and find the magneto-gravity (MG) waveforms whose peaks correspond to the above resonance argument. This diagonalization allows us to determine how the mode energies and excitation rates depend on the strength of the background magnetic field.

The resonance we identify has also been recognized as possibly playing a role elsewhere in the Sun. In particular, the feeding of energy into this resonance has been proposed  as a mechanism for heating the solar corona. Indeed the analysis of these authors is very similar to the one we present here, with the main difference being their neglect of gravitational forces (as is appropriate for applications to the solar atmosphere, but not for the radiative zone).

Indeed, we believe the study of the implications for neutrinos of the Alfvén resonant layers is well worth pursuing, particularly since density fluctuations themselves can have effects, without invoking the existence of large magnetic moments. Two other lines of argument also encourage such a study. First, previous experience along these lines  argues that neutrinos are most sensitive to density fluctuations near the depths at which MSW resonant oscillations take place, $`r0.3R_{}`$, where the resonances we derive indeed can arise. Second, the density spikes we here find are much narrower than are generic helioseismic $`g`$-wave profiles, possibly permitting Alfvén waves to cause effects which garden-variety $`g`$-waves cannot . We present a discussion of neutrino propagation through these waves in Ref. , where we find potentially interesting prospects for detectable effects if radiative-zone magnetic fields lie in the range 10 – 100 kG.

We organize our presentation as follows. In section 2 we set up the hydrodynamic problem in the presence of gravity and magnetic fields, and give the total set of MG equations whose solution is required. We next derive the linearized MHD equations, which we solve by relating all physical quantities to the radial component of the perturbation magnetic field amplitude. In section 3, we apply reasonable boundary conditions to set up an eigenvalue problem which determines the dispersion relation for the discrete eigenfrequency spectrum, $`\omega _n(B_0)`$. For a simple geometry we identify the labels, $`n`$, and spectrum of the MG eigenwaves. In section 4 we focus on the resonance layer position and compute how its position and width vary from mode to mode. We also here present results for the size of the maximum density excursions at the resonant points as a function of magnetic field strength. Finally, section 5 discusses some observational constraints on the wave frequencies and amplitudes.

## 2 The master equation for MG waves

In this section we briefly review how magnetic fields influence the equations of hydrostatic equilibrium on which helioseimic analyses are based. To this end we first set up the relevant MHD equations, and linearize them to identify the differential operator whose spectrum dictates the allowed MG wave frequencies.

Two equations of fluid dynamics are not changed by the presence of a magnetic field – at least within the regime of present interest. The first of these is the continuity equation, expressing conservation of mass:

$$\frac{d\rho }{dt}+\rho u=0,$$
(1)
where $`d/dt=/t+𝐯`$ is the usual convective derivative, with $`𝐯`$ representing the fluid velocity, and $`\rho `$ and $`p`$ are the fluid’s mass density and pressure. The variable $`u=\mathrm{div}𝐯`$ need not vanish if the fluid is compressible 111This variable $`u`$ was the main focus of our earlier asymptotic analysis of MG waves .

The second unchanged equation expresses energy conservation for a polytropic medium. This equation states:

$$\frac{dp}{dt}\gamma \frac{p}{\rho }\frac{d\rho }{dt}=(\gamma 1)Q,$$
(2)
where $`\gamma =c_p/c_V`$ is given by the ratio of heat capacities, and $`Q`$ is the sum of all energy density sources and losses. In principle, electromagnetic processes enter into Eq. (2) through their contributions to $`Q`$, but for ideal MHD we neglect both the heat conductivity and viscosity contributions to energy losses, as well as the ohmic dissipation, $`Q=j^2/\sigma _{cond}`$, where $`\sigma _{cond}`$ is the electrical conductivity.

The magnetic fields enter more directly into the Euler equation (conservation of momentum), which is of the form

$$\rho \frac{d𝐯}{dt}=p+\rho 𝐠+\frac{1}{4\pi }\left[\mathrm{rot}𝐁\times 𝐁\right],$$
(3)
with the local force of gravity given by $`\rho 𝐠`$. The local acceleration due to gravity is related to the Newtonian potential by $`𝐠=\varphi `$, with $`\varphi `$ given by the Poisson equation $`^2\varphi =4\pi G\rho `$. As usual, $`G`$ here denotes Newton’s constant. The last term of equation (3) expresses the contribution of the Lorentz force to the local momentum budget.

Finally, the system is completed by Faraday’s equation,

$$\frac{𝐁}{t}=\mathrm{rot}[𝐯\times 𝐁]+\frac{c^2}{4\pi \sigma _{cond}}^2𝐁,$$
(4)
that governs the time evolution of the magnetic field. Here $`\sigma _{cond}`$ represents the fluid’s conductivity. In what follows we take the plasma to be an ideal conductor, meaning we take $`\sigma _{cond}`$ to be large enough to neglect the last term in this last equation. Clearly these equations reduce to those of standard Helioseismology Models (SHM) in the limit $`B0`$.

### 2.1 Linearized form

Our interest is in small-amplitude oscillations about a background configuration, and so at this point we split all variables into background and fluctuating quantities, $`A=A_0+A^{}`$, with the $`A^{}`$ denoting small fluctuations about the background value $`A_0`$. Since we take our background configuration to be static, $`𝐯_0=0`$, we regard $`𝐯`$ to be pure fluctuation, but suppress the prime on $`𝐯`$ in the interests of notational simplicity.

The system of linearized equations obtained in this way is

$`{\displaystyle \frac{\rho ^{}}{t}}+(𝐯)\rho _0+\rho _0u=0,`$
$`\rho _0{\displaystyle \frac{𝐯}{t}}+p^{}𝐠\rho ^{}{\displaystyle \frac{1}{4\pi }}\left[\mathrm{rot}𝐁^{}\times 𝐁_0\right]=0,`$
$`{\displaystyle \frac{p^{}}{t}}+(𝐯)p_0+\gamma p_0u=0,`$
$`{\displaystyle \frac{𝐁^{}}{t}}=\mathrm{rot}\left(𝐯\times 𝐁_0\right),`$ (5)
where, as described above, $`𝐯`$, $`\rho ^{}`$, $`p^{}`$ and $`𝐁^{}`$ are the small Euler perturbations.

In addition to our previously-mentioned assumptions of ideal conductance – i.e. $`\sigma _{cond}\mathrm{}`$ – and static backgrounds – $`𝐯_0=0`$, with background quantities time independent – these equations make three further assumptions. First, they adopt the Cowling approximation, which amounts to the neglect of perturbations of the gravitational potential, (i.e.: $`\varphi ^{}=0`$). Second, they assume the fluctuations are adiabatic, $`Q^{}=0`$, as is satisfied with good accuracy by high frequency oscillations. Finally, they assume the background field to satisfy $`\mathrm{rot}𝐁_0=0`$.

The backgrounds $`\rho _0,p_0`$ etc., must themselves satisfy the hydrodynamic equations, which for negligible (or constant) magnetic fields reduce to the usual equations of hydrostatic equilibrium, $`p_0=\rho _0𝐠`$. Finally, we suppose there to be no energy sources or sinks in the background, $`Q_0=0`$. Under these circumstances the background density $`\rho _0`$ is as given by Standard Solar Models (SSM) .

If we also assume $`\gamma `$ to be constant – as is a good approximation within the whole radiative zone – we can derive from the system of linearized equations (5) an evolution equation for the velocity field $`𝐯`$:

$$\frac{^2𝐯}{t^2}=(𝐠𝐯)+(\gamma 1)𝐠u+c_s^2uv_A^2\frac{}{t}\left[\frac{𝐁_\mathrm{𝟎}}{B_0}\times \mathrm{rot}𝐛\right],$$
(6)
where $`𝐛=𝐁^{}/B_0`$ and $`v_A=B_0/\left(4\pi \rho _0\right)^{1/2}`$ defines the Alfvén speed.

### 2.2 Geometry

To proceed further we next choose the geometry of the solution we shall seek. Our interest is in a magnetic field which is perpendicular to the local density gradients and gravitational fields. Our interest is also in the deep interior of the radiative zone, but not right down to the solar center.

For these purposes we can take $`𝐁_0`$ to be both constant and approximately uniform. It then suffices to work within an approximately rectangular, rather than cylindrical, geometry. In these circumstances it is convenient to choose a Cartesian coordinate system whose z-axis is the ‘radial’ direction, as defined by (but opposite in direction to) the local gravitational acceleration, $`𝐠`$. With this choice we take $`z(0,R_{})`$, where $`z=0`$ represents the solar center and $`z=R_{}`$ denotes the solar surface. We are actually mainly interested in the radiative zone, for which $`z0.7R_{}`$.

Since we have in mind magnetic fields perpendicular to physical gradients and gravitational fields, we take $`B_0`$ to lie along the x-axis and the gradients to lie along the z-axis, leading to the components $`𝐁_0=(B_0,0,0)`$, and $`𝐠=(0,0,g(z))`$.

### 2.3 Isolating variables

Returning to our original problem, we wish to solve the linear MHD equations without resorting to the WKB approximation. Our first goal to this end is to eliminate all variables in terms of a single one, which we choose to be the $`z`$-component of the perturbed magnetic field, $`b_z(z)`$. The spectrum of MHD modes is then found from the eigenvalues of the linear ordinary differential equation satisfied by $`b_z`$.

This is most easily done by transforming $`𝐯`$ and $`𝐛`$ to Fourier space in the $`x`$ and $`y`$ directions, such as in $`𝐯=𝐯(z)\mathrm{exp}\left(i(\omega tk_xxk_yy)\right)`$. In this way we find that the velocity perturbations are given by the system

$`v_x`$ $`=`$ $`{\displaystyle \frac{ik_x(c_s^2ugv_z)}{\omega ^2}},`$ (7)
$`v_y`$ $`=`$ $`{\displaystyle \frac{1}{\omega ^2}}\left[ik_y(c_s^2ugv_z)i\omega v_A^2\mathrm{rot}_z𝐛\right],`$ (8)
$`v_z`$ $`=`$ $`{\displaystyle \frac{1}{\omega ^2}}[{\displaystyle \frac{}{z}}(gv_z)+(\gamma 1)guc_s^2{\displaystyle \frac{u}{z}}+`$ (10)
$`+i\omega v_A^2\mathrm{rot}_y𝐛].`$
The equation for the magnetic field perturbation $`𝐛(z)`$ is then derived from the last equation of eqs. (5),

$$i\omega 𝐛=ik_x𝐯\frac{𝐁_0}{B_0}u,$$
(11)
and by the definition of the compressibility, $`u`$:

$$u=ik_xv_x+ik_yv_y+\frac{v_z}{z}.$$
(12)
In our later applications we are most interested in the density perturbation given in terms of the other quantities by

$$i\omega \frac{\rho ^{}}{\rho _0}=\frac{\omega }{k_xH}b_z+u.$$
(13)
Substituting the expressions for $`v_{x,y}`$ from eqs. (7) and (12), and finding the magnetic field component $`b_y`$ from the Maxwell equation $`\mathrm{div}𝐛=0`$ we find expressions for the dynamical variables in terms of $`b_z`$. In so doing we also use two further approximations. We first assume the adiabatic parameter, $`\gamma `$, to be constant, and we also use the low-frequency approximation, $`\omega ^2k_x^2c_s^2`$. (This latter assumption has the effect of filtering out the acoustic $`p`$-modes from our analysis.) We find in this way:

$`b_x(z)={\displaystyle \frac{ik_x}{k_{}^2}}{\displaystyle \frac{db_z(z)}{dz}},`$
$`b_y(z)={\displaystyle \frac{ik_y}{k_{}^2}}{\displaystyle \frac{db_z(z)}{dz}},`$
$`v_z(z)={\displaystyle \frac{\omega }{k_x}}b_z(z),`$
$`v_x(z)={\displaystyle \frac{ik_x\omega }{k_{}^2}}\left[{\displaystyle \frac{gb_z(z)}{k_xc_s^2}}+\left(1+{\displaystyle \frac{\omega ^2}{k_{}^2c_s^2}}\right){\displaystyle \frac{1}{k_x}}{\displaystyle \frac{db_z(z)}{dz}}\right],`$
$`v_y(z)={\displaystyle \frac{ik_y\omega }{k_{}^2}}\left[{\displaystyle \frac{gb_z(z)}{k_xc_s^2}}+\left(1+{\displaystyle \frac{\omega ^2}{k_{}^2c_s^2}}\right){\displaystyle \frac{1}{k_x}}{\displaystyle \frac{db_z(z)}{dz}}\right],`$
$`{\displaystyle \frac{p^{}(z)}{p_0}}=i{\displaystyle \frac{\gamma \omega ^2}{k_{}^2c_s^2k_x}}{\displaystyle \frac{db_z(z)}{dz}},`$
$`{\displaystyle \frac{\rho ^{}(z)}{\rho _0}}={\displaystyle \frac{(+i)}{k_xH}}\left({\displaystyle \frac{\gamma 1}{\gamma }}b_z(z)+{\displaystyle \frac{\omega ^2H}{k_{}^2c_s^2}}{\displaystyle \frac{db_z(z)}{dz}}\right),`$ (14)
where $`k_{}^2=k_x^2+k_y^2`$. Using these, expression (12) for $`u`$ becomes

$$u(z)=\frac{\omega ^3}{k_xc_s^2}\left(\frac{g}{\omega ^2}b_z(z)\frac{1}{k_{}^2}\frac{db_z(z)}{dz}\right)0.$$
(15)
Note that we retain here in the transverse velocity components, $`v_{x,y}`$, terms which are subdominant in the quantity $`\omega ^2/k_{}^2c_s^21`$, since these are needed to obey the compressibility condition, Eq. (12) or (15).

## 3 The eigenvalue problem

We may now identify the equation satisfied by $`b_z(z)`$ itself, whose solution determines all of the other variables through the expressions derived above. We find $`b_z`$ to be determined as the solution to the following second order linear ordinary differential equation,

$`\left(1{\displaystyle \frac{k_x^2v_A^2}{\omega ^2}}\right){\displaystyle \frac{d^2b_z(z)}{dz^2}}{\displaystyle \frac{N^2}{g}}{\displaystyle \frac{db_z(z)}{dz}}+`$
$`+k_{}^2\left({\displaystyle \frac{k_x^2v_A^2}{\omega ^2}}1+{\displaystyle \frac{N^2}{\omega ^2}}\right)b_z(z)=0,`$ (16)
where $`N`$ denotes the Brunt-Väisälä frequency, as defined by

$$N^2(z)=g(z)\left(\frac{1}{\gamma p_0}\frac{dp_0}{dz}\frac{1}{\rho _0}\frac{d\rho _0}{dz}\right).$$
(17)
Equation (16) is one of our main results, and the remainder of the paper is devoted to the construction of its solutions.

### 3.1 Qualitative Analysis

Eq. (16) describes the interaction of adiabatic gravity modes with magnetic fields (magneto-gravity modes), in the approximation of low frequencies, $`\omega ^2k_x^2c_s^2`$. Before pressing on, it is instructive to examine the qualitative properties of the solutions to Eq. (16), since these capture the results we will obtain from a more detailed analysis.

In the limit $`B_00`$ (and so $`v_A0`$), Eq. (16) reduces to the standard evolution equation for ‘pure’ helioseismic $`g`$-modes, which is usually expressed in terms of the variable $`v_z`$ (which is related to $`b_z`$ through the relation $`v_z=\omega b_z/k_x`$). Since more detailed analyses – see Fig. 1 – show that $`N`$ rises from zero at the solar center, remains approximately constant $`NN_0`$, through the radiative zone, and then falls to zero again at the bottom of the convective zone, these $`g`$-modes can be thought of as the eigen-modes of oscillations inside the cavity in between the two regions where $`N`$ goes through zero. For a given wave frequency, $`\omega <N_0`$, the turning points of this cavity are given by the condition $`\omega N`$. For smaller $`\omega `$ the lower turning point gets closer to the solar center, $`z0`$, and the upper one gets slightly closer to the bottom of the convective zone (CZ).

Conversely, if gravity is turned off ($`N0`$), then Eq. (16) describes Alfvén waves, which oscillate with frequency $`\omega =k_xv_A`$ and propagate along the magnetic field lines. Notice that since $`v_A\rho ^{1/2}`$ this frequency grows with $`z`$, since the density of the medium falls.

Keeping both magnetic and gravitational fields introduces qualitatively new behavior, as may be seen mathematically because equation (16) acquires a new singular point which occurs when the coefficient of the second-derivative term vanishes. Since this singularity appears at the Alfvén frequency,

$$\omega =k_xv_A$$
(18)
it can be interpreted as being due to a resonance between the $`g`$-modes and Alfvén waves. This resonance turns out to occur at a particular radius because the Alfvén frequency varies with radius while the $`g`$-mode frequency is independent of radius (see Fig. 1). The resonance occurs where the growing Alfvén mode frequency crosses the frequency of one of the $`g`$-modes, and the resulting waveforms would be expected to vary strongly at these points.

The resonance can occur inside the radiative zone if the Alfvén frequency climbs high enough to cross a $`g`$-mode frequency before reaching the top of the radiative zone. Since $`v_AB_0`$, whether this occurs or not depends on the central field value, $`B_0`$. Our main task is to find the dependence of the spectrum of eigen-modes on magnetic field, and to elucidate when such singular resonances can appear within the Sun. Of particular interest for observational purposes is to know how the density and sound speed vary near the resonances, especially when the resonance region occurs near the locations of resonant neutrino oscillations.

The presence of a resonance at a particular radius, $`z=z_r`$, shifts the position of the upper turning points of the corresponding magneto-gravity wave towards deeper regions of the Sun. This leads to several new effects. First, the shortening of the cavity due to the presence of magnetic field causes the eigen-frequencies to depend on the magnetic field value. Second, there can be energy transfer between $`g`$-modes and Alfvén waves, within the narrow singular resonance layer, leading to the corresponding eigen-frequencies acquiring imaginary parts. Third, the nearer the upper boundary of MHD cavity is to the solar center, the stronger the $`g`$-modes are confined to the solar core.

More details of these waveforms can be obtained from a WKB-type analysis of the master Eq. (16). Near the solar center where $`N^20`$, if $`(k_xv_A)^2/\omega ^21`$ one obtains the exponential solution $`b_ze^{\pm k_{}z}`$. (Which combination of these solutions appears is fixed by boundary conditions, e.g. if $`b_z(0)=0`$ at the solar center we have $`b_z\mathrm{sinh}(k_{}z)`$.) One finds similar behaviour for the solution above the singular resonant layer, $`z>z_r`$, up to the top of the radiative zone. This exponential growth happens because of the exponential decrease of density (and so exponential increase in $`v_A`$) with $`z`$. The requirement for complex frequencies arises from the demand that the solution remain regular in the narrow Alfvén resonance layer.

### 3.2 Background Configurations

In order to acquire analytic solutions to eq. (16), we must specify a background profile for $`\rho _0(z)`$ and $`g(z)`$. We model $`\rho _0(z)`$ as an exponential density profile, $`\rho _0(z)=\rho _ce^{z/H}`$ with $`H`$ constant, and take $`g`$ to be approximately constant. We now argue these to be reasonable approximations well inside the radiative zone, and not too close to the solar center.

The variation of background quantities may be characterized by the scale heights, $`H_A^1:=|(1/A_0)(dA_0/dz)|`$, for $`A=\rho ,p,T`$ etc.. Solar models  show that the density scale height is the shortest in the radiative zone, $`H_g,H_T>H_\rho `$, and that $`H_\rho `$ takes the constant value $`H_\rho =H=0.1R_{}`$ to a good approximation, except near the solar center and solar surface. On these grounds we have that the gravitational equation, $`dg/dz=4\pi G\rho _0`$, implies $`g=4\pi G\rho _cH\left(1e^{z/H}\right)`$, which is approximately constant for $`z`$ much larger than $`H`$. For constant $`g`$ the hydrostatic equation, $`dp_0/dz=\rho _0g`$, then implies $`p_0=p_c\mathrm{exp}(z/H)`$ which also implies constancy of the temperature for gases obeying the ideal gas equation state, $`p_0/\rho _0=RT_0/\mu `$. Here $`R`$ is the gas constant, $`R=k_B/m_u`$, $`\mu `$ is the molecular weight measured in atomic mass units $`m_um_p`$, and $`k_B`$ is Boltzmann’s constant. In terms of these quantities the density scale height becomes: $`H=c^2/g`$ (= constant) , with $`c=\sqrt{RT_0/\mu }`$ being the isothermal sound velocity, which is related to the adiabatic sound velocity, $`c_s`$, by $`c_s^2=\gamma c^2`$. All of these quantities are constants within the approximations we are using.

Notice that our approximations that $`\gamma `$, $`H=0.1R_{}`$ and $`g`$ are all constant also imply approximate constancy for the Brunt-Väisälä frequency, $`N`$, within the radiative zone.

### 3.3 Solution Using Hypergeometric Functions

Changing variables to $`\xi =k_x^2v_A^2/\omega ^2=\xi _ce^{z/H}`$ we recast Eq. (16) as

$`(1\xi )\xi ^2{\displaystyle \frac{d^2b_z(\xi )}{d\xi ^2}}+\xi \left(1\xi {\displaystyle \frac{\gamma 1}{\gamma }}\right){\displaystyle \frac{db_z(\xi )}{d\xi }}+`$
$`+K_{}^2(\kappa ^2+\xi 1)b_z(\xi )=0,`$ (19)
where we introduce the dimensionless parameters $`K_{}^2=k_{}^2H^2`$, $`\kappa ^2=N^2/\omega ^2`$.

Since the equation (19) has three regular singular points – at $`\xi =0,1`$ and $`\mathrm{}`$ – it can be put into standard hypergeometric form,

$$\xi (1\xi )\frac{d^2Y}{d\xi ^2}+[c(a+b+1)\xi ]\frac{dY}{d\xi }abY=0,$$
(20)
through a change of the dependent variable $`b_z=\xi ^\sigma Y(\xi )`$. This permits an explicit solution for $`b_z`$ in terms of Gaussian hypergeometric functions, $`F(a,b;c;\xi )={}_{2}{}^{}F_{1}^{}(a,b;c;\xi )`$. Here the hypergeometric coefficients, $`a,b`$ and $`c`$, are given by

$$a=\sigma +K_{},b=\sigma K_{},c=2\sigma +\gamma ^1,$$
(21)
where the complex index, $`\sigma `$, is defined by the expression

$$2\sigma =\frac{\gamma 1}{\gamma }q,$$
(22)
with

$$q1c=\left[\left(\frac{\gamma 1}{\gamma }\right)^24K_{}^2(\kappa ^21)\right]^{1/2}.$$

General Solutions
Two linearly-independent solutions to the Hypergeometric equation, (20), are

$$Y_1=F(a,b;c;\xi ),Y_2=\xi ^{1c}F(ac+1,bc+1;2c;\xi ).$$

These solutions must be analytic for all $`\xi 0,1`$ and $`\mathrm{}`$, but they may acquire branch points at these three regular singular points depending on the values of $`a`$, $`b`$ and $`c`$. In the present case because $`Re(cab)=\gamma ^1>0`$, the solutions are bounded at $`\xi =1`$ as long as neither $`c`$ nor $`2c`$ is zero or a negative integer. Moreover, the second fundamental solution, $`Y_2`$, is finite at $`\xi =0`$ if $`Re(1c)=Req>0`$.

We shall find that the waves of interest have complex frequencies, which we parameterize as $`\omega =\omega _1(1+id)`$ with $`\omega _1`$ and $`d`$ real, with $`|d|1`$. Complex $`\omega `$ implies, from (22), that the parameters $`q`$ (and so also $`c`$) are also complex, $`q=q_1+iq_2`$, with real and imaginary parts:

$$q=q_1+iq_2\sqrt{w_1+iw_2}=$$

$$\frac{1}{\sqrt{2}}\left[\sqrt{\sqrt{w_1^2+w_2^2}+w_1}+i\mathrm{sgn}(d)\sqrt{\sqrt{w_1^2+w_2^2}w_1}\right],$$
(23)
where $`w_{1,2}`$ are given by:

$$w_1=\left(\frac{\gamma 1}{\gamma }\right)^2+4K_{}^2\left[1\frac{N^2(1d^2)}{\omega _1^2(1+d^2)^2}\right],$$

$$w_2=\frac{8dN^2K_{}^2}{\omega _1^2(1+d^2)^2}.$$

It is clear from these equations that one always has $`Req=q_1>0`$.

The general solution for $`b_z(z)`$ is therefore given by the linear combination of the independent solutions to Eq. (20):

$`b_z(z)=D_1\xi ^\sigma F(a,b;c;\xi )+`$
$`+D_2\xi ^{1\sigma 1/\gamma }F(1+ac,1+bc;2c;\xi ),`$ (24)
where one combination of the integration constants, $`D_1/D_2`$ is obtained by imposing boundary conditions at $`\xi _c=\xi (z=0)`$ and $`\xi _{}=\xi (z=z_{})`$, where $`z_{}0.7R_{}`$ denotes the top of the radiative zone.

Resonance Conditions
Although the solutions are bounded at the singular point $`\xi =1`$, they are generally not analytic there. The position $`z_r`$ determined by the condition

$$\text{Re}\xi (z_r)=(k_xv_A(z_r)/\omega _1)^2=1$$
(25)
therefore defines the position of the $`g`$-mode/Alfvén resonance. This resonance occurs within the radiative zone only if $`\xi _c<1`$ and $`\xi _{}>1`$. In this section we ask what this requires for the properties of the wave.

Expressed in terms of the central quantities, the Alfvén resonance condition is $`Re\xi (z_r)=1`$, or $`Re\xi _c=(k_xv_{Ac}/\omega _1)^2=e^{z_r/H}`$. Since $`0z_rz_{}0.7R_{}`$ this condition limits the values of $`K_{}`$ for which resonances can occur. Using the numerical values $`H=7\times 10^9`$ cm, $`v_{Ac}=(B_0/43.4\mathrm{G})`$ cm s-1 (obtained using the central density $`\rho _c=150`$ g cm-3), $`N=10^3`$ s-1 (obtained using the mean acceleration $`g=0.2`$ km s-2 and the ideal-gas adiabatic parameter $`\gamma =5/3`$), we find the following limits on the transverse wave number $`K_{}K_x=k_xH`$:

$$1\frac{7\times 10^6(\omega _1/N)}{K_x(B_0/43.4\mathrm{G})}30,$$
(26)
where the lower (upper) limit corresponds to the choice $`z_r=0`$ ($`z_r=z_{}0.7R_{}`$).

For frequencies, $`\omega _1N`$, and moderate central magnetic fields, $`B_0100G`$, we see that the existence of a resonance requires the wave number $`K_{}`$ to be large, $`K_{}1`$. Smaller $`K_{}`$ necessarily requires either larger magnetic fields or extremely low frequencies or both. Since we cannot, within the scope of our approximations, reduce the angular frequency lower than the solar angular rotation frequency (for which the Solar rotation period is $`T_{}27`$ days), a reasonable lower limit for $`\omega _1`$ might be $`\omega _110^5`$ s-1 (corresponding to a period $`2\pi /\omega _18`$ days). Even for such small frequencies one finds from Eq. (26) the condition $`3\times 10^4K_{}10^3`$, as required.

Boundary Conditions
At the solar center ($`z=0`$ and so $`\xi =\xi _c`$) we use the boundary condition which would have arisen as a smoothness condition in cylindrical coordinates:

$$b_z(0)=k_xv_z(0)/\omega =0.$$
(27)
This requires the integration constants to satisfy

$$\frac{D_2}{D_1}=\xi _c^q\frac{F(a,b;c;\xi _c)}{F(1+ac,1+bc;2c;\xi _c)}.$$
(28)
This boundary condition is sometimes called the ‘total reflection’ condition because of the observation that it implies (for real $`\omega `$) that the reflection coefficient at $`z=0`$ equals unity: $`R=|D_2/D_1|^2=1`$. To see why this is so notice that real $`\omega `$ implies $`q`$ is pure real or pure imaginary. Within the domain of the approximation $`k_{}^2c_s^2/\omega ^21`$ it is pure imaginary, so we write $`qi\beta `$ for some real $`\beta `$. With this choice we have $`\sigma (\gamma 1)/(2\gamma )+i\beta /2`$, and so

$$a\sigma +K_{},b\sigma K_{},\text{and}c1+i\beta ,$$

which leads to the relations:

$`1+ac`$ $``$ $`{\displaystyle \frac{\gamma 1}{2\gamma }}+K_{}{\displaystyle \frac{i\beta }{2}}=a^{},`$
$`1+bc`$ $``$ $`{\displaystyle \frac{\gamma 1}{2\gamma }}K_{}{\displaystyle \frac{i\beta }{2}}=b^{},`$
$`2c`$ $``$ $`1i\beta =c^{}.`$ (29)
From these conditions follows the conclusion $`R=1`$.

We next turn to the boundary condition at the top of the radiative zone, $`z_{}0.7R_{}`$. Since our interest is in resonant waves which lie deep in the radiative zone, we will assume $`\xi _{}=\xi (z_{})1`$ when applying this boundary condition.

The behaviour of $`b_z`$ for $`\xi >1`$ is found using Eq. (62) from Appendix A together with Eq. (24):

$`{\displaystyle \frac{b_z}{D_1}}`$ $`=`$ $`[A_1()^cA_3{\displaystyle \frac{D_2}{D_1}}]()^a\xi ^K_{}\times `$ (30)
$`\times F(a,1c+a;1+ab;{\displaystyle \frac{1}{\xi }})+`$
$`+[A_2()^cA_4{\displaystyle \frac{D_2}{D_1}}]()^b\xi ^{+K_{}}\times `$
$`\times F(1c+b,b;1a+b;{\displaystyle \frac{1}{\xi }}),`$
where

$`A_1={\displaystyle \frac{\mathrm{\Gamma }(c)\mathrm{\Gamma }(ba)}{\mathrm{\Gamma }(b)\mathrm{\Gamma }(ca)}},A_2={\displaystyle \frac{\mathrm{\Gamma }(c)\mathrm{\Gamma }(ab)}{\mathrm{\Gamma }(a)\mathrm{\Gamma }(cb)}},`$
$`A_3={\displaystyle \frac{\mathrm{\Gamma }(2c)\mathrm{\Gamma }(ba)}{\mathrm{\Gamma }(bc+1)\mathrm{\Gamma }(1a)}},`$
$`A_4={\displaystyle \frac{\mathrm{\Gamma }(2c)\mathrm{\Gamma }(ab)}{\mathrm{\Gamma }(ac+1)\mathrm{\Gamma }(1b)}}.`$ (31)
For $`\xi 1`$, the two terms in Eq. (30) behave as $`\xi ^{\pm K_{}}e^{\pm K_{}z/H}`$, and so either fall or grow exponentially as functions of $`z`$. We take the absence of the growing behaviour to be our boundary condition for $`\xi =\xi _{}`$. This requires vanishing of the whole factor within brackets in the second line of Eq. (30), thereby providing a second condition which must be satisfied by the ratio $`D_2/D_1`$:

$$\frac{D_2}{D_1}=()^c\frac{\mathrm{\Gamma }(c)\mathrm{\Gamma }(1b)\mathrm{\Gamma }(1+ac)}{\mathrm{\Gamma }(2c)\mathrm{\Gamma }(a)\mathrm{\Gamma }(cb)}.$$
(32)
Given this boundary condition, $`b_z`$ falls exponentially for $`z>z_r`$, $`b_z\xi ^K_{}=\xi _0^K_{}e^{K_{}z/H}`$, as given by the first line in Eq.(30).

### 3.4 Eigenfrequencies

The consistency of the two conditions, eqs. (32) and (28), for $`D_2/D_1`$ gives the eigenvalue condition which fixes the dispersion relation for the resonant waves:

$`\xi _c^q{\displaystyle \frac{F(a,b;c;\xi _c)}{F(1+ac,1+bc;2c;\xi _c)}}=`$
$`=()^c{\displaystyle \frac{\mathrm{\Gamma }(c)\mathrm{\Gamma }(1b)\mathrm{\Gamma }(1+ac)}{\mathrm{\Gamma }(2c)\mathrm{\Gamma }(a)\mathrm{\Gamma }(cb)}},`$ (33)
from which we read off the MHD frequency, $`\omega (k_x,k_y,B_0,n)`$, where the integer $`n`$ is the mode number for the eigenfunction of interest. The solutions which are obtained in this way have the important property that they are complex, $`\omega =\omega _1(1+id)`$, a feature which is also present in the absence of gravitational loading, and which has been exploited in other contexts in attempts to explain the nature of coronal heating .

We note in passing that, in contrast to the master equation in (16), Eq. (3.4) does not reduce to the standard $`g`$-mode spectrum in the limit $`B_00`$. It does not do so because of our assumption that $`N`$ be approximately constant. Although this approximation suffices for the study of MG waves, it is too crude to capture the resonant frequencies of $`g`$-modes. It cannot do so because $`g`$-modes are trapped inside the radiative zone precisely because $`N`$ is not constant, but instead drops to zero at the top of the radiative zone. For MG waves this dropping of $`N`$ is not crucial because its role in trapping the modes is instead played by the resonant layer at $`z=z_r`$.

The eigenvalue condition can be written in a more transparent way by taking advantage of some properties of hypergeometric functions. Recall that our interest is in the regime $`K_{}N/\omega _11`$, and so the parameter $`q`$, as defined by Eq. (22), is mostly imaginary, $`q=i\beta `$, since $`|d|1`$, with $`\beta =2K_{}N/\omega 1`$. This implies that $`\sigma =(\gamma 1)/(2\gamma )+i\beta /2`$ is also large inasmuch as $`|\sigma |1`$, and so must also be $`a`$ and $`b`$. This observation allows us to use the Watson asymptotic form for hypergeometric functions (Appendix A. Some Hypergeometric Facts) – see Appendix A – which applies for large values of $`a`$ and $`b`$.

Using this, together with liberal use of the convolution identity $`\mathrm{\Gamma }(z)\mathrm{\Gamma }(1z)=\pi /\mathrm{sin}(\pi z)`$, allows us to rewrite the eigenvalue condition, Eq. (3.4), as a simple transcendental equation:

$$2i\beta \mathrm{ln}\eta =2i\pi n+\mathrm{ln}\left|\mathrm{tan}\frac{\pi }{\gamma }\right|\pm i\pi \left(\frac{1\gamma }{\gamma }\right),$$
(34)
where the upper (or lower) sign on the right-hand-side is chosen if the imaginary part of the frequency is positive: $`d>0`$ (or negative, $`d<0`$). The variable $`\eta `$ is defined by $`\eta =\xi _c^{1/2}\sqrt{1/\xi _c1}`$, and depends on the magnetic field $`𝐁_0`$ through the Alfvén velocity $`v_{Ac}`$, since $`\xi _c=(k_xv_{Ac}/\omega _1)^2`$. The integer $`n=0,\pm 1,\pm 2,\mathrm{}`$ is the mode number.

The eigenvalue equation is easier to work with if re-expressed slightly. Solving for $`\xi _c(\eta )`$ we find $`\xi _c^{1/2}=2\eta /(1+\eta ^2)=k_xv_{Ac}/\omega `$, and so

$$\beta =\left(\frac{4HN}{\alpha v_{Ac}}\right)\frac{\eta }{1+\eta ^2}$$

where $`\alpha =k_x/k_{}1`$. If we now change variables to $`\chi `$ defined by $`\eta =e^\chi `$ we find

$$\chi =\frac{\alpha v_{Ac}}{4NH}\left[2\pi ni\mathrm{ln}\left|\mathrm{tan}\frac{\pi }{\gamma }\right|\pm \pi \left(\frac{1\gamma }{\gamma }\right)\right]\mathrm{cosh}\chi .$$
(35)
The variable $`\chi `$ is convenient because physical quantities have simple expressions in terms of it. In particular, the phase velocity is $`v_{ph}=\omega _1/k_x=v_{Ac}Re\mathrm{cosh}\chi `$; the frequency’s imaginary part is determined by $`d=Im\mathrm{cosh}\chi /Re\mathrm{cosh}\chi `$, and the position of the resonant layer is $`z_r=2H\mathrm{ln}\left[Re\mathrm{cosh}\chi \right]`$.

Notice that Eq. (35) takes the form $`\chi =A\mathrm{cosh}\chi `$. For real $`A`$ this equation has two real roots for $`\chi `$ if $`A<A_{max}`$ and none if $`A>A_{max}`$, where $`A_{max}`$ is defined as the maximum value taken by the function $`\chi /\mathrm{cosh}\chi `$ (which numerically occurs for $`\chi =\chi _{max}=1.1997)`$. A similar result holds for $`A`$ and $`\chi `$ close to real, such as happens when $`\omega =\omega _1(1+id)`$, and so solutions with $`d`$ small only exist up to a maximum mode number, $`|n|n_{max}`$, where

$$n_{max}=\left(\frac{2NH}{\pi \alpha v_{Ac}}\right)\frac{\chi _{max}}{\mathrm{cosh}\chi _{max}}.$$

Solutions with large $`d`$ must be discarded because they lie beyond the scope of the linearized approximations within which we work.

Figures 2 and 3 present our numerical solution of the eigenvalue spectrum obtained from Eq. (35). Figure 2 plots $`v_{ph}(n)/c_s=\omega _1/(c_sk_x)`$ vs mode number $`n`$ for various magnetic fields, $`B_0`$. Figure 3 similarly plots $`d`$ against mode number for the same magnetic fields. In both plots the parameter $`\alpha =K_x/K_{}`$ is taken to be unity, $`\alpha =1`$. The result with $`\alpha 1`$ is trivially found by rescaling $`B_0`$ because Eq. (35) depends only on the product $`\alpha B_0`$. In both figures a dashed line is used if the resonance of interest occurs above the top of the radiative zone (and so outside the domain of many of our approximations).

Both of the figures 2 and 3 show two branches of solutions up to a maximum mode number, as expected. Their dependence on $`n`$ can also be understood analytically if we take $`\chi \xi _c1`$, and so $`\eta 1`$. In this limit we obtain from (34) the following approximate equation for the MHD spectrum:

$`{\displaystyle \frac{\omega _1}{N}}={\displaystyle \frac{2K_{}}{\pi n}}\left(z_r/2H+\mathrm{ln}2\right),`$ (36)
$`d={\displaystyle \frac{\mathrm{ln}|\mathrm{tan}(\pi /\gamma )|}{2\pi n\left(1(z_r/2H+\mathrm{ln}2)^1\right)}},`$
which leads to the asymptotic solution

$`{\displaystyle \frac{\omega _1}{N}}={\displaystyle \frac{2K_{}}{\pi n}}\mathrm{ln}\left({\displaystyle \frac{4NH}{\alpha v_{Ac}}}{\displaystyle \frac{1}{\pi n}}\right),`$
$`d={\displaystyle \frac{\mathrm{ln}\left|\mathrm{tan}(\pi /\gamma )\right|}{2\pi n}}.`$ (37)
These modes correspond to the branches of the figures 2 and 3 for which both $`v_{ph}/c_s`$ and $`d`$ fall with $`n`$. The other branch can be found by using the approximation $`|\chi |1`$ and hence $`\mathrm{cosh}\chi 1+\chi ^2/2`$. The spectrum in this case is

$`\omega _1={\displaystyle \frac{\alpha K_{}v_{Ac}}{H}},`$
$`d=\left({\displaystyle \frac{\alpha v_{Ac}}{4NH}}\right)^2\pi n\mathrm{ln}\left|\mathrm{tan}{\displaystyle \frac{\pi }{\gamma }}\right|.`$ (38)
Note that for this branch $`v_{ph}(n)/c_s`$ is independent of the mode number, $`n`$, and $`d`$ grows with $`n`$, as is also seen in the figure.

In addition to requiring the resonance to occur in the radiative zone (the solid line in the figures), the validity of our approximations also demand the frequency to not be smaller than $`10^5`$ s-1 due to our neglect of the 27-day period solar rotation. Inspection of Fig. 2 shows that these two conditions are consistent with one another for a reasonably large range of modes only for magnetic fields larger than a kG or so.

## 4 The Resonant Properties

Several striking features emerge from the previous section’s analysis of the hydrodynamic eigenvalues and eigenfunctions. The main two of these are:

1. The eigenfrequencies predicted by Eq. (34) are complex, implying both damped and exponentially-growing modes, and
2. The eigenmodes are not smooth as functions of $`z`$ about the singular resonant point, $`z=z_r(n,k_x,k_y)`$, whose position depends on the quantum numbers of the mode in question.

These features lie at the root of the surprisingly large effects which are produced by even moderate magnetic fields. We have argued that both of these features can be traced to a resonance between $`g`$-modes and Alfvén waves, which arises due to the existence of gradients in the background density, $`\rho _0(z)`$.

In this section we expand on these two properties, to more precisely pin down their origin. In particular, we expect the development of complex frequencies to indicate a hydrodynamical instability. This is much like the imaginary part which $`g`$-waves develop as one passes into the convective zone from the radiative zone. For $`g`$-modes this imaginary part arises because buoyancy no longer acts as a restoring force in the convective zone, but instead becomes destabilizing. Indeed the instability in this case drives the convection, which is why the radiative zone ends. We expect a similar instability to arise for the Alfvén/$`g`$-mode resonance.

### 4.1 Instability

The necessity for complex frequencies is the most disturbing feature about the above eigenvalue problem, because it implies mode functions which grow exponentially in amplitude with time. This signals an instability in the physics which pumps energy into these resonances, and this section aims to discuss the nature of this instability, and the physical interpretation to be assigned to the imaginary part, $`d`$.

Exponentially-growing instabilities within an approximate (e.g. linearized) analysis reflect a system’s propensity to leave the small-field regime, on which the validity of the approximate analysis relies. The question becomes: where does the instability lead, and what previously-neglected effects stabilize the runaway behaviour?

In the present instance we shall see that the normal modes become strongly peaked near the resonant radii, and that the energy flow near these radii is directed along the resonance plane, as would be true for a pure Alfvén wave. This leads us to the following picture. Since helioseismic waves are likely generated by the turbulence at the bottom of the convective zone, it is natural to imagine starting the system with a regular helioseismic $`g`$-mode and asking how it evolves. In this case the resonance allows the energy in this mode to be funnelled into the Alfvén wave, and so to be channelled along the magnetic field lines away from the solar equatorial plane. This process drives us beyond the range of our approximations by making the wave amplitude large near the resonant surface, but also by driving energy out of the equatorial plane where the rectangular analysis applies.

We arrive in this way at an interpretation in which the instability signals the excitation of Alfvén waves from $`g`$-modes due to their resonant mixing. In this picture the imaginary part of the frequency, $`\mathrm{\Gamma }=Im\omega =d\omega _1`$, describes the rate at which the Alfvén mode is excited in this process. Since the mixing and excitation are weak-field hydrodynamic phenomena it is reasonable that this rate can be computed using only the linearized magneto-hydrodynamical equations.

Once excited, the mode amplitudes near the resonances grow until the energy in them becomes dissipated by effects which are not captured by the approximate discussion we present here. The rate for this dissipation must grow as the mode amplitude grows, until it equals the production rate, $`\mathrm{\Gamma }`$, at which point a steady state develops and the mode growth is stabilized. Because the resonant mode grows until its damping rate equals its production rate, it suffices to know the mode’s growth rate in order to determine the overall on-resonance amplitude.

How much more can be said about the amplitude of the modes depends on whether or not the mode stabilizes at small enough amplitudes to allow the hydrodynamic approximation to the production rate, $`\mathrm{\Gamma }=d\omega _1`$, to be accurate. For instance, if the mode only stabilizes once it is large enough to require a nonlinear analysis, then the final production and dissipation rates may be very different from the linearized expressions derived above. If, on the other hand, the modes saturate at comparatively small amplitudes by dissipating energy into non-hydrodynamical modes (such as by Landau damping), then it can happen that the stabilized mode amplitude is not large enough to significantly change the linearized prediction for its production rate.

In what follows we use the imaginary part of $`\omega `$ to determine the resonant mode’s amplitude at the position of the resonance. There is no loss in doing so to the extent that we take $`Im\omega `$ as a parameter which is not related to the mode quantum numbers by the hydrodynamical equations of motion. On occasion we shall also use the explicit formulae for $`Im\omega `$ which is predicted by the hydrodynamics, and in this case we are assuming that the mode stabilization occurs at small enough amplitudes to not invalidate the assumption that this is a good approximation to the mode production rate.

### 4.2 Spatial dependence

We next address the singular shape of the wave fronts at the resonant points, $`z=z_r`$. Mathematically, this singularity arises because the interval of physical interest, $`0zz_{}`$, corresponds to the interval $`0<\xi _c\xi \xi _{}`$, and this contains the singular point $`\xi =1`$ provided $`\xi _c<1`$ and $`\xi _{}>1`$. The hypergeometric functions can – and do – develop singularities at $`\xi =1`$, and it is this singular behaviour which we now examine.

Using Eq. (Appendix B. Asymptotic Forms) of the Appendix B one finds as $`\xi 1`$ (from either side):

$$b_z|\xi 1|^{1\nu }+\text{nonsingular}\text{and}\frac{db_z}{dz}|\xi 1|^\nu ,$$
(39)
with $`\nu =(\gamma 1)/\gamma `$ (giving $`\nu =\frac{2}{5}`$ for an ideal gas, for which $`\gamma =\frac{5}{3}`$). For $`0<\nu <1`$ we see that $`b_z`$ remains finite but $`db_z/dz\mathrm{}`$ as $`\xi 1`$. As we shall see, once the eigenfrequency’s imaginary part is included $`db_z/dz`$ does not actually diverge as $`\xi 1`$, but instead adopts a systematically large resonant line-shape.

Using this information in eqs. (14) allows the asymptotic form for all other quantities to be determined at the resonant point. In particular, since $`v_z`$ does not involve $`db_z/dz`$, it also remains finite on resonance, unlike all other components of the fluid velocity. This indicates that it is the transverse velocity perturbations (and so also the transverse kinetic energy of the oscillations) which increase most dramatically at the resonant plane.

Positions of resonances
The equation for the position of resonances can be found by using the resonance condition

$$\text{Re}\xi =\text{Re}\xi _ce^{z/H}=\left(\frac{k_xv_{Ac}}{\omega _1(n)}\right)^2e^{z_r(n)/H}=1$$
(40)
and Eq. (36). It takes the following transcendental form

$$z_s(n)=H\left[\frac{K_x}{K_{}}\left(\frac{\pi |n|v_{Ac}}{c_s}\right)\frac{\gamma }{\sqrt{\gamma 1}}e^{z_s(n)/2H}\mathrm{ln}4\right].$$
(41)
Taking into account the asymptotic form for the real part of the eigenfrequency, $`\omega _1(n)`$, as given by Eq. (3.4), leads to the following asymptotic expression for the dependence of the resonant layer position on the node number, $`n`$:

$$z_r(n)=2H\mathrm{ln}\left[\frac{4NH}{\alpha v_{Ac}}\frac{1}{2\pi n}\mathrm{ln}\left(\frac{4NH}{\alpha v_{Ac}}\frac{1}{\pi n}\right)\right],$$
(42)
here

$$v_{Ac}=\frac{B_0}{\sqrt{4\pi \rho _c}}=\left(\frac{B_0}{43.4\text{G}}\right)\text{cm s}^1$$
(43)
is the Alfvén velocity at the solar center.

Fig. 4 plots the position, $`z_r(n)`$, predicted by Eq. (42) for the case of longitudinal wave propagation ($`k_y=0`$) and for magnetic fields in the range $`B_0=`$10G–1MGauss. (The same result for an obliquely-propagating wave with $`k_y0`$ is produced by a larger value for $`B_0`$.) This plot is the analog of Figs. 5a,b in .

Spacing of resonances
Knowing the position of the resonances also permits us to determine the distance between them. This quantity is relevant to the propagation of particles like neutrinos through the resonant waves. The spacing is:

$$z_r(n+1)z_r(n)=\left(\frac{K_x}{K_{}}\right)\frac{\pi Hv_{Ac}\gamma }{c_s\sqrt{\gamma 1}}e^{z_r(n)/2H}.$$
(44)
This dependence of this quantity on mode number, $`n`$, is shown in Fig. 5 for longitudinal wave propagation, $`k_y=0`$, and for different magnetic field values. From this figure we see, in the region $`z_r0.3R_{}`$, that the distance between layer positions grows with magnetic field and with distance from the solar center. Both of these features were already seen in preliminary WKB calculations of the resonances (see Figs. 7a, 7b of ref. ), and our present numerics also confirm that the spacing is approximately proportional to $`B_0`$ for $`z0.3R_{}`$, that was seen earlier.

Eq. (44) can be simplified using the expression, Eq. (42), for the resonance position:

$$z_r(n+1)z_r(n)\frac{2H}{n},$$
(45)
where the inverse mode number $`n^1`$ is proportional to the background magnetic field $`B_0`$. This is in accordance with the spectrum Eq. (36), for which the phase velocity $`v_{ph}=\omega _1/k_x`$ equals the Alfvén speed at the resonance position, $`\omega _1/k_x=v_A(z_r,B_0)B_0`$.

Numerically, it is noteworthy that the spacing between resonances can be in the ballpark of hundreds of kilometers. This is significant because this is close to the resonant oscillation length, $`l_{osc}^{res}`$, for $`E\mathrm{MeV}`$ neutrinos, if – as now seems quite likely – resonant LMA oscillations are responsible for explaining the solar neutrino problem since

$$l_{osc}^{res}=\frac{250\mathrm{km}(\mathrm{E}/\mathrm{MeV})}{\mathrm{\Delta }m_5^2\mathrm{sin}2\theta },$$
(46)
where $`\mathrm{\Delta }m_5^2=\mathrm{\Delta }m^2/10^5\text{eV}^2`$. Repeatedly perturbing neutrinos over distance scales comparable to $`l_{osc}^{res}`$ has long been known to be a prerequisite for disturbing the standard MSW picture of oscillations in the solar medium. This raises the possibility – recently explored in more detail in ref.  – that $`g`$-mode/Alfvén resonances can alter neutrino propagation. If so, the observation of resonant oscillations of solar neutrinos may provide some direct information about the properties of the MG waves we discuss here.

Widths of the resonance layers
We next examine the energy flow in the vicinity of a resonance. This gives some information about the nature of the instability which is indicated by the complex mode frequencies. We find here that there is a large energy flow along the magnetic field lines, perpendicular to the $`z`$ direction. We use the dimension of the region in which this energy flows as a measure of the width of the resonant layer.

In the linearized regime the energy flux, $`𝐒`$, has both mechanical and electromagnetic (Poynting) components, $`𝐒=𝐒_{\mathrm{mech}}+𝐒_{\mathrm{em}}`$, with

$`𝐒_{\mathrm{mech}}`$ $`=`$ $`p^{}𝐯,`$ (47)
$`𝐒_{\mathrm{em}}`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}\left\{\left[𝐁^{}\times [𝐯\times 𝐁_0]\right]\right\}.`$
In order of magnitude, the relative size of these contributions is $`B_0B^{}/(4\pi p^{})`$ and this turns out to be small. This smallness is easy to see to the extent that $`B^{}/B_0p^{}/p_0`$, since then $`B_0B^{}/(4\pi p^{})B_0^2/(4\pi p_0)`$, which is small by virtue of the background inequality $`B_0^2/(4\pi )=\gamma (v_A^2/c_s^2)p_0p_0`$.

Since the mechanical energy flux dominates, its direction is simply given by the direction of $`𝐯`$ near the resonance. This may be determined using the asymptotic form as $`\xi 1`$ in the expressions for $`p^{}`$ and $`v_z`$ in eqs. (14). Notice that $`v_z`$ approaches a finite value in this limit, while $`v_x,v_y`$ and $`p^{}`$ all diverge as $`\xi 1`$. This shows that the energy flow, like the fluid flow, is mainly parallel to the resonant plane when evaluated near a resonance.

Using the asymptotic expressions, $`db_z(z)/dz|1\xi |^\nu `$, and $`b_z(z)|1\xi |^{1\nu }+𝒞`$, with constants $``$, $`𝒞`$ and $`\nu =(\gamma 1)/\gamma 0.4`$, we obtain $`v_z|b_z||1\xi |^{1\nu }`$, and

$`p^{},v_x,v_y\left|{\displaystyle \frac{db_z}{dz}}\right|{\displaystyle \frac{||}{|1\xi |^\nu }}=`$ (48)
$`={\displaystyle \frac{||}{|1\xi _r+\xi _r2id|^\nu }}={\displaystyle \frac{||}{[(1\xi _r)^2+4\xi _r^2d^2|^{\nu /2}}},`$
where we have used our definitions – $`\omega =\omega _1(1+id)`$, with $`|d|1`$, and so $`\xi =\xi _r(12id)`$, with $`\xi _r=\left(k_xv_{Ac}/\omega _1\right)^2e^{z_r/H}`$ – to write the real and imaginary parts explicitly. The energy flow near the resonant plane is similarly given by the time average $`𝐒=p^{}𝐯=\frac{1}{2}\mathrm{Re}(p^{}𝐯^{})`$, where the symbol $`\mathrm{}`$ means average over time. Combining the results for $`p^{}`$ and $`𝐯`$ we see that near resonance

$`S_z`$ $``$ $`{\displaystyle \frac{1}{[(1\xi _r)^2+4\xi _r^2d^2]^{\nu /2}}},`$
$`S_x,S_y`$ $``$ $`{\displaystyle \frac{1}{[(1\xi _r)^2+4\xi _r^2d^2]^\nu }}`$ (49)
and so, for instance, $`S_z`$ reaches half of its resonant value when $`\xi =\xi _\pm =1\pm \delta \xi `$ where $`[(1\xi _\pm )^2+4\xi _\pm ^2d^2]^{\nu /2}=2(2d)^\nu `$.

On the other hand, for $`z`$ well beyond the resonant zone, where $`\xi >1`$, we know that $`b_z(z)`$ is given by Eq. (30), which states $`b_z(z)\xi ^K_{}F(a,a^{};1+2K_{};\xi ^1)`$. By virtue of the condition $`a^{}=1c+a`$ we see that this is a real function (in the limit $`d0`$). Consequently, the explicit factors of ‘$`i`$’ in eqs. (14) imply that quantities like the pressure perturbation, $`p^{}(z)`$, and the transverse velocities $`v_{x,y}(z)`$ are pure imaginary. It therefore follows that the $`z`$-component of the energy flux vanishes in the region $`z>z_r`$, outside of the resonance, because $`S_z=\frac{1}{2}\mathrm{Re}(p^{}v_z^{})=0`$, while the transverse components are nonzero in this region, but are falling off exponentially.

This peaking, then vanishing, of $`S_z`$ motivates our choice of the width of the $`S_z(z)`$ resonant line-shape as a measure of the resonance layer width. From this condition we can easy find the results for the widths $`z_\pm `$ (or, equivalently, $`\xi _\pm `$):

$$\delta \xi ^2=(2|d|)^2\left(2^{2/\nu }1\right),$$
(50)
giving our final expression for the width of the Alfvén singular layer:

$$\delta z=z_+z_{}=H\mathrm{ln}\left(\frac{\xi _+}{\xi _{}}\right)4H|d|2^{\gamma /(\gamma 1)}.$$
(51)
Substituting the expression for $`d`$ in terms of mode number, using Eq. (3.4), this becomes

$$\mathrm{\Delta }z_r(n)\frac{2H}{\pi n}2^{\gamma /(\gamma 1)}\mathrm{ln}\left|\mathrm{tan}\frac{\pi }{\gamma }\right|.$$
(52)
Together with Eq. (45) we then get the ratio of the distance between resonances to the width of each resonance,

$$\frac{z_r(n+1)z_r(n)}{\mathrm{\Delta }z_r(n)}=\frac{\pi }{2^{\gamma /(\gamma 1)}\mathrm{ln}\mathrm{tan}\frac{\pi }{\gamma }}0.5.$$
(53)
Figure 6 plots the density profiles of two neighboring resonances in the region $`z_r0.3R_{}`$. In this region their width is approximately twice their separation as measured by the distance between peaks. That is not surprising since the shape of the density profile is the same as in Eq. (48) for other MHD perturbations (see below, Eq. (55)).

Density Profiles
We now give the explicit expression for the density profile predicted for our wave solutions. The basic expression is obtained through the direct substitution of our expression for $`b_z(z)`$ – keeping in mind eqs. (24), (28) and (3.3) – into the last of eqs. (14). We find in this way the following expression for $`\rho ^{}/\rho _0`$, divided by the arbitrary wave normalization, $`D_1`$:

$`{\displaystyle \frac{1}{D_1}}\left({\displaystyle \frac{\rho ^{}(z)}{\rho _0(z)}}\right)=`$
$`={\displaystyle \frac{i\xi _c^{\sigma _1}}{k_xH}}\{({\displaystyle \frac{\gamma 1}{\gamma }}+{\displaystyle \frac{\sigma _1\omega ^2}{k_{}^2c_s^2}})\left({\displaystyle \frac{\xi }{\xi _c}}\right)^{\sigma _1}F(a_1,b_1;c_1;\xi )`$
$`\left({\displaystyle \frac{\gamma 1}{\gamma }}+{\displaystyle \frac{\sigma _2\omega ^2}{k_{}^2c_s^2}}\right)\left({\displaystyle \frac{\xi }{\xi _c}}\right)^{\sigma _2}F(a_2,b_2;c_2;\xi ){\displaystyle \frac{F(a_1,b_1;c_1;\xi _c)}{F(a_2,b_2;c_2;\xi _c)}}`$
$`+{\displaystyle \frac{\xi _c\omega ^2}{k_{}^2c_s^2}}[{\displaystyle \frac{(\sigma _1^2K_{}^2)}{1+i\beta }}\left({\displaystyle \frac{\xi }{\xi _c}}\right)^{\sigma _1+1}F(a_1+1,b_1+1;c_1+1;\xi )`$
$`{\displaystyle \frac{(\sigma _2^2K_{}^2)}{1i\beta }}\left({\displaystyle \frac{\xi }{\xi _c}}\right)^{\sigma _2+1}F(a_2+1,b_2+1;c_2+1;\xi )\times `$
$`\times {\displaystyle \frac{F(a_1,b_1;c_1;\xi _c)}{F(a_2,b_2;c_2;\xi _c)}}]\},`$ (54)
where $`\sigma _1`$ is given by Eq. (22) for $`q=i\beta `$, $`\sigma _2=(\gamma 1)/2\gamma i\beta /2`$; $`a_{1,2}=\sigma _{1,2}+K_{}`$, $`b_{1,2}=\sigma _{1,2}K_{}`$; $`c_{1,2}=1\pm i\beta `$ (see Eq. (3.3)).

The behaviour of the density perturbation near the resonant point, $`z_r`$, takes a simpler form for the Lagrange perturbation

$`\delta \rho (\stackrel{}{r})=\rho ^{}(\stackrel{}{r})+\delta \stackrel{}{r}\rho _0=`$
$`=[\rho ^{}(z)(\delta z/H)\rho _0(z)]e^{i(\omega tk_xxk_yy)},`$ (55)
where the Eulerian perturbation $`\rho ^{}(z)`$ is given by Eq. (4.2). In linear MHD the Lagrangian density perturbation is directly constrained by mass conservation, Eq. (1), $`i\omega \rho _0^1\delta \rho (z)=u`$, where $`u=\mathrm{div}𝐯`$ is the compressibility, given by Eq. (15). One obtains in this way the following Lagrange density perturbation near the resonance, $`zz_r(n)`$:

$$\frac{1}{D_1}\frac{\delta \rho (z)}{\rho _0}C\left[e^{(zz_r(n))/H}(1+2id)1\right]^{(\gamma 1)/\gamma },$$
(56)
where

$$C=\frac{\mathrm{\Gamma }(\gamma ^1)}{\mathrm{\Gamma }(\gamma ^1)}\left(\frac{\gamma v_{ph}(n)}{i\sqrt{\gamma 1}c_s}\right)^{2/\gamma }\left(\frac{\gamma 1}{4\gamma }+\frac{\gamma ^2K_{}^2}{\gamma 1}\right).$$

Here the phase velocity, $`v_{ph}(n)=\omega _1(n)/k_x`$, is given by the solution to the eigen-spectrum, $`\omega =\omega _1(1+id)`$, such as obtained in Eq. ((3.4)), and we include the unknown normalization factor $`D_1`$.

This quantity $`\left|\delta \rho /\rho _0\right|`$ is plotted in Fig. 6 for two neighbouring resonances near $`z=0.3R_{}`$ for the two values of the magnetic field, $`B_0=`$10kG and $`B_0=`$100kG. We see that the presence of the resonances introduces a series of density excursions, located at the resonant positions. The disappearance of the resonant density spikes in the zero-field limit is most easily seen by noticing that in the limit $`B_00`$ the position $`z_r`$ becomes negative, $`z_r=H\mathrm{ln}4<0`$ (see Eq. (41)), or occurs outside the physical region $`0zz_{}`$.

## 5 Possibilities for Observations

Given the strength of the resonant effects on the helioseismic wave profiles, we next explore some of their observable properties to determine how they might be detected, and to what extent their effects are consistent with current helioseismic observations.

### 5.1 Measurements of Mode Frequencies?

Many helioseismic waves have been detected, but to date all of these are pressure-driven $`p`$-modes rather than the usually-lower-frequency $`g`$-modes. The latter are more difficult to detect because they are so strongly attenuated within the convective zone.

Because the central Alfvén frequency is so small compared with the other helioseismic frequencies, we only expect resonances between Alfvén waves and those helioseismic waves which penetrate most deeply into the solar interior, and which have the lowest frequencies. This is why the resonance we have found involves mixing only with $`g`$-modes, and not $`p`$-modes. It follows that so long as observations are restricted to $`p`$-modes, the resonances we describe here will have negligible implications for the comparison between measured $`p`$-wave frequencies and solar models.

Things become more interesting should $`g`$-modes be observed, since then a comparison between observations and predicted frequencies could depend more strongly on the solar magnetic field. Figure (2) (or Eqs. (3.4) and (42)) show how the frequencies of the resonant modes depend on magnetic field. Unfortunately, frequencies this low are likely to be a challenge to detect for quite some time.

Of course, a proper comparison requires the inclusion of the full (spherical) geometry of the problem, as well as temperature and gravity gradients, solar rotation, and other factors which we ignore here. Our goal in this paper is merely to point out the potential for larger-than-expected effects, and to obtain a rough first determination of the sensitivity to magnetic fields.

### 5.2 Constraints on Mode Amplitudes

Mode amplitude remains undetermined within the linear approximations assumed here, and is ultimately controlled by the mechanisms which stimulate and damp the waves of interest in the solar environment. Amplitude can nevertheless be characterized in any of three equivalent ways: the size of the mode’s contribution to the mass-density profile, the size of the fluid velocity, $`v_z`$, at the solar surface, and the total mode energy. In principle, our solutions can be used to express any two of these in terms of the third, as functions of the background magnetic field.

A sufficiently detailed comparison of the resulting expressions – should this become possible – would then allow one to obtain the magnetic field strength. We illustrate this by calculating here the energy content of the resonant MG waves.

Mode Energy
The mode energy is a useful way to quantify the mode amplitude. Here we provide formulae which make this connection explicit. The mode energy is most simply computed by using the virial theorem, which relates the total energy to the kinetic energy averaged over a single period of oscillations: $`E=2E_{\mathrm{kin}}`$:

$$E_{\mathrm{kin}}=\frac{1}{2}d^3𝐱\rho _0[\text{Re}𝐯(𝐱,t)]^2.$$
(57)
For modes having a small imaginary frequency, $`|d|1`$, one obtains the following simple approximate formula for the above average

$`[\text{Re}𝐯(𝐱,t)]^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left|𝐯(z)\right|^2`$
$`+{\displaystyle \frac{1}{4}}[(\text{Im}v_i(z))^2(\text{Re}v_i(z))^2]`$
$`\times \mathrm{sin}[2(k_xx+k_yy)]+\text{O}(d),`$
The second term in this expression is typically negligible once it is integrated over $`dx`$ and $`dy`$.

Using our low frequency approximation, $`\omega ^2k_x^2c_s^2`$, we find the following result for the squared velocity, $`|𝐯(z)|^2=|v_z|^2+|v_x|^2+|v_y|^2`$:

$`|𝐯(\xi )|^2=v_{ph}^2\{|b_z(\xi )|^2+{\displaystyle \frac{\alpha ^2}{K_x^2}}({\displaystyle \frac{|b_z(\xi )|^2}{\gamma ^2}}`$ (59)
$`{\displaystyle \frac{1}{\gamma }}[b_z^{}(\xi )\xi {\displaystyle \frac{db_z(\xi )}{d\xi }}+b_z(\xi )\left(\xi {\displaystyle \frac{db_z(\xi )}{d\xi }}\right)^{}]+\left|\xi {\displaystyle \frac{db_z(\xi )}{d\xi }}|^2\right)\}`$
We next integrate this result over the solar volume. In our simplified geometry we do so by bounding the transverse directions with a circle of radius $`R_{}`$, having cross-sectional area $`\pi R_{}^2`$. This gives the basic expression for the energy per unit length as $`dE/dz=(\pi /4)R_{}^2\rho _0(z)|𝐯(z)|^2`$.

A useful expression is obtained by eliminating the unknown normalization, $`D_1`$ between eqs. (59) and (4.2). This gives a direct, but complicated, relation between the mode energy and the size of the peak density excursion on resonance. To this end we denote the function in brackets in Eq. (59) by $`A(\xi )`$. With this choice we have $`|𝐯|^2=v_{ph}^2(n)A(\xi ,n)`$, with $`v_{ph}^2(n)=\omega ^2/k_x^2=v_A^2(z_r)=B_0^2e^{z_r(n)/H}/(4\pi \rho _c)`$. Note that this last expression uses the resonance condition, $`\xi =(k_xv_A(z_r))^2/\omega ^2=1`$, for the mode eigenfrequencies. If we similarly denote the right-hand-side of Eq. (4.2) by $`B(\xi )`$, so $`\rho ^{}/\rho _0=D_1B(\xi )`$, then we find

$`E_{MHD}(n)=E_0\left({\displaystyle \frac{B_0}{1\text{Gauss}}}\right)^2\times {\displaystyle \frac{1}{(\text{Re}B(\xi =1,n))^2}}\times `$
$`\times \left({\displaystyle \frac{\text{Re}\rho ^{}(z_r(n))}{\rho _0(z_r(n))}}\right)^2\times {\displaystyle _{\xi _c}^\xi _{}}{\displaystyle \frac{d\xi }{\xi ^2}}{\displaystyle \frac{A(\xi ,n)}{D_1^2(n)}},`$ (60)
where for convenience we assume $`D_1`$ to be real. In this expression we have used the explicit background density distribution, $`\rho _0=\rho _ce^{z/H}`$, and changed variables from $`z`$ to $`\xi `$. The lower limit of integration is $`\xi _c=e^{z_r(n)/H}`$ which takes values $`1>\xi _c10^3`$ for a resonance between the solar center and the bottom of the convective zone. We take as upper limit $`\xi _{}=e^{(0.7R_{}z_r)/H}1`$ for the same choice of the Alfvén resonance positions.

The energy scale, $`E_0`$, in Eq. (60) denotes the combination $`E_0=(1/16)R_{}^2H\times (1\text{Gauss})^2`$. Using $`(1\text{Gauss})^2=1`$ erg/cm-3 gives $`E_02.1\times 10^{30}`$ erg. From this we can see that, for magnetic fields $`B_0=10^510^4`$ Gauss, density spikes as large as $`\rho ^{}/\rho _0110`$ per cent will involve energies $`E10^{36}`$ erg, provided $`A/D_1^2`$ and $`B`$ are order unity, which is reasonable for the maximum g-mode energy.

Density Perturbations and Neutrino Oscillations
The g-mode/Alfvén resonances could give rise density perturbations with a large amplitude. Are such density excursions constrained by observations? At first sight one might think so, because the success of the comparison of standard solar models with helioseismic observations constrains the solar density profile to be within 1 per cent of solar model predictions. Since these models ignore magnetic fields, one might expect magnetic-field-induced effects like those we consider here to be ruled out at this level.

Unfortunately, deconvoluting solar properties from observed wave profiles is at present only possible subject to smoothness assumptions. As a result current constraints on deviations of $`\rho `$ from solar model predictions only apply to profiles which do not vary on scales shorter than around 1000 km . It is for this reason that present analyses do not yet exclude the existence of such resonances at an interesting level. Observable effects are only possible if the distance between spikes is comparable to the neutrino oscillation length at the position where neutrino conversions occur. Furthermore, their amplitude must be of order a few percent in order to produce observable effects.

If these density spikes are sufficiently large, they could have implications for neutrino oscillations, and we have argued  that this may be the first place where they are potentially detectable. Besides identifying whether the required density profiles can be obtained using reasonable magnetic fields, a proper analysis must also ask whether the resulting modification of the signal can be identified within the considerable uncertainties which are inherent in any neutrino measurement.

The sensitivity of solar neutrino oscillations to noise in the solar interior has been re-examined in ref. , using the best current estimates of neutrino properties. The outcome is that the measurement of neutrino properties at KamLAND now provides potentially important information on fluctuations in the solar environment on scales $`l_{osc}^{res}`$ (see Eq. (46)) to which standard helioseismic constraints are currently insensitive.

Performing a combined fit of KamLAND and solar data it has also been shown how the determination of neutrino oscillation parameters strongly depends on the magnitude of solar density fluctuations.

The Alfvén/$`g`$-mode resonances turn out to have the right spacing and position for magnetic fields in the range of $`B_0=10100`$ kG. Depending on the size of the magnetic fields which are ultimately found in the radiative zone, this perhaps opens a completely new and interesting observational window on the solar interior.

## 6 Discussion

Here we have presented an analytic discussion of the influence of magnetic fields on helioseismic waves in a simplified geometry. We have argued that this geometry captures the physics not too close to the solar center ($`rH0.1R_{}`$) within the radiative zone, provided that the background magnetic field is perpendicular to the background density gradient. In particular it would apply to a dipole magnetic field configuration near the dipole’s equatorial plane.

We find that within our approximations sufficiently large magnetic fields can cause significant changes to the profiles of helioseismic $`g`$-waves, while not appreciably perturbing helioseismic $`p`$-waves. The comparatively large $`g`$-wave effects arise because of a resonance which occurs between the $`g`$-modes and magnetic Alfvén waves in the solar radiative zone. Although the radiative-zone magnetic fields required to produce observable effects are larger than are often considered – more than a few kG – they are not directly ruled out by any observations. Noting and characterizing the basic features of such resonance constitutes our main result.

Although the density profiles at their maxima could have amplitudes as large as a few percent or more on resonance, we do not believe the corresponding radiative-zone magnetic fields can yet be ruled out by comparison with helioseismic data, since the density excursions are sufficiently narrow (hundreds of kilometers) as to evade the assumptions which underlie standard helioseismic analyses.

For magnetic fields in the 10 kG range, the best hopes for detection of the resonant waves may be through their influence on neutrino propagation. This influence essentially arises because the presence of strong density variations affects the solar neutrino survival probability, with a corresponding change in the resulting solar neutrino fluxes. As shown in ref. , the measurement of neutrino properties at KamLAND provides new information about fluctuations in the solar environment on correlation length scales close to 100 km, to which standard helioseismic constraints are largely insensitive. It has also been shown how the determination of neutrino oscillation parameters from a combined fit of KamLAND and solar data depends strongly on the magnitude of solar density fluctuations.

Since the resonances rely on the condition that $`𝐁\rho `$, there are several magnetic-field geometries to which our analysis might apply, of which we consider here two illustrative extremes. Suppose first that, in spherical coordinates $`(r,\theta ,\varphi )`$, we imagine $`B_r0`$ but $`B_\theta 0`$. Then the field is always perpendicular to a radial density gradient and the resonance we find might be expected to arise in all directions as one comes away from the solar center. In this case the solar $`g`$-modes would tend to be trapped behind the resonance, and so are kept away from the solar surface even more strongly than is normally expected. This would make the prospects for their eventual detection very poor.

Alternatively, if the magnetic field has more of a dipole form it might be imagined that the condition $`𝐁\rho `$ only holds near the solar equatorial plane, and not near the solar poles. In this case a more detailed calculation is needed in order to determine the resulting wave form. This kind of geometry could have interesting consequences for the solar neutrino signal, because in this case the deviation from standard MSW analyses only arises for neutrinos which travel sufficiently close to the solar equatorial plane. Given the roughly 7-degree inclination of the Earth’s orbit relative to the plane of the solar equator, there is a possibility of observing a seasonal dependence in the observed solar neutrino flux. Since the presence of the MG resonance tends to decrease the MSW effect, the prediction would be that the observed rate of solar electron-neutrino events is maximized when the Earth is closest to the solar equatorial plane (December and June) and is minimized when furthest from this plane (March and September).

## Acknowledgements

This work was supported by Spanish grant BFM2002-00345, by European RTN network HPRN-CT-2000-00148, by European Science Foundation network grant N. 86, by Iberdrola Foundation (VBS), by INTAS YSF grant 2001/2-148 and MECD grant SB2000-0464 (TIR). C.B.’s research is supported by grants from NSERC (Canada) and FCAR (Quebec). VBS, NSD and TIR were partially supported by the RFBR grant 00-02-16271. C.B. would like to thank hospitality in Valencia during part of this work.

## Appendix A. Some Hypergeometric Facts

We here record some properties of Hypergeometric functions which are useful in obtaining our results. The series solution to the Hypergeometric equation, Eq. (20), which is analytic throughout the unit disk $`|\xi |<1`$ is:

$$F(a,b;c;\xi )=1+\frac{ab}{c}\left(\frac{\xi }{1!}\right)+\frac{a(a+1)b(b+1)}{c(c+1)}\left(\frac{\xi ^2}{2!}\right)+\mathrm{},$$
(61)
for $`c0,1,2,\mathrm{}`$. The function is defined elsewhere in the complex $`\xi `$ plane by analytic continuation from this series expression.

A convenient way to describe the properties of $`F`$ for $`|\xi |>1`$ is to interchange the interior and exterior of the unit circle using the mapping $`\xi 1/\xi `$. Because this is a particular element of the Hypergeometric group, which maps the three singular points of Eq. (20) ($`\xi =0,1,\mathrm{}`$) among themselves, it has a simple representation on the functions $`F(a,b;c;\xi )`$, given explicitly by the following identity :

$`{\displaystyle \frac{F(a,b;c;\xi )}{\mathrm{\Gamma }(c)}}=`$ (62)
$`={\displaystyle \frac{\mathrm{\Gamma }(ba)}{\mathrm{\Gamma }(b)\mathrm{\Gamma }(ca)}}(\xi )^aF(a,1c+a;1b+a;\xi ^1)+`$
$`+{\displaystyle \frac{\mathrm{\Gamma }(ab)}{\mathrm{\Gamma }(a)\mathrm{\Gamma }(cb)}}(\xi )^bF(b,1c+b;1a+b;\xi ^1).`$
This identity often can be usefully combined with the series expression, Eq. (61) to obtain asymptotic forms as $`|\xi |\mathrm{}`$.

The series solutions given by eqs. (61), (62) are calculated numerically very slowly from both sides in the vicinity of the Alfvén resonance, $`\xi 1`$ and $`\xi 1`$ respectively. For these regions we have used other representations of the Hypergeometric functions . For $`\xi 1`$:

$`{\displaystyle \frac{F(a,b;c;\xi )}{\mathrm{\Gamma }(c)}}={\displaystyle \frac{\mathrm{\Gamma }(cab)}{\mathrm{\Gamma }(ca)\mathrm{\Gamma }(cb)}}\times `$
$`F(a,b;a+bc+1;1\xi )`$
$`+{\displaystyle \frac{\mathrm{\Gamma }(a+bc)}{\mathrm{\Gamma }(a)\mathrm{\Gamma }(b)}}(1\xi )^{cab}\times `$
$`\times F(ca,cb;cab+1;1\xi ),`$ (63)
and for $`\xi 1`$:

$`{\displaystyle \frac{F(a,b;c;\xi )}{\mathrm{\Gamma }(c)}}={\displaystyle \frac{\mathrm{\Gamma }(cab)}{\mathrm{\Gamma }(ca)\mathrm{\Gamma }(cb)}}\times `$
$`F(a,ac+1;a+bc+1;1\xi ^1)`$
$`+{\displaystyle \frac{\mathrm{\Gamma }(a+bc)}{\mathrm{\Gamma }(a)\mathrm{\Gamma }(b)}}(1\xi )^{cab}\xi ^{ac}\times `$
$`F(ca,1a;cab+1;1\xi ^1).`$ (64)
Another asymptotic form of the Hypergeometric function which proves useful for our purposes is the second Watson’s form for large values of the parameters $`a,b`$ and finite $`c`$. We use this in its following variant :

$`F(a+\lambda ,b\lambda ;c;{\displaystyle \frac{1}{2}}{\displaystyle \frac{z}{2}})=`$
$`={\displaystyle \frac{\mathrm{\Gamma }(1b+\lambda )\mathrm{\Gamma }(c)}{\mathrm{\Gamma }(1/2)\mathrm{\Gamma }(cb+\lambda )}}2^{a+b1}(1u)^{c+\frac{1}{2}}\times `$
$`\times (1+u)^{cab\frac{1}{2}}\times \lambda ^{1/2}\left[u^{(\lambda b)}+e^{\pm i\pi (c1/2)}u^{(\lambda +a)}\right]`$
$`\times \left(1+O(\lambda ^1)\right).`$ (65)
Here the argument $`\xi `$ is connected with the auxiliary variable $`z=12\xi `$ and corresponds to $`u=[z+\sqrt{z^21}]^1`$ (if $`Imz<0`$) and to $`u=[z+\sqrt{z^21}]`$ (if $`Imz>0`$) with the same rule for the choice of upper (lower) sign in the exponential factor within brackets.

The correction $`O(\lambda ^1)1`$ turns out to be important for the case $`\lambda =K_{}1`$ which can be used for high frequency spectrum $`\omega _1N`$ (although we do not present this result here). We find the following correction in :

$`O={\displaystyle \frac{1}{2\lambda }}\times {\displaystyle \frac{R+Te^{\pm i\pi (c\frac{1}{2})}u^{ab+2\lambda }}{1+e^{\pm i\pi (c\frac{1}{2})}u^{ab+2\lambda }}},`$
$`R={\displaystyle \frac{1}{2}}{\displaystyle \frac{L+Mu^1+Nu^2}{1u^2}},`$
$`T={\displaystyle \frac{1}{2}}{\displaystyle \frac{L+Mu+Nu^2}{1u^2}},`$
$`L=(a+b2c+1)^2a+b{\displaystyle \frac{1}{2}},`$
$`M=2(a+b1)(a+b2c+1),`$
$`N=(a+b1)^2+ab+{\displaystyle \frac{1}{2}}.`$ (66)
We use this expression obtaining the low frequency spectrum (34), $`\omega _1(n)N`$, and elsewhere in the derivation of approximate forms.

## Appendix B. Asymptotic Forms

The identities of the preceding Appendix A can be used to derive the asymptotic forms in the limit of large $`\xi `$ for the solutions $`Y_1`$ and $`Y_2`$ discussed in the main text. For these purposes we use the following more convenient basis of solutions to Eq. (20):

$`Y_3=\xi ^bF(b+1c,b;a+bc+1;1{\displaystyle \frac{1}{\xi }}),`$
$`Y_4=\xi ^b({\displaystyle \frac{1}{\xi }}1)^{cab}\times `$
$`\times F(ca,1a;cab+1;1{\displaystyle \frac{1}{\xi }}).`$ (67)
The expansion of our original solutions, $`Y_1`$ and $`Y_2`$, in terms of these new ones has the form $`Y_1=N_1Y_3+N_2Y_4`$ and $`Y_2=N_1^{}Y_3+N_2^{}Y_4`$, with coefficients $`N_{1,2}`$ and $`N_{1,2}^{}`$ given by

$`N_1={\displaystyle \frac{\mathrm{\Gamma }(c)\mathrm{\Gamma }(cab)}{\mathrm{\Gamma }(ca)\mathrm{\Gamma }(cb)}},`$
$`N_2={\displaystyle \frac{\mathrm{\Gamma }(c)\mathrm{\Gamma }(a+bc)}{\mathrm{\Gamma }(a)\mathrm{\Gamma }(b)}},`$
$`N_1^{}={\displaystyle \frac{\mathrm{\Gamma }(2c)\mathrm{\Gamma }(cab)}{\mathrm{\Gamma }(1a)\mathrm{\Gamma }(1b)}},`$
$`N_2^{}={\displaystyle \frac{\mathrm{\Gamma }(2c)\mathrm{\Gamma }(a+bc)}{\mathrm{\Gamma }(a+1c)\mathrm{\Gamma }(b+1c)}}.`$ (68)
The large-$`\xi `$ limit is now easily obtained using the result

$$F(a,b;c;1)=\frac{\mathrm{\Gamma }(c)\mathrm{\Gamma }(cab)}{\mathrm{\Gamma }(ca)\mathrm{\Gamma }(cb)},$$

which is finite provided $`Re(cab)>0`$ and $`c0,1,2,\mathrm{}`$ . Notice that these conditions are satisfied for the choices of $`a,b`$ and $`c`$ obtained in the main text. Applying these expressions to $`Y_3`$ and $`Y_4`$ then gives, in the limit $`\xi \mathrm{}`$:

$`Y_3`$ $``$ $`\xi ^b{\displaystyle \frac{\mathrm{\Gamma }(1+a+bc)\mathrm{\Gamma }(ab)}{\mathrm{\Gamma }(a)\mathrm{\Gamma }(a+1c)}},`$
$`Y_4`$ $``$ $`\xi ^b(1)^{cab}{\displaystyle \frac{\mathrm{\Gamma }(1ab+c)\mathrm{\Gamma }(ab)}{\mathrm{\Gamma }(1b)\mathrm{\Gamma }(cb)}}.`$ (69)
Using the previously-given expression, Eq. (68), for the coefficients $`a,b`$ and $`c`$, the limiting form for original solutions $`Y_{1,2}(\xi )`$ finally become:

$`Y_1\xi ^b{\displaystyle \frac{\mathrm{\Gamma }(c)\mathrm{\Gamma }(ab)}{\mathrm{\Gamma }(a)\mathrm{\Gamma }(cb)}}l_1,`$
$`Y_2\xi ^b{\displaystyle \frac{\mathrm{\Gamma }(2c)\mathrm{\Gamma }(ab)}{\mathrm{\Gamma }(1b)\mathrm{\Gamma }(1+ac)}}l_2,`$ (70)
with coefficients $`l_{1,2}`$ given by the elementary functions

$`l_1={\displaystyle \frac{\mathrm{sin}\pi (ca)()^{cab}\mathrm{sin}\pi b}{\mathrm{sin}\pi (cab)}},`$
$`l_2={\displaystyle \frac{\mathrm{sin}\pi a()^{cab}\mathrm{sin}\pi (cb)}{\mathrm{sin}\pi (cab)}}.`$ (71)
Notice that these quantities satisfy the relation $`l_2/l_1=e^{i\pi c}`$.