© ESO, 2011

1. Introduction

The magnetic field in the solar atmosphere tends to accumulate in various entities such as sunspots, pores, spiculae, coronal loops, plumes, etc. The idea that separated magnetic structures can act as channels for wave propagation was recognised a long time ago (see, e.g. Wilson 1980; Roberts 1981a,b; Edwin & Roberts 1983; Roberts et al. 1984; Cally 1986). Stresses created in the convective zone could travel in the form of, e.g. waves, some of them being able to reach even the solar corona. Numerous pieces of evidence for waves in various solar structures have been adduced since high resolution space satellite observation became available. Theoretical models on wave propagation in structured plasmas (e.g. Roberts 1981a,b; Edwin & Roberts 1983; Roberts et al. 1984) were developed within the framework of magnetohydrodynamics (MHD) at a time when observations were sporadic and mainly in the radio-frequency domain (see, e.g. Kai & Takayanagi 1973; Aschwanden et al. 1991). Nowadays theoretical predictions are continuously confirmed by high-resolution observations in (extreme) ultraviolet (for excellent reviews on solar atmospheric waves see, e.g. Nakariakov & Verwichte 2005; Banerjee et al. 2007, and references therein).

With the advent of high resolution observations, coronal seismology emerged as one of the newest branches of modern solar physics, where advances are achieved in diagnosing the structure, the magnitude and sub-resolution structure of the magnetic fields, densities, scale heights, heating functions, etc. (Nakariakov et al. 1999; Ballai 2007; Ruderman et al. 2008) with the help of observable properties of waves and oscillations.

In the theory of MHD waves special attention has been paid to guided waves propagating in various magnetic structures. As a result of propagation in geometrically restricted magnetic “channels”, waves become dispersive with their dispersion relation given in terms of hyperbolic/trigonometric functions for magnetic slabs (Roberts 1981b) and Bessel functions in magnetic cylinders (Edwin & Roberts 1983). Mathematically tractable forms for dispersion relations were obtained in the long wavelength (i.e. for wavelengths much longer than the transversal geometrical size of the waveguide) and short wavelength limits. In general most previous studies considered that the magnetic field is along the density interface. In connection with prominence oscillations Joarder & Roberts (1992) and Oliver et al. (1993) determined the propagation characteristics of waves when the magnetic field was perpendicular (or inclined) to the density jump.

Although significant advances in the theory of guided MHD waves were achieved for waves propagating in a magnetic slab (i.e. Cartesian geometry) and magnetic flux tubes (i.e. cylindrical geometry), very little is known on the properties of guided waves in spherical geometry. The very low number of studies on waves in spherical geometry (compared with other geometries) could be attributed to the fact that for a long time no spherical large scale waves were observed, but also to the mathematical description of wave propagation in this geometry being rather cumbersome. Neglecting the magnetic field, the hydrodynamical approach of spherical waves and oscillations (applicable to, e.g. photospheric p-modes) is well developed and the global oscillations of the solar photosphere are studied within the framework of spherical hydrodynamics (see, e.g. Cox & Guili 1968; Unno et al. 1989; Gough & Toomre 1991, etc.). The lack of observation of global waves (in the sense that the extent of the propagating front is so large that the curvature of the Sun must be considered) has been overcome by recent high resolution observation of global waves in the solar atmosphere by the SoHO/EIT, TRACE, STEREO space satellites. In general, large scale global waves in the solar atmosphere are generated by eruptive events such as flares and/or coronal mass ejections (CMEs) (Biesecker et al. 2002). Such global disturbances were observed propagating in almost every layer of the solar atmosphere, e.g. MDI waves in the photosphere (Kosovichev & Zharkova 1998), Moreton waves in the chromosphere (Moreton & Ramsey 1960), EIT and X ray waves in the corona (Thompson et al. 1999; Narukage et al. 2002).

Any discussion of EIT waves requires some degree of simplification. It is inevitable that a large-scale disturbance, such as an EIT wave, is influenced by the plasma environment. Although most of the investigations suppose a “smooth” environment, in reality EIT waves encounter regions of different densities, active regions with enhanced magnetic field where the waves are reflected/refracted (see, e.g. Terradas & Ofman 2004). In addition, in the light of investigations carried out recently, the nature of these disturbances is far from being elucidated. Several models have been proposed to explain the nature, generation and propagation of these coronal disturbances. It has been suggested that EIT waves are not even real waves, but instead that they are rearrangements of the magnetic structures during eruption of a CME and that this produces currents and pressure enhancements observed as brightenings (e.g. Delanée 2000; Chen et al. 2002). Attrill et al. (2007a,b) proposed that the diffuse EIT coronal bright fronts are due to driven magnetic reconnections between the skirt of the expanding CME magnetic field and the favourably orientated quiet Sun magnetic field. According to this latter model the propagation process of the front consists of a sequence of reconnection events. This scenario was recently reconsidered by Delannée (2009) who concluded that the magnetic topology in the quiet Sun is not adequate for such a series of reconnections.

The interpretation of EIT waves as fast magnetoacoustic waves (FMW) propagating in the quiet Sun perpendicular to the background magnetic field has been put forward by e.g. Wang (2000), Warmuth et al. (2004), Ballai et al. (2005), Long et al. (2008), Patsourakos et al. (2009). In reality it is very likely that both models (magnetic and wave) may coexist and they are just manifestations of the same physical processes at different moments in their evolution. Recent studies (e.g. Temmer et al. 2009; Cohen et al. 2009) suggest that EIT waves are actually a coupling between the CME expansion and MHD waves. When the CME expansion stops the MHD wave component continues to propagate freely in the coronal plasma. The present study assumes that the plasma environment is such that the magnetic field is in the radial direction at all times, however a more realistic approach would be the model where the magnetic field is changing in direction and strength, this assumption requiring a more detailed numerical analysis. A further insight into this aspect would be the investigation of the interaction of an EIT wave with 3D null-points with the purpose for analysing the changes in the properties of EIT waves when interacting with reconnecting field lines. These topics will be the subject of a future investigation.

Analytical approaches aimed at describing the propagation of spherical waves have been carried out in the past by, e.g. Lou (1995a,b, 1996a,b) where the signature of radially propagating magnetoacoustic waves was studied for solar and stellar applications. The radial propagation of weakly nonlinear dissipative spherical linearly polarized Alfvén waves in coronal holes was considered by Nakariakov et al. (2000). Solving numerically the spherical analog of the scalar Cohen–Kulsrud–Burgers equation, these authors showed that the dominant behaviour of the wave is changing with distance, re-confirming similar conclusions found by Lou (1996a,b) for compressional waves. According to the results by Nakariakov et al. (2000), an initial geometrical amplification is followed by a wave breaking and finally an enhanced dissipation.

All these previous studies of wave propagation in spherical geometry described waves propagating in the radial direction. Here our goal is to find an approximative description for waves which propagate in the tangential direction, i.e. waves for which the change in the radial direction is smaller than the change in the direction perpendicular to the radius. In the present paper we start with the supposition that global waves are propagating at a spherical interface in the presence of a radially expanding mean field under the effect of gravity. Although the existence of such density jumps is not observationally obvious (as in the case of coronal loops), however a first approximation of such density interfaces could be the transition region where the density has a drop of nearly two orders of magnitude (see, e.g. the VAL IIIC model by Vernazza et al. 1981). In our analysis we will suppose that global waves observed in the low solar corona are surface gravity waves propagating at this density interface.

The paper is structured in the following way. In Sect. 2 we introduce the fundamental equations and mathematical/physical considerations used later in our analysis. Our equations will be simplified in Sect. 3. considering that the propagation occurs at short distances above the solar surface and that these distances are small compared to the radius of the Sun and gravitational scale-height. The dispersion relation of waves will be derived in Sect. 3 and the variation of the frequency of waves is studied numerically for a series of parameters. Theoretical results are applied to study the propagation characteristics of EIT waves, as a possible application to surface gravity modes at the spherical interface in Sect. 4. Finally, our major results and conclusions are summed up in the Conclusions.

2. Mathematical derivations

The description of wave propagation in spherical geometry with realistic solar magnetic field configurations requires rather complicated mathematical steps often obscuring the physical meaning of the results. Therefore in our analysis we will try to avoid, as much as possible, complicated mathematical approaches. Our investigation will be made considerably easier by assuming that the magnetic field has a monopolar radial configuration and we neglect the effect of rotation which means that the background is spherically symmetric. It should be pointed out that a monopolar magnetic field is a rather simplistic approach but can reasonably approximate the distribution of the magnetic field near the solar minimum. In addition we suppose that the atmosphere of the Sun will remain in hydrostatic and isothermal equilibrium for the whole temporal and spatial domain of interest. As a result all equilibrium quantities (magnetic field, density, pressure) depend solely upon the radial coordinate and the background magnetic field is of the form B₀ = A_B / r²ê_r where is the magnetic flux at the solar surface (up to a factor of 4π) and ê_r is the unit vector in the radial direction. The mathematical derivation used in the set-up of the problem follows the pioneering work in this field by Lou (1995a,b; 1996a,b).

In the isothermal atmosphere we will prescribe the background density, ρ₀(r), and pressure, p₀(r), which together with the ratio of specific heats, γ, will determine the acoustic sound speed c₀ = (γp₀ / ρ₀)^1 / 2, a constant quantity here. The perturbed pressure and density will be denoted by p and ρ, respectively, while the magnetic field perturbation and velocity vectors have the components b = (b_r,b_θ,b_φ) and v = (v_r,v_θ,v_φ).

Since we will distinguish between radial and horizontal dynamics, we split the compression term, ∇·v, into the sum of the radial, Δ_r, and horizontal (tangential), Δ_h, parts where

${\mathrm{\Delta }}_{\mathrm{r}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{{\mathit{\partial }}^{}}{\mathit{\partial }{\mathit{r}}^{}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}{\mathit{v}}_{\mathrm{r}}\mathrm{\right)}\mathit{,}$

${\mathrm{\Delta }}_{\mathrm{h}}\mathrm{=}\frac{\mathrm{1}}{\mathit{r}\sin \mathit{\theta }}\frac{{\mathit{\partial }}^{}}{\mathit{\partial }{\mathit{\theta }}^{}}\left(\sin \mathit{\theta }{\mathit{v}}_{\mathit{\theta }}\right)\mathrm{+}\frac{\mathrm{1}}{\mathit{r}\sin \mathit{\theta }}\frac{{\mathit{\partial }}^{}{\mathit{v}}_{\mathit{\phi }}}{\mathit{\partial }{\mathit{\phi }}^{}}\mathrm{·}$In the forthcoming calculations will also have recourse to the horizontal spherical Laplacian operator

${\mathrm{\nabla }}_{\mathrm{h}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}\sin \mathit{\theta }}\frac{{\mathit{\partial }}^{}}{\mathit{\partial }{\mathit{\theta }}^{}}\left(\sin \mathit{\theta }\frac{{\mathit{\partial }}^{}}{\mathit{\partial }{\mathit{\theta }}^{}}\right)\mathrm{+}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\theta }}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{\phi }}^{\mathrm{2}}}\mathrm{·}$The presence of gravity gives rise to the Brunt- (BV) frequency which will modulate the propagation domain of possible modes and it is defined (for equilibrium quantities depending on the radial direction only) as

${\mathit{N}}^{\mathrm{2}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{0}}}{{\mathit{r}}^{\mathrm{2}}}\left(\frac{\mathrm{1}}{\mathit{\gamma }}\frac{\mathrm{d}\ln {\mathit{p}}_{\mathrm{0}}}{\mathrm{d}\mathit{r}}\mathrm{-}\frac{\mathrm{d}\ln {\mathit{\rho }}_{\mathrm{0}}}{\mathrm{d}\mathit{r}}\right)$where M₀ is the mass of the spherical object contained inside the radius r. According to the standard definition the motion will be oscillatory if N² > 0, i.e.

$\frac{\mathrm{d}\ln {\mathit{p}}_{\mathrm{0}}}{\mathrm{d}\mathit{r}}\mathit{/}\frac{\mathrm{d}\ln {\mathit{\rho }}_{\mathrm{0}}}{\mathrm{d}\mathit{r}}\mathit{>}\mathit{\gamma }$which constitutes the Schwarzschild criterion for convective stability. Since the BV frequency is sensitive to the sharp changes in density and temperature in the transition region, the rapid variation of the BV frequency will act as a reflection filter for many modes propagating in the radial direction. For the particular type of the atmosphere we will consider in the present paper (isothermal) the motion of plasma elements is always stable. We should mention that Ershkovich & Israelevich (2000) investigated the effect of the magnetic field on the value of the BV frequency, but -as the authors point out- these changes appear for cases when the structure of the magnetic field is more complex than the one used here, i.e. when the background magnetic field has a dipolar configuration.

A similar set of equations, but with constant gravitational acceleration, was used by Cally (1983) to study the propagation of waves in stratified plasmas in the context of sunspots. Miles & Roberts (1992) and Miles et al. (1992) also considered waves in gravitationally stratified plasmas, but in Cartesian geometry with a magnetic field along the density interface, while Lou (1995a,b, 1996a,b) extended the analysis of waves propagating in the radial direction in stratified plasmas by considering spherical waves in an unbounded medium permeated by a radially expanding magnetic field.

Let us express all perturbations in the system of MHD equations written for a spherical geometry in the form

$\mathit{f}\mathrm{\left(}\mathit{r,\theta ,\phi ,t}\mathrm{\right)}\mathrm{=}\mathit{f̂}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\exp \mathrm{\left(}\mathrm{i}\mathit{\omega t}\mathrm{\right)}\mathit{,}$where is the real part of a spherical harmonic with l being the spherical degree. Spherical harmonics have the property that

${\mathrm{\nabla }}_{\mathit{h}}^{\mathrm{2}}{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{l}\mathrm{(}\mathit{l}\mathrm{+}\mathrm{1}\mathrm{)}}{{\mathit{r}}^{\mathrm{2}}}{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathit{.}$The assumption we made initially referring to the spherically symmetric background means that the spherical harmonics have an m-degeneracy with a given degree l for the small (i.e. linear) perturbations. Denoting k² = l(l + 1), the system of equations describing the dynamics of waves is transformed into $\mathrm{i}\mathit{\omega }\mathit{\rho ̂}\mathrm{+}\frac{{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{r}}^{\mathrm{2}}}\frac{{\mathrm{d}}^{}}{\mathrm{d}{\mathit{r}}^{}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathit{v̂}{\mathit{r}}^{\mathrm{)}}\mathrm{+}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\Delta ̂}\mathrm{h}\mathrm{+}\mathit{v̂}\mathit{r}\frac{{\mathrm{d}}^{}{\mathit{\rho }}_{\mathrm{0}}}{\mathrm{d}{\mathit{r}}^{}}\mathrm{=}\mathrm{0}\mathit{,}$(1)$\mathrm{(}{\mathrm{\nabla }}_{\mathrm{h}}\mathrm{·}{b}\mathrm{h}\mathrm{)}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{{\mathrm{d}}^{}}{\mathrm{d}{\mathit{r}}^{}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathit{b̂}{\mathit{r}}^{\mathrm{)}}\mathit{,}$(2)$\mathrm{i}\mathit{\omega }\mathit{b̂}\mathit{r}\mathrm{=}\mathrm{-}{\mathit{B}}_{\mathrm{0}}\mathrm{\Delta ̂}\mathrm{h}\mathit{,}$(3)$\mathrm{i}\mathit{\omega }{\mathit{\rho }}_{\mathrm{0}}\mathit{v̂}\mathit{r}\mathrm{=}\mathrm{-}\frac{{\mathrm{d}}^{}\mathit{p̂}}{\mathrm{d}{\mathit{r}}^{}}\mathrm{-}\frac{\mathit{\rho ̂}\mathit{GM}}{{\mathit{r}}^{\mathrm{2}}}\mathit{,}$(4)$\mathrm{i}\mathit{\omega }{\mathit{\rho }}_{\mathrm{0}}\mathrm{\Delta ̂}\mathit{h}\mathrm{=}\frac{{\mathit{k}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}\mathit{p̂}\mathrm{+}\frac{{\mathit{B}}_{\mathrm{0}}}{\mathit{\mu }}\frac{{\mathit{k}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}\mathit{b̂}\mathit{r}\mathrm{-}\frac{{\mathit{B}}_{\mathrm{0}}}{\mathit{\mu }{\mathit{r}}^{\mathrm{2}}}\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{\mathit{r}}^{\mathrm{2}}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathit{b̂}{\mathit{r}}^{\mathrm{)}}\mathit{,}$(5)$\mathrm{i}\mathit{\omega }\mathit{p̂}\mathrm{+}\mathit{v̂}\mathit{r}\frac{{\mathrm{d}}^{}{\mathit{p}}_{\mathrm{0}}}{\mathrm{d}{\mathit{r}}^{}}\mathrm{+}\frac{\mathit{\gamma }{\mathit{p}}_{\mathrm{0}}}{{\mathit{r}}^{\mathrm{2}}}\frac{{\mathrm{d}}^{}}{\mathrm{d}{\mathit{r}}^{}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathit{v̂}{\mathit{r}}^{\mathrm{)}}\mathrm{+}\mathit{\gamma }{\mathit{p}}_{\mathrm{0}}\mathrm{\Delta ̂}\mathrm{h}\mathrm{=}\mathrm{0}$(6)where G is the gravitational constant, M the solar mass and the horizontal component of the magnetic field perturbation has been eliminated by means of the divergence-free condition imposed on the magnetic field perturbation. For simplicity here we assumed that the perturbations are adiabatic. Our gravitational acceleration (MG / r²) decreases with radial distance r. For simplicity, the “hat” symbol will be dropped in the forthcoming calculations. Supposing real values for wavenumbers, the frequency ω will be complex. As a result outward propagating waves will have a growing amplitude with time while those waves which propagate towards the centre of the sphere will have a decreasing amplitude with time. For the particular choice of time dependence employed in the present study (Fourier analyzing with respect to time and supposing that all perturbations oscillate at the same frequency ω) we cannot perform an initial value analysis. However, in an earlier study, Douglas & Ballai (2007) investigated the problem for a simplified atmosphere. In their numerical investigation a sudden energy release could generate initially a non-harmonic disturbance which later decayed into a wave-like propagating feature. They also showed that at the location of the generation of global waves a strong pressure gradient is developed which drives a siphon-like flow oriented towards exterior. This flow could account for the dense and cold plasma (coming from the chromosphere) filling up the dimming region formed after the lift-off of a CME.

In general the description of wave propagation in plasmas where the background equilibrium quantities have a spatial dependence is cumbersome and analytical solutions are not easy to find. Remarkable advances have been achieved in numerical simulations of wave propagation but in this case the variation of perturbations over the domain of interest is not clear and many fine aspects are in general omitted. Serious drawbacks of numerical investigations include the nature of the boundaries and the numerical diffusion which can dramatically change the nature of the investigated physics. Here we will attempt to find approximative formulae to describe the propagation of waves and results will be numerically represented in conjunction with solar applications. The key ingredient in our quest will be the proximity of the domain of interest to the solar surface at distances which are much smaller than any spatial scale involved in the problem (radius of the Sun and gravitational scale-height).

In order to tackle the problem from a mathematical point of view we need to specify the radial dependence of the background density and pressure. Given the particular choice of the background magnetic field, it will not play any role in the equilibrium of forces in the radial direction, so the equilibrium is fully described by the relation

$\frac{{\mathrm{d}}^{}{\mathit{p}}_{\mathrm{0}}}{\mathrm{d}{\mathit{r}}^{}}\mathrm{=}\mathrm{-}\frac{{\mathit{\rho }}_{\mathrm{0}}\mathit{GM}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{·}$Supposing ρ_⊙ to be the density at the solar surface, the radial dependence of the equilibrium density in an isothermal atmosphere is ${\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{\odot }}\exp \left[\mathrm{\Lambda }\left(\frac{\mathrm{1}}{\mathit{r}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{R}}_{\mathrm{\odot }}}\right)\right]\mathit{,}$(7)where is the gravitational length scale in an isothermal atmosphere, in which , m_H, k_B and T₀ are the mean molecular weight, proton mass, Boltzmann constant and temperature, respectively. For typical low coronal temperatures (T₀ = 10⁶ K, ), Λ ≈ 11.05   R_⊙. The background pressure p₀(r), with initial condition p₀(R_⊙) = p_⊙, is similarly given by

${\mathit{p}}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{p}}_{\mathrm{\odot }}\exp \left[\mathrm{\Lambda }\left(\frac{\mathrm{1}}{\mathit{r}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{R}}_{\mathrm{\odot }}}\right)\right]\mathrm{·}$Since the height dependence of the equilibrium pressure and density is the same, the sound speed will be independent of height. At the same time, the Alfvén speed will depend on the radial distance and its form is given by

${\mathit{v}}_{\mathrm{A}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{B}}_{\mathrm{\odot }}{\mathit{R}}_{\mathrm{\odot }}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}\mathrm{(}\mathit{\mu }{\mathit{\rho }}_{\mathrm{\odot }}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\exp \left[\frac{\mathrm{\Lambda }}{\mathrm{2}}\left(\frac{\mathrm{1}}{{\mathit{R}}_{\mathrm{\odot }}}\mathrm{-}\frac{\mathrm{1}}{\mathit{r}}\right)\right]\mathrm{·}$With background pressure and density defined it is now easy to show that the combination of Eqs. (1), (4) and (6) leads to $\left(\mathrm{i}\mathit{\omega }{\mathit{\rho }}_{\mathrm{0}}\mathrm{-}\frac{\mathrm{i}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}{\mathrm{\Lambda }}^{\mathrm{2}}{\mathit{p}}_{\mathrm{0}}}{\mathit{\gamma \omega }{\mathit{r}}^{\mathrm{4}}}\right){\mathit{v}}_{\mathit{r}}\mathrm{=}\mathrm{-}\frac{{\mathrm{d}}^{}\mathit{p}}{\mathrm{d}{\mathit{r}}^{}}\mathrm{-}\frac{\mathrm{\Lambda }\mathit{p}}{\mathit{\gamma }{\mathit{r}}^{\mathrm{2}}}\mathrm{·}$(8)The above equation will allow us to determine the pressure perturbation in terms of the radial component of velocity using the integrating factor technique (see Appendix). Once the pressure perturbation is found, we can also determine the density perturbation in terms of the radial velocity perturbation, by considering Eqs. (1) and (6) to arrive at $\mathit{\rho }\mathrm{=}\frac{\mathit{p}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{i}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{\Lambda }{\mathit{p}}_{\mathrm{0}}}{\mathit{\omega }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}}{\mathit{v}}_{\mathit{r}}\mathit{.}$(9)Another relation which will be used later when deriving the dispersion relation is the radial component of the magnetic field perturbation. Combining Eqs. (3) and (5) we arrive at $\frac{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{\mu }{\mathit{r}}^{\mathrm{2}}}\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{\mathit{r}}^{\mathrm{2}}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}{\mathit{b}}_{\mathit{r}}\mathrm{\right)}\mathrm{-}\frac{{\mathit{k}}^{\mathrm{2}}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{\mu }{\mathit{r}}^{\mathrm{2}}}{\mathit{b}}_{\mathit{r}}\mathrm{+}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\omega }}^{\mathrm{2}}{\mathit{b}}_{\mathit{r}}\mathrm{=}\frac{{\mathit{k}}^{\mathrm{2}}{\mathit{B}}_{\mathrm{0}}}{{\mathit{r}}^{\mathrm{2}}}\mathit{p}\mathit{.}$(10)Similarly to Lou (1995a,b, 1996a,b) we introduce a new dimensionless variable x = 2r / Λ, so that

$\frac{{\mathrm{d}}^{}}{\mathrm{d}{\mathit{r}}^{}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{\Lambda }}\frac{{\mathrm{d}}^{}}{\mathrm{d}{\mathit{x}}^{}}\mathrm{·}$Since we are interested in dynamics near the solar surface, it is clear that the new variable, x, is small. In terms of this dimensionless spatial parameter Eqs. (8) and (9) become $\frac{\mathrm{d}\mathit{p}}{\mathrm{d}\mathit{x}}\mathrm{+}\frac{\mathrm{2}}{\mathit{\gamma }{\mathit{x}}^{\mathrm{2}}}\mathit{p}\mathrm{=}\mathrm{-}\frac{\mathrm{i\Lambda }\mathit{\omega }{\mathit{\rho }}_{\mathrm{0}}}{\mathrm{2}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\delta }}^{\mathrm{2}}}{{\mathit{x}}^{\mathrm{4}}}\right){\mathit{v}}_{\mathit{r}}\mathit{,}$(11)and $\mathit{\rho }\mathrm{=}\frac{\mathit{p}}{{\mathit{c}}_{\mathrm{S}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{4}\mathrm{i}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{p}}_{\mathrm{0}}}{{\mathit{c}}_{\mathrm{S}}^{\mathrm{2}}\mathit{\omega }{\mathit{x}}^{\mathrm{2}}\mathrm{\Lambda }}{\mathit{v}}_{\mathit{r}}\mathit{,}$(12)where we have introduced the notation

$\mathit{\delta }\mathrm{=}\frac{{\mathit{c}}_{\mathrm{S}}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathit{\gamma \omega }}\left(\frac{\mathrm{4}}{\mathrm{\Lambda }}\right)\mathrm{·}$In addition, for small values of x, Eq. (10) reduces to ${\mathit{\rho }}_{\mathrm{0}}{\mathit{\omega }}^{\mathrm{2}}{\mathit{b}}_{\mathit{r}}\mathrm{=}\frac{\mathrm{4}{\mathit{k}}^{\mathrm{2}}{\mathit{B}}_{\mathrm{0}}}{{\mathit{x}}^{\mathrm{2}}{\mathrm{\Lambda }}^{\mathrm{2}}}\mathit{p}\mathit{.}$(13)For arbitrary values of x, Lou (1996a) reduced the MHD equations to two coupled ordinary differential equations for p and b_r which were solved in the limits x ≪ 1 and x ≫ 1. In an interface-free atmosphere, near the solar surface (small values of x), and in the low frequency limit (i.e. frequencies smaller than the BV frequency) v_r is already known to the leading order (see, e.g. Lou 1996a) as ${\mathit{v}}_{\mathit{r}}\mathrm{=}\mathrm{\pm }\mathrm{\left(}\mathrm{𝒜}\mathit{,}\mathrm{ℬ}\mathrm{\right)}{\mathit{c}}_{\mathrm{s}}\mathit{k}\frac{\mathrm{1}}{\sqrt{\mathit{x}}}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathit{x}}\mathrm{\mp }\frac{\mathrm{i}\mathit{k\delta }}{\mathrm{2}{\mathit{x}}^{\mathrm{2}}}\right)\mathit{,}$(14)where and ℬ are constants of integration used later in our derivation. These solutions are obtained using a combination of series techniques and a limiting argument for small x. The upper signs denote a wave travelling inward; conversely the lower signs signify a wave propagating outward.

We should note that in the expression of v_r imported from previous studies there is no information on the magnetic field, so in this approximation (i.e. close to the solar surface) waves are behaving like acoustic- gravity waves.

3. Dispersion relation of spherical surface waves

In order to derive the dispersion relation for spherical surface waves, let us suppose that there is a density jump situated at the dimensionless distance x₀ in the radial direction above the solar surface. A schematic picture of the equilibrium configuration is shown in Fig. 1.

The radial variation of density in terms of the dimensionless quantity x is given as

${\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ {\mathit{\rho }}_{\mathrm{\odot }}\exp \left(\frac{\mathrm{2}}{\mathit{x}}\mathrm{-}\frac{\mathrm{2}}{{\mathit{x}}_{\mathrm{\odot }}}\right)\mathit{,}& \mathit{x}\mathit{<}{\mathit{x}}_{\mathrm{0}}\mathit{,}\\ \mathit{\alpha }{\mathit{\rho }}_{\mathit{i}}\exp \left(\frac{\mathrm{2}}{\mathit{x}}\mathrm{-}\frac{\mathrm{2}}{{\mathit{x}}_{\mathrm{0}}}\right)\mathit{,}& \mathit{x}\mathit{>}{\mathit{x}}_{\mathrm{0}}\mathit{,}\end{array}$

where α ∈ (0,1) represents the “strength” of the density jump with α = 1 corresponding to the situation where the discontinuity is absent while α = 0 represents a case in which the domain of interest is limited to x₀. In the above relation ρ_i is the density at the interface defined as

${\mathit{\rho }}_{\mathrm{i}}\mathrm{=}{\mathit{\rho }}_{\mathrm{\odot }}\exp \left({\frac{\mathrm{2}}{\mathit{x}}}_{\mathrm{0}}\mathrm{-}\frac{\mathrm{2}}{{\mathit{x}}_{\mathrm{\odot }}}\right)\mathrm{·}$Figure 2 shows the variation of the background density around the region of the surface of discontinuity for a range of α between 0.1 and 1.

The introduction of a density interface mandates the introduction of surface boundary conditions, specifically that the continuities of both momentum and energy are satisfied across the surface. Based on the derivation of Lou (1996a) the radial components of the velocity in the two domains now take the form

${\mathit{v}}_{\mathit{r}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \mathrm{𝒜}{\mathit{c}}_{\mathrm{S}\mathrm{1}}\mathit{k}\frac{\mathrm{1}}{\sqrt{\mathit{x}}}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathit{x}}\mathrm{-}\frac{\mathrm{i}\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}{\mathit{x}}^{\mathrm{2}}}\right)\mathit{,}& \mathit{x}\mathit{<}{\mathit{x}}_{\mathrm{0}}\mathit{,}\\ \mathrm{-}\mathrm{ℬ}{\mathit{c}}_{\mathrm{S}\mathrm{2}}\mathit{k}\frac{\mathrm{1}}{\sqrt{\mathit{x}}}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathit{x}}\mathrm{+}\frac{\mathrm{i}\mathit{k}{\mathit{\delta }}_{\mathrm{2}}}{\mathrm{2}{\mathit{x}}^{\mathrm{2}}}\right)\mathit{,}& \mathit{x}\mathit{>}{\mathit{x}}_{\mathrm{0}}\mathit{,}\end{array}$where the subscripts “1” and “2” denote values beneath and above the interface, respectively, and and ℬ are the constants of integration mentioned earlier. The relations between c_S1 and c_S2, and δ₁ and δ₂, are

$\mathit{\alpha }{\mathit{c}}_{\mathrm{S}\mathrm{2}}^{\mathrm{2}}\mathrm{=}{\mathit{c}}_{\mathrm{S}\mathrm{1}}^{\mathrm{2}} \mathrm{and} \mathit{\alpha }{\mathit{\delta }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{=}{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}\mathit{.}$

Continuity of the radial momentum across the surface allows us to express in terms of ℬ, giving $\mathrm{𝒜}\mathrm{=}\mathrm{-}\mathrm{ℬ}\frac{{\mathit{c}}_{\mathrm{S}\mathrm{2}}}{{\mathit{c}}_{\mathrm{S}\mathrm{1}}}\exp \left(\frac{\mathrm{i}\mathit{k}}{\mathrm{2}{\mathit{x}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{(}{\mathit{\delta }}_{\mathrm{1}}\mathrm{+}{\mathit{\delta }}_{\mathrm{2}}\mathrm{)}\right)\mathit{.}$(15)Due to the physical setup of the problem (contact discontinuity), apart from the radial component of the velocity which is continuous at the interface, all other quantities will be also continuous except the density. It can be easily shown that a continuity imposed on, e.g. the radial component of the magnetic field perturbation would automatically imply the continuity of the pressure perturbation, etc.

Equation (13) can be used to derive the expression of the radial component of the magnetic field perturbation in the form ${\mathit{b}}_{\mathit{r}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\mathrm{16}{\mathit{k}}^{\mathrm{2}}{\mathit{A}}_{\mathit{B}}}{{\mathit{\rho }}_{\mathrm{0}}{\mathit{x}}^{\mathrm{4}}{\mathrm{\Lambda }}^{\mathrm{4}}{\mathit{\omega }}^{\mathrm{2}}}\mathit{p}\mathit{.}$(16)Obviously, in the two domains all quantities must be evaluated according to the profile of the background density and radial component of the velocity.

The continuity of b_r(x) at x = x₀ requires that the condition ${\mathit{p}}_{\mathrm{1}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{{\mathit{p}}_{\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\mathit{\alpha }}$(17)is satisfied, where the expressions for the kinetic pressure are given in the appendix. Combining the continuity of the radial displacement and the radial component of the magnetic field perturbation we obtain that the dispersion relation of waves propagating along the density interface at small radial distances is given by $\mathit{\omega }\mathrm{=}\mathrm{\pm }\frac{\mathrm{1}}{{\mathit{x}}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathit{\gamma }}\left(\frac{\mathrm{4}}{\mathrm{\Lambda }}\right)\sqrt{{\mathit{c}}_{\mathrm{S}\mathrm{1}}{\mathit{c}}_{\mathrm{S}\mathrm{2}}}\mathit{.}$(18)It is obvious that the dispersion relation described by Eq. (18) is acoustic in nature despite the presence of the radially expanding background magnetic field. This sort of behaviour was expected given the approximative form of the radial velocity given by Eq. (14). The dispersion relation derived above is also independent of the spherical degree l.

Equation (18) will be used in the next section to draw conclusions on the propagation of EIT waves modelled as fast magnetoacoustic-gravity waves propagating in the low corona. The waves we plan to study here are magnetoacoustic-gravity waves propagating perpendicularly to the radial magnetic field near the solar surface. Before applying the dispersion relation to study the propagation characteristics of global waves in the solar corona, let us make one more remark. Stable oscillations will occur if the frequency of oscillations is smaller than the BV frequency. In terms of the dimensionless spatial length, the BV frequency will be given by

${\mathit{N}}_{\mathrm{BV}}\mathrm{=}\frac{\mathrm{4}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{c}}_{\mathrm{S}}}{\mathit{\gamma }\mathrm{\Lambda }{\mathit{x}}^{\mathrm{2}}}\mathit{,}$which will be estimated at the interface. The results of numerical investigations of the following section (frequency with respect to other parameters) will always be checked to be less than the corresponding BV frequency.

4. Application: EIT waves as surface gravity waves propagating at the density interface

Global waves have been known since the early 1960s. Although it is still unknown how the release of energy and energized particles will transform into waves, today it is widely accepted that these disturbances are similar to the circularly expanding bubble-like shocks after atomic bomb explosions or shock waves which follow the explosion of a supernova. Thanks to the available observational facilities, global waves have been observed in a range of wavelengths in different layers of the solar atmosphere. In the corona, a flare or CME can generate an EIT wave (Thompson et al. 1999), first seen by the SOHO/EIT instrument, or an X-ray wave seen in SXT (Narukage et al. 2002). There is still a vigorous debate about how these fifferent global waves are connected (if they are, at all). Co-spatial and co-temporal investigations of various global waves have been carried out but without a final widely accepted result being reached.

EIT waves are observed to propagate in the quiet Sun with speeds of 250–400 km s^-1 at an almost constant altitude (Thompson et al. 1999; Wills-Davey & Thompson 1999; Klassen et al. 2000). At a later stage in their propagation EIT waves can be considered as freely propagating wavefronts which are observed to interact with coronal loops (see, e.g. Ballai et al. 2008). Using TRACE/EUV 195 Å observations, Ballai et al. (2005) have shown that EIT waves (seen in this wavelength) are waves with average periods of the order of 400 s. Since, at this height, the magnetic field can be considered vertical, EIT waves were interpreted as fast MHD waves (also supported by other observations, e.g. Long et al. 2008; Patsourakos et al. 2009). A great deal of understanding of the properties of EIT waves and their role in generating oscillations of coronal loops came after extensive 3D numerical simulations of MHD wave propagation in the solar corona (see, e.g. Ofman & Thompson 2002; Cohen et al. 2009; Schmidt & Ofman 2010).

There is little known about the propagation characteristics of EIT waves, mainly because the instruments used for their observation (the EIT instruments onboard SOHO and STEREO) have poor temporal resolution (although full disk imagers) while the EUV imager onboard TRACE has a very limited field of view, despite its excellent temporal resolution. Once EIT waves were observed to propagate in a rather localised way in the low corona, the theoretical description of global spherical waves became a necessity.

EIT waves have been observed to propagate at relatively constant altitude corresponding to 0.08   R_⊙ = 56000 km above the solar surface (Mann et al. 1999), which means that in our analysis x₀ = 0.195, i.e. our small x approximation is fully justified.

Although we started with a model which included terms describing magnetic field and gravity, the dispersion relation given by Eq. (18) contains no information about the magnetic field, i.e. near the solar surface these waves behave as acousto-gravity waves.

The dispersion relation is studied numerically for a range of plasma configurations. The density contrast at the interface, α, the location of the density interface, x₀, and the sound speed in the region “2” enter the dispersion relation as free parameters, so our dispersion curves will be investigated in terms of these quantities.

Our first aim is to investigate the effect of varying the intensity of the density jump, from a very slight jump (values of α close to 1) to an almost total jump (small values of α), and the sound speed in the second region (varied here from 130 to 430 km s^-1) on the possible frequencies of global waves propagating along the density interface with the interface being situated at the dimensionless distance x₀ = 0.195 above the solar surface (see Fig. 3).

The solutions were obtained with the condition that the allowed solution is below the corresponding BV frequency and above the acoustic cut-off. The first feature to notice is the absence of solutions for relatively large density contrast and relatively small sound speed. For a smaller value of the density jump, the frequency of the waves is observed to increase with increasing sound speed c_S2. The largest value of the frequency (2.86 mHz) is obtained for small values of density jump and and large values of c_S2. The obtained maximum frequency corresponds to a period of about 36.5 min.

The second case under investigation is the study of the solution of the dispersion relation keeping the sound speed in the region above the interface constant at 300 km s^-1, but now varying the strength of the density contrast and the location of the density interface which is now in the range (in dimensionless units) 0.18 to 0.27 (see Fig. 4).

The largest value of the frequency is attained for relatively small value of the density jump and for the density interface closer to the solar surface.

Before we conclude, let us make a comment. In all investigated cases we obtained frequencies less than the value observed by Ballai et al. (2005) using results from the TRACE/EUV telescope. This discrepancy could arise from the rather simplistic model we used for our investigation. Describing the magnetic field by a radially expanding r-dependent quantity imposes serious limitations; the realistic magnetic configuration in the quiet Sun is much more complicated. A more realistic height dependence for density and pressure would also change considerably the obtained frequency domain.

5. Conclusions

The theory of guided waves along magnetic and density structures in the solar atmosphere developed earlier for Cartesian and cylindrical geometries was extended to include the study of wave propagation along a spherical density interface in the solar corona. The plasma was supposed to be a gravitationally stratified isothermal fluid permeated by a radially expanding magnetic field. Due to the complexity of the mathematical description, our analysis was restricted to a domain close to the solar surface and the extent of the domain of interest was considered short compared with the radius of the Sun and the gravitational scale-height.

The theory developed in the present study was intended to mathematically and numerically describe the propagation of recently observed global coronal waves (also known as EIT waves) in the low solar corona treating them as fast magnetoacoustic waves propagating along a spherical density interface. The knowledge of propagation characteristics of EIT waves is essential as global waves can be used for high accuracy global coronal seismology of the quiet Sun.

The dispersion relation was derived using an approximative method, and the solutions were studied in terms of three parameters (strength of the density jump, sound speed in the region above the interface, and the location of the interface). All investigated cases show that solutions of the dispersion relation become acceptable once the density contrast between the two regions is tending towards values close to 1, representing the case of very weak density jump. The solutions found numerically were bounded from above by the Brunt- (BV) frequency. Due to the approximations used in our analysis, the waves described by the dispersion relation (18) are acousto-gravity in nature.

Given the discrepancy between the measured and derived values of frequencies, the results obtained here also mean that there are still some fundamental effects which need to be added to our working model.

In terms of magnetoseismological purposes, the dispersion relation we have obtained allows the inverse problem to be solved, i.e. once propagation characteristic of waves are observed, a straightforward reconstruction of the plasma parameters can be found.

The results of the present study will be expanded in the near future to the investigation of wave propagation in more complicated structures such as the plasma layer (or spherical shells), to include more realistic physics such as non-ideal effects, and make use of ever more accurate approximations. We will also try to study the conditions which need to be satisfied in order to trap EIT waves in the magnetic spherical shell in the low solar corona.

Acknowledgments

M.D. acknowledges the support of STFC. E. F.-D. acknowledges the support of the EC Solaire Network (MRTN-CT-2006-035484). I.B. was supported by The National University Research Council Romania (CNCSIS-PN-II/531/2007). E.F.-D. and I.B. acknowledge the support by the Hungarian Science Research Fund (OTKA) under grant No. K67746 and  K81421. Part of the present study was carried out while I.B. was a visitor of Eötvös University. I.B. gratefully acknowledges the kind hospitality of his Host.

References

Appendix A: Obtaining solutions for the pressure perturbation in the x  ≪  1 limit

The Appendix is devoted to the mathematical background calculations employed when calculating the dependence of gas pressure on the dimensionless quantity x in the limit of propagation near the solar surface, i.e. x ≪ 1. The starting point of our analysis is Eq. (11), which is a first order linear differential equation to be solved for pressure, p(x). It is easy to show that the pressure p(x) is given by $\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{i}\mathit{\omega }{\mathit{\rho }}_{\mathrm{\odot }}\mathrm{\Lambda }{\mathrm{e}}^{\mathrm{2}\mathit{/}\mathit{\gamma x}}{\mathit{c}}_{\mathrm{S}\mathrm{1}}\mathit{k}}{\mathrm{2}{\mathrm{e}}^{\mathrm{2}\mathit{/}{\mathit{x}}_{\mathrm{\odot }}}}{\mathit{I}}_{\mathrm{int}}$(A.1)for the internal region and $\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\mathrm{i}\mathit{\omega \alpha }{\mathit{\rho }}_{\mathit{i}}\mathrm{\Lambda }{\mathrm{e}}^{\mathrm{2}\mathit{/}\mathit{\gamma x}}{\mathit{c}}_{\mathrm{S}\mathrm{2}}\mathit{k}}{\mathrm{2}{\mathrm{e}}^{\mathrm{2}\mathit{/}{\mathit{x}}_{\mathrm{0}}}}{\mathit{I}}_{\mathrm{ext}}$(A.2)in the region x > x₀. In the last two relations I_int and I_ext are given by ${\mathit{I}}_{\mathrm{int}}\mathrm{=}\mathrm{\int }{\mathrm{e}}^{\mathit{a}\mathit{/}\mathit{x}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{x}}^{\mathrm{4}}}\right){\mathit{x}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\exp \left[\mathrm{-}\frac{\mathrm{i}\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}{\mathit{x}}^{\mathrm{2}}}\right]\mathrm{d}\mathit{x}$(A.3)${\mathit{I}}_{\mathrm{ext}}\mathrm{=}\mathrm{\int }{\mathrm{e}}^{\mathit{a}\mathit{/}\mathit{x}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\delta }}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{x}}^{\mathrm{4}}}\right){\mathit{x}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\exp \left[\frac{\mathrm{i}\mathit{k}{\mathit{\delta }}_{\mathrm{2}}}{\mathrm{2}{\mathit{x}}^{\mathrm{2}}}\right]\mathrm{d}\mathit{x}$(A.4)where a = 1 − 2 / γ.

The two integrals cannot be evaluated exactly, however approximate solutions can be found taking into account that x ≪ 1. For illustration we will approximate the first integral while the second one can be computed in analogous manner.

Let us write the first integral as ${\mathit{I}}_{\mathrm{int}}\mathrm{=}\mathrm{\int }\frac{{\mathrm{e}}^{\mathit{a}\mathit{/}\mathit{x}}}{{\mathit{x}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{x}}^{\mathrm{4}}}\right)\left[\cos \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}{\mathit{x}}^{\mathrm{2}}}\right)\mathrm{-}\mathrm{i}\sin \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}{\mathit{x}}^{\mathrm{2}}}\right)\right]\mathrm{d}\mathit{x,}$(A.5)which can be split into two expressions according to the terms in the last square bracket, so that ${\mathit{I}}_{\mathrm{int}}^{\mathrm{1}}\mathrm{=}\mathrm{\int }\frac{{\mathrm{e}}^{\mathit{a}\mathit{/}\mathit{x}}}{{\mathit{x}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{x}}^{\mathrm{4}}}\right)\cos \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}{\mathit{x}}^{\mathrm{2}}}\right)\mathrm{d}\mathit{x}\mathit{.}$(A.6)Now let us introduce a new variable, u, so that u = 1 / x. Since x ≪ 1, it would mean that u is very large. In the new variable the last integral becomes ${\mathit{I}}_{\mathrm{int}}^{\mathrm{1}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\int }\frac{{\mathrm{e}}^{\mathrm{au}}}{{\mathit{u}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{u}}^{\mathrm{4}}\mathrm{)}\cos \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}{\mathit{u}}^{\mathrm{2}}\right)\mathrm{d}\mathit{u}\mathit{.}$(A.7)We introduce a new variable so that which means that the integral (A.7) becomes $\begin{array}{ccc}{\mathit{I}}_{\mathrm{int}}^{\mathrm{1}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\int }\frac{{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}}{{\mathit{y}}^{\mathrm{5}\mathit{/}\mathrm{4}}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{y}}^{\mathrm{2}}\mathrm{)}\cos \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\mathrm{d}\mathit{y}\\ & \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\int }\frac{{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}}{{\mathit{y}}^{\mathrm{5}\mathit{/}\mathrm{4}}}\cos \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\mathrm{d}\mathit{y}\mathrm{+}\frac{{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{\int }{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}{\mathit{y}}^{\mathrm{3}\mathit{/}\mathrm{4}}\cos \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\mathrm{d}\mathit{y}\\ & \mathrm{=}& \end{array}$(A.8)Now we use integration by parts to approximate the first integral and write $\begin{array}{ccc}{\mathit{I}}_{\mathrm{int}}^{\mathrm{11}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}\frac{{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}}{{\mathit{y}}^{\mathrm{5}\mathit{/}\mathrm{4}}}\sin \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\\ & & \end{array}$(A.9)Since the variable y is very large, the first term in the above square bracket dominates over the second term. Applying once more integration by parts we have

${\mathit{I}}_{\mathrm{int}}^{\mathrm{11}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}\frac{{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}}{{\mathit{y}}^{\mathrm{5}\mathit{/}\mathrm{4}}}\sin \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\mathrm{-}\frac{\mathit{a}}{{\mathit{k}}^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}}\frac{{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}}{{\mathit{y}}^{\mathrm{7}\mathit{/}\mathrm{4}}}\cos \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)$$\mathrm{+}\frac{\mathit{a}}{{\mathit{k}}^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{\int }\left[\frac{\mathit{a}}{\mathrm{2}{\mathit{y}}^{\mathrm{9}\mathit{/}\mathrm{4}}}\mathrm{-}\frac{\mathrm{7}}{\mathrm{4}{\mathit{y}}^{\mathrm{11}\mathit{/}\mathrm{4}}}\right]{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}\cos \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\mathrm{d}\mathit{y}\mathit{.}$(A.10)Using again the approximation of large y it is obvious that the first term of the above square bracket is larger than the second term, so we are going to neglect the second term. The integral on the right-hand side of Eq. (A.10) can be taken to the left-hand side so

$\mathrm{\int }{\mathit{a}}^{\mathit{a}\sqrt{\mathit{y}}}\cos \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\left[\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}{\mathit{y}}^{\mathrm{5}\mathit{/}\mathrm{4}}}\mathrm{-}\frac{{\mathit{a}}^{\mathrm{2}}}{{\mathit{k}}^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{y}}^{\mathrm{9}\mathit{/}\mathrm{4}}}\right]\mathrm{d}\mathit{y}\mathrm{\approx }$$\mathrm{-}\frac{\mathrm{1}}{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}\frac{{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}}{{\mathit{y}}^{\mathrm{5}\mathit{/}\mathrm{4}}}\sin \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\mathrm{-}\frac{\mathit{a}}{{\mathit{k}}^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}}\frac{{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}}{{\mathit{y}}^{\mathrm{7}\mathit{/}\mathrm{4}}}\cos \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\mathit{.}$(A.11)Now comparing the terms on the left-hand side of the above equation it is clear that the second term of the square bracket can be neglected, so the term remaining on the left-hand side of Eq. (A.11) is exactly the starting form of , therefore ${\mathit{I}}_{\mathrm{int}}^{\mathrm{11}}\mathrm{\approx }\mathrm{-}\frac{\mathrm{1}}{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}\frac{{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}}{{\mathit{y}}^{\mathrm{5}\mathit{/}\mathrm{4}}}\sin \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\mathrm{-}\frac{\mathit{a}}{{\mathit{k}}^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}}\frac{{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}}{{\mathit{y}}^{\mathrm{7}\mathit{/}\mathrm{4}}}\cos \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\mathit{.}$(A.12)For the second integral, we will need to apply the same integration by parts which leads to

${\mathit{I}}_{\mathrm{int}}^{\mathrm{12}}\mathrm{\approx }\frac{{\mathit{\delta }}_{\mathrm{1}}}{\mathit{k}}{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}{\mathit{y}}^{\mathrm{3}\mathit{/}\mathrm{4}}\sin \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\mathrm{+}\frac{\mathit{a}}{{\mathit{k}}^{\mathrm{2}}}{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}{\mathit{y}}^{\mathrm{1}\mathit{/}\mathrm{4}}\cos \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\mathrm{·}$Now turning to the second integral of Eq. (A.5) a similar procedure can be employed as above to obtain

${\mathit{I}}_{\mathrm{int}}^{\mathrm{2}}\mathrm{\approx }\frac{\mathrm{1}}{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}\frac{{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}}{{\mathit{y}}^{\mathrm{5}\mathit{/}\mathrm{4}}}\cos \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\mathrm{+}\frac{\mathit{a}}{{\mathit{k}}^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}}\frac{{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}}{{\mathit{y}}^{\mathrm{7}\mathit{/}\mathrm{4}}}\sin \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)$$\mathrm{-}\frac{{\mathit{\delta }}_{\mathrm{1}}}{\mathit{k}}{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}{\mathit{y}}^{\mathrm{3}\mathit{/}\mathrm{4}}\cos \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\mathrm{+}\frac{\mathit{a}}{{\mathit{k}}^{\mathrm{2}}}{\mathrm{e}}^{\mathit{a}\sqrt{\mathit{y}}}{\mathit{y}}^{\mathrm{1}\mathit{/}\mathrm{4}}\sin \left(\frac{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}}\mathit{y}\right)\mathit{.}$(A.13)Adding together the above integrals (in the sense of we obtain the approximate value of the integral defined by Eq. (A.3) which gives the x-dependence of the pressure in the x < x₀ region. A similar mathematical calculation is needed for the integral defined by Eq. (A.4) describing the x-dependence of the pressure in the x > x₀ region.

Returning to the original variable (x), the approximate expressions for the pressure in the two regions become

${\mathit{p}}_{\mathrm{1}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{\approx }\mathrm{𝒜}\frac{\mathrm{i\Lambda }\mathit{\omega }{\mathit{\rho }}_{\mathrm{\odot }}{\mathrm{e}}^{\mathrm{2}\mathit{/}\mathit{\gamma x}}{\mathit{c}}_{\mathrm{S}\mathrm{1}}{\mathrm{e}}^{\mathit{a}\mathit{/}\mathit{x}}{\mathit{x}}^{\mathrm{5}\mathit{/}\mathrm{2}}}{\mathrm{2}\mathit{k}{\mathit{\delta }}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{2}\mathit{/}{\mathit{x}}_{\mathrm{\odot }}}}\exp \left[\mathrm{-}\frac{\mathrm{i}\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}{\mathrm{2}{\mathit{x}}^{\mathrm{2}}}\right]$ $\mathrm{\times }\left[\mathrm{i}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{x}}^{\mathrm{4}}}\right)\mathrm{+}\frac{\mathit{ax}}{\mathit{k}{\mathit{\delta }}_{\mathrm{1}}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\delta }}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{x}}^{\mathrm{4}}}\right)\right]\mathit{,}$(A.14)

${\mathit{p}}_{\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{\approx }\mathrm{ℬ}\frac{\mathrm{i\Lambda }\mathit{\omega \alpha }{\mathit{\rho }}_{\mathit{i}}{\mathrm{e}}^{\mathrm{2}\mathit{/}\mathit{\gamma x}}{\mathit{c}}_{\mathrm{S}\mathrm{2}}{\mathrm{e}}^{\mathit{a}\mathit{/}\mathit{x}}{\mathit{x}}^{\mathrm{5}\mathit{/}\mathrm{2}}}{\mathrm{2}\mathit{k}{\mathit{\delta }}_{\mathrm{2}}{\mathrm{e}}^{\mathrm{2}\mathit{/}{\mathit{x}}_{\mathrm{0}}}}\exp \left[\frac{\mathrm{i}\mathit{k}{\mathit{\delta }}_{\mathrm{2}}}{\mathrm{2}{\mathit{x}}^{\mathrm{2}}}\right]$$\mathrm{\times }\left[\mathrm{i}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\delta }}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{x}}^{\mathrm{4}}}\right)\mathrm{-}\frac{\mathit{a}}{\mathrm{2}\mathit{k}{\mathit{\delta }}_{\mathrm{2}}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\delta }}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{x}}^{\mathrm{4}}}\right)\right]\mathrm{·}$(A.15)These expressions will be used when equating the radial components of the magnetic field perturbation (Eqs. (16) and (17)) in

order to derive the dispersion relation of waves propagating near the solar surface.

All Figures