# Numerical simulation of prominence oscillations

## 1 Introduction

Solar prominence oscillations have been the subject of both observational and theoretical papers for the past 35 years. One of the first studies (Ramsey & Smith 1966) concerned observations of global oscillations of disk filaments with periods of $`6^m`$ to $`40^m`$, that were interpreted by Hyder (1966) as predominantly vertical motions and by Kleczek & Kuperus (1969) as predominantly horizontal motions.

Tsubaki (1988) published a review on oscillation studies of limb prominences. Most of these oscillations pertain to Doppler shifts in the spectra of part of a prominence. The associated mass flows are essentially parallel to the photosphere, both longitudinal and transverse to the prominence main axis. The observed periods range from $`160^s`$ to $`82^m`$, with velocity amplitudes in the range of 0.2–3 km/s.

Zhang et al. (1991), Zhang & Engvold (1991) and Thompson & Schmieder (1991) studied disk filaments and found locally periods of $`2.5^m22^m`$, with velocity amplitudes of $`0.51.25`$ km/s. As in the observations by Ramsey & Smith these oscillations may have both horizontal and vertical components.

Since Tsubaki’s review more limb studies have been performed by Mashnich & Bashkirtsev (1990), Suematsu et al. (1990), Bashkirtsev & Mashnich (1993), Mashnich et al. (1993), Balthasar et al. (1993), Balthasar & Wiehr (1994), Park et al. (1995), Sütterlin et al. (1997) and Molowny-Horas et al. (1997). These observations all confirm and extend previous results. There is now enough evidence to suggest the following tentative classification for prominence oscillations (see also Bashkirtsev & Mashnich 1993).

* very short periods: $`P30^s`$ (Balthasar et al. 1993), perhaps due to fast waves propagating along flux tubes (Roberts et al. 1984).
* short periods: $`P3^m10^m`$, at least some of which are related to photospheric or chromospheric forcing with periods of $`3^m`$ and $`5^m`$ (Balthasar et al. 1986, Zhang et al. 1991).
* intermediate periods: $`P10^m40^m`$, which are probably genuine eigenmodes of the prominence (Ramsey & Smith 1966, Balthasar et al. 1988, Bashkirtsev & Mashnich 1993). Many prominences show nearly the same oscillation period each time they are perturbed.
* long periods: ($`P40^m114^m`$), which may be related to chromospheric forcing (Balthasar et al. 1988).

We point out that the longest data set is about 7 hours long and the best temporary resolution is a couple of seconds.

In particular the observed $`10^m40^m`$ eigenmodes are interesting as they are damped and thus loose energy, perhaps due to some interaction with the ambient corona (Ramsey & Smith 1966, Kleczek & Kuperus 1969). Typically, the quality factor $`Q=\pi T_{\mathrm{damp}}/P<6`$, with $`T_{\mathrm{damp}}`$ the damping time of the oscillation. This implies that three to four oscillations can be observed after the impulsive perturbation of the prominence. Understanding the damping mechanisms will give more insight in prominence dynamics and may yield an extra diagnostic tool for prominence and ambient coronal plasma parameters.

Prominence oscillations have been theoretically modelled by many authors. Some model the prominence as a harmonic oscillator with an ad-hoc damping term: Hyder (1966) used viscous effects, Kleczek & Kuperus (1969) used emission of sound waves, van den Oord & Kuperus (1992), Schutgens (1997a) and van den Oord et al. (1998) used emission of Alfvén waves. Other authors construct a simple MHD equilibrium and study oscillations thereof using the linear adiabatic MHD equations (Oliver et al. 1992, 1993, Oliver & Ballester 1995, 1996, Joarder & Roberts 1992ab, 1993, Joarder et al. 1997). Generally the latter approach yields marginally stable oscillations (no damping), although Joarder & Roberts (1992a) and Joarder et al. (1997) claim that leaky waves are possible solutions to their equations. Apparently, little effort has gone into studying the damping mechanisms themselves.

Schutgens (1997ab) and van den Oord et al. (1998) recently studied the global equilibrium of prominences treating the evolution of the magnetic field in a self-consistent way. The equation of motion for the prominence (approximated as a line current) was solved simultaneously with the Maxwell equations for the electro-magnetic fields. Their results show that the Alfvén travel time between prominence and photosphere $`\tau `$ is an important time scale of the system. In particular, van den Oord et al. found that it is only possible to obtain stable prominence equilibria by taking damping mechanisms into account. Hence, damping is not just necessary to describe prominence oscillations quantitatively correct, but it is an essential ingredient of a prominence equilibrium.

In this paper, we study prominence oscillations, and in particular the effect of the ambient coronal plasma on the prominence motion. We use the Versatile Advection Code (VAC) to solve numerically the isothermal MHD equations in two dimensions. VAC has been developed by G. Tóth (1996, 1997) and is capable of solving a variety of hydrodynamical and magneto-hydrodynamical problems in one, two and three dimensions using a host of numerical methods. In Sect. 2 a simple analytical model for prominence equilibrium and dynamics is discussed that will be used for comparison with our numerical results. In section 3 we describe the isothermal equations that were solved numerically, the methods used, the grid structure and the initial conditions. In Sect. 4 the simulations are described. Using the analytical model mentioned before, these simulations are interpreted in terms of the physical processes involved. A summary and the conclusions can be found in Sect. 5. All dimensional variables in this paper are in rational MKSA units, unless stated differently.

## 2 The line current approximation

We briefly recapitulate the line current approximation for prominence equilibrium and dynamics. This approximation will serve as a guideline for discussing our numerical results. The prominence is approximated by a straight, infinitely thin and long, line current $`I_0>0`$ at a height $`y_0`$ parallel to the photosphere. Along the prominence we assume invariance. The effect of the massive photosphere on the prominence magnetic field is modelled through a mirror current $`I_0`$ at depth $`y_0`$ below the surface of the Sun (Kuperus & Raadu 1974, van Tend & Kuperus 1978, Kaastra 1985, Schutgens 1997a). The momentum equations governing the global prominence dynamics are (ignoring gravity) 

$`\sigma \ddot{x}`$ $`=`$ $`I_0\left[B_{\mathrm{mir}}^y(x,y)+B_{\mathrm{cor}}^y(x,y)\right]\nu _x\dot{x},`$   
$`\sigma \ddot{y}`$ $`=`$ $`I_0\left[B_{\mathrm{mir}}^x(x,y)+B_{\mathrm{cor}}^x(x,y)\right]\nu _y\dot{y}.`$ (1)   
where $`\sigma `$ is the longitudinal mass density of the oscillating structure, $`𝑩_{\mathrm{cor}}`$ is the coronal arcade field in which the prominence is located and $`𝑩_{\mathrm{mir}}`$ is the field due to the mirror current.

The interaction between the moving filament and the ambient coronal plasma gives rise to viscous effects (Hyder 1966) and emission of magneto-acoustic waves (Kleczek & Kuperus 1967, van den Oord & Kuperus 1992) that act as damping mechanisms. These are heuristically modelled through $`\nu _x`$ and $`\nu _y`$, damping constants that are in reality determined by the flow field around the prominence. Since viscosity of the coronal plasma is negligible, we concentrate on the emission of magneto-acoustic waves. An approximation for the damping constants can be found by considering linear motion with constant velocity $`v`$ of a solid body through a homogeneous plasma. If the cross section of the body perpendicular to its motion is $`A`$, the body transfers $`2\rho _{\mathrm{cor}}cAv`$ momentum per unit time onto the plasma (Landau & Lifschitz 1989, p. 256). Here $`c`$ is the characteristic wave speed of the plasma with density $`\rho _{\mathrm{cor}}`$. A factor $`2`$ is added since both the front- and backside of the object transfer momentum. Hence, the damping constants have the form

$$\nu =2\rho _{\mathrm{cor}}cA.$$
(2)   
The coronal arcade is generated by a magnetic line dipole $`M_\mathrm{d}>0`$, a depth $`H_\mathrm{d}`$ below the photosphere

$`B_{\mathrm{cor}}^x(x,y)`$ $`=`$ $`{\displaystyle \frac{\mu _0M_\mathrm{d}}{\pi }}{\displaystyle \frac{x^2(y+H_\mathrm{d})^2}{\left(x^2+(y+H_\mathrm{d})^2\right)^2}},`$   
$`B_{\mathrm{cor}}^y(x,y)`$ $`=`$ $`{\displaystyle \frac{2\mu _0M_\mathrm{d}}{\pi }}{\displaystyle \frac{x(y+H_\mathrm{d})}{\left(x^2+(y+H_\mathrm{d})^2\right)^2}}.`$ (3)   
The mirror current’s field at the location of the filament, in the quasi-stationary approach, is given by 

$`B_{\mathrm{mir}}^x(x,y)`$ $`=`$ $`{\displaystyle \frac{\mu _0}{4\pi }}{\displaystyle \frac{I_0}{y}},`$   
$`B_{\mathrm{mir}}^y(x,y)`$ $`=`$ $`0.`$ (4)   
The fact that $`B_{\mathrm{mir}}^y=0`$ is a direct consequence of the mirror current mirroring the motion of the filament (to conserve the photospheric flux) and the quasi-stationary field assumption. When this assumption is dropped and the magnetic fields evolve dynamically according to Maxwell’s laws, the expressions for $`𝑩_{\mathrm{mir}}`$ become far more complicated and in particular $`B_{\mathrm{mir}}^y0`$ (Schutgens 1997a, van den Oord et al. 1998).

Assuming quasi-stationary field evolution ($`v_\mathrm{A}\mathrm{}`$), prominences are in stable equilibrium provided they are on the symmetry axis of the arcade ($`y=0`$) and at a height $`y_0<H_\mathrm{d}`$. Furthermore, the current should attain the value

$$I_0=\frac{4y_0M_\mathrm{d}}{(y_0+H_\mathrm{d})^2}.$$
(5)   
Note that this prominence equilibrium has an inverse polarity topology (see also Fig. 5). In fact, it corresponds to a Kuperus-Raadu prominence (see van Tend & Kuperus 1978).

If one linearizes around this equilibrium, the equations of vertical and horizontal motion decouple (due to the symmetry of the coronal field) and both have the form of the familiar damped harmonic oscillator. The solutions are oscillations, characterized by frequencies $`\omega `$ and damping rates $`\delta `$:

$`\omega _x`$ $`=`$ $`\left(\mathrm{\Omega }_x^2{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{\nu _x}{\sigma }}\right)^2\right)^{\frac{1}{2}}\mathrm{\Omega }_x^2={\displaystyle \frac{8\mu _0M_\mathrm{d}^2}{\pi \sigma }}{\displaystyle \frac{y_0}{(y_0+H_\mathrm{d})^5}},`$ (6)   
$`\delta _x`$ $`=`$ $`{\displaystyle \frac{\nu _x}{2\sigma }},`$ (7)   
$`\omega _y`$ $`=`$ $`\left(\mathrm{\Omega }_y^2{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{\nu _y}{\sigma }}\right)^2\right)^{\frac{1}{2}}\mathrm{\Omega }_y^2={\displaystyle \frac{4\mu _0M_\mathrm{d}^2}{\pi \sigma }}{\displaystyle \frac{H_\mathrm{d}y_0}{(y_0+H_\mathrm{d})^5}},`$ (8)   
$`\delta _y`$ $`=`$ $`{\displaystyle \frac{\nu _y}{2\sigma }}.`$ (9)   
where $`\mathrm{\Omega }_x`$ and $`\mathrm{\Omega }_y`$ are the quasi-stationary frequencies of the system without damping.

If one drops the quasi-stationary assumption, the same equilibrium is found, but its stability then also depends on the value of the coronal Alfvén speed . In general, the solutions to the equations, which are again decoupled, are damped or growing (!) harmonic oscillations (Schutgens 1997ab, van den Oord et al. 1998).

## 3 Numerical methods

We solve the time-dependent ideal isothermal MHD equations in two dimensions using the Versatile Advection Code (VAC) developed by one of us (G. Tóth). Descriptions of this numerical code can be found in Tóth (1996, 1997). The photosphere coincides with the $`y=0`$ plane. In conservative form the equations for density $`\rho `$, mass flux $`\rho 𝒗`$ and magnetic field $`𝑩`$ are ($`p`$ is the thermal pressure)

$`_t\rho `$ $`+\left(𝒗\rho \right)`$ $`=0,`$   
$`_t(\rho 𝒗)`$ $`+\left(𝒗\rho 𝒗𝑩𝑩\right)+(p+B^2/2)`$ $`=0,`$   
$`_t𝑩`$ $`+\left(𝒗𝑩𝑩𝒗\right)`$ $`=0.`$ (10)   
Note that we ignore gravity (see Sect. 5 for a discussion). These equations must be solved together with an equation of state $`p=c_\mathrm{s}^2\rho `$ and the condition that $`𝑩=0`$. Here $`c_\mathrm{s}`$ is the sound speed, a free parameter of the system, which is constant throughout the numerical domain. The magnetic field unit is chosen such that the current density satisfies $`𝑱=\times 𝑩`$ (i.e. $`\mu _0=1`$), all other units are rational MKSA.

The equations are discretized on the same grid and solved using a FCT (Flux Corrected Transport) scheme. Since there are no discontinuities in the solution, FCT performs well. Similar results can be obtained using a TVD (Total Variation Diminishing) Lax-Friedrichs scheme but, for the problem at hand, FCT is more efficient. For a comparison of different methods see Tóth & Odstrčil (1996). A projection scheme (Brackbill & Barnes 1980) is used to keep the magnetic divergence small. Typically $`|𝑩|<10^3B/L`$, where $`B`$ and $`L`$ are characteristic strength and length scale of the magnetic field.

We use a Cartesian grid of $`190\times 155`$ cells ($`3\times \mathrm{\hspace{0.17em}10}^9`$ m by $`1.5\times \mathrm{\hspace{0.17em}10}^9`$ m) that is strongly distorted, with the highest resolution at the location of the prominence and a much smaller resolution near the coronal edges of the computational domain (see Fig. 1). In a region surrounding the actual prominence the grid spacing is constant. Typically, the grid spacing increases with 10% from cell to cell outside this region. As a consequence the difference in resolution at the prominence and at the far coronal edges can amount to a factor 100 (see Fig. 2). Near the boundaries the grid spacing is again kept constant.

The boundary conditions are implemented using two layers of ghost cells around the physical part of the grid. The actual boundary is located between the inner ghost cells and the outermost cells of the physical grid. The solution on the physical part of the grid can be advanced using fluxes calculated from the ghost and physical cells around the boundary. The prescription for the ghost cells depends on the specific boundary condition and may be constant or depend on the solution in the adjacent physical cells.

The photosphere is much denser than the corona and is therefore strongly reflecting: there should be no mass flux or energy flux across it. We therefore choose the plasma density symmetric around the photospheric boundary ($`\rho ^{\mathrm{ghost}}=\rho ^{\mathrm{physical}}`$) and vertical momentum anti-symmetric ($`\rho v_y^{\mathrm{ghost}}=\rho v_y^{\mathrm{physical}}`$). The other momentum component $`\rho v_x`$ is chosen anti-symmetric as well, since we assume no flows along the photosphere (‘no slip’ condition). The magnetic field is tied to the dense photosphere and the photospheric flux ($`B_y`$) is conserved. Magnetic waves should be reflected. The initial field solution is stored in computer memory and subtracted from the advanced solution. The ‘linearized’ field solution thus obtained is chosen to be symmetric or anti-symmetric in the photospheric boundary: $`B_x^{\mathrm{ghost}}=B_x^{\mathrm{ghost},0}+(B_x^{\mathrm{physical}}B_x^{\mathrm{physical},0})`$ and $`B_y^{\mathrm{ghost}}=B_y^{\mathrm{ghost},0}(B_y^{\mathrm{physical}}B_y^{\mathrm{physical},0})`$. In this way, one obtains a good reflection of waves for small perturbations.

The other (coronal) boundaries do not coincide with a physical boundary and should be as open as possible. Since any choice of boundary condition will always generate some reflection, which we want to avoid, we decided to place these boundaries at large distances from the prominence so that twice the wave crossing time (prominence–coronal boundary) is longer than the simulation time. Also, the coarsening of the grid damps the outward moving waves and the reflection is minimized. At the location of the coronal boundaries we prescribe fixed $`𝑩^{\mathrm{ghost}}`$, but copy $`\rho `$ and $`\rho 𝒗`$ from the physical part of the grid to the ghost cells.

The initial configuration is a superposition of three different analytical MHD equilibria. The global coronal structure is a potential arcade, given by Eq. (3). Since we ignore gravity, the plasma density in this arcade is constant. To this equilibrium we add the fields and plasma densities of two current carrying flux tubes, one above, the other below the photosphere. Both flux tubes are located on the polarity inversion line of the arcade. The flux tube at a height $`y_0`$ above the photosphere represents the prominence. Its total axial current is $`I_0`$. The flux tube $`y_0`$ below the photosphere represents the mirror current and has a total current $`I_0`$.

The flux tube equilibrium is derived starting from a current profile (see also Forbes 1990)

$$j_z(r)=\{\begin{array}{cc}j_0\hfill & rr_0\mathrm{\Delta }r_0,\\ j_0\mathrm{sin}^2\left(\frac{\pi }{2}\frac{rr_0}{\mathrm{\Delta }r_0}\right)\hfill & r_0\mathrm{\Delta }r_0<rr_0,\\ 0\hfill & r>r_0.\end{array}$$
(11)   
The radius of the tube $`r_0`$ is larger than $`\mathrm{\Delta }r_0`$, the size of the region in which the current density drops off to zero. The total axial current in the flux tube is

$$I_0=2\pi _0^{r_0}r^{}j_z(r^{})𝑑r^{}$$
(12)   
and the azimuthal field component follows from Stokes’ theorem

$$B_\varphi (r)=\{\begin{array}{cc}\frac{\mu _0}{r}_0^rr^{}j(r^{})𝑑r^{}\hfill & rr_0,\\ \frac{\mu _0I_0}{2\pi r}\hfill & r>r_0.\end{array}$$
(13)   
The current is confined to $`r<r_0`$. Beyond this range, the field is potential. The pinching field $`B_\varphi `$ has to be balanced by gas pressure

$$p(r)=_{r_0}^rB_\varphi (r^{})j(r)𝑑r,$$
(14)   
where it is assumed that the flux tube proper carries no mass outside $`r>r_0`$. The longitudinal mass density of the flux tube equals

$$\sigma _0=2\pi _0^{r_0}r^{}p(r^{})/c_\mathrm{s}^2𝑑r^{},$$
(15)   
while the total longitudinal mass density of the prominence is (due to the superposition of equilibria)

$$\sigma =\sigma _0+\pi \rho _{\mathrm{cor}}r_0^2.$$
(16)   
In a global sense, equilibrium is obtained when Eq. (5) is satisfied. However, especially near the outer edge of the prominence no force balance exists. It is necessary to let the initial configuration relax to a numerical equilibrium, before studying oscillations. As it turns out, force imbalance is mostly due to the gradients in density and field within the flux tube (i.e. numerical discretization errors) and a stable equilibrium is readily obtained. To prevent the flows at the outer edge of the prominence from becoming too large during the relaxation phase, and cause an instability, an artificial damping term $`_t\rho 𝒗=\alpha \rho 𝒗`$ is used. This term is switched off during subsequent oscillation studies.

Once a numerical equilibrium has been found, we perturb it and study the resulting oscillations. In all cases, the perturbation is caused by instantaneously adding momentum to the prominence.

To check the validity of our results, we performed oscillation simulations for the same physical parameters, but different numerical parameters. In particular we changed the grid resolution ($`\mathrm{\Delta }x=5\times 10^5`$ or $`10^6`$ m inside the prominence), the amount of stretching of the grid (0%, 5% and 10%), the duration of the relaxation ($`T=\mathrm{9\hspace{0.17em}000}`$ or $`T=\mathrm{20\hspace{0.17em}000}`$ s) and the constraint on $`|𝑩|`$ ($`<10^4,10^3B/L`$) and found that the results are similar. Using $`0\%`$ stretching implies that the grid is rather small in its physical size, due to computational limitations. Reflection at the coronal boundaries will then limit the usefulness of the simulations to the first $`4000`$ s for $`\rho _{\mathrm{cor}}=10^{12}`$ kg/m<sup>3</sup>.

In addition, the effect of the photospheric boundary condition was studied. A perfectly reflecting photosphere does not allow a net momentum or energy flux. Also, the magnetic flux distribution is constant. The implementation of the photospheric boundary conditions automatically ensures zero momentum flux and constant magnetic flux. However, since the magnetic field is not symmetric in $`y=0`$, a small Poynting flux is present. Consider a box-shaped surface around the prominence. The coronal ‘walls’ of this surface each have a shortest distance to the prominence center of $`2.5\times 10^7`$ m (see Fig. 6) and reach down to the photosphere, along which we close the box. The total time-averaged Poynting flux through the photosphere is typically $`10^310^2`$ times smaller than the coronal flux, for vertical oscillations. For horizontal oscillations, the coronal flux is small in essence (due to the anti-symmetry) and provides no reasonable yardstick (of course one could still use the flux through a single ‘wall’, but this is bound to lead to similar results as for vertical oscillations). Furthermore, we point out that the equilibrium height of the prominence during oscillation studies does not change by more than 0.06%, and often even less. We therefore feel that the photospheric boundary is sufficiently well represented in our numerical boundary conditions.

## 4 Oscillation studies

Before we study oscillations in prominences, we first have to compute a stable, numerical prominence equilibrium.

For the coronal arcade, we choose the following parameter values: $`H_\mathrm{d}=4\times \mathrm{\hspace{0.17em}10}^8`$ m, $`M_\mathrm{d}=10^{20}`$ Am and plasma-density $`\rho _{\mathrm{cor}}=10^{12}`$ kg/m<sup>3</sup>. The two flux tubes are given by: $`r_0=1.5\times \mathrm{\hspace{0.17em}10}^7`$ m, $`\mathrm{\Delta }r_0=10^7`$ m and $`j_0=2\times \mathrm{\hspace{0.17em}10}^4`$ A/m<sup>2</sup>. Hence the longitudinal mass density of the flux tube proper is $`\sigma _0=1.28\times \mathrm{\hspace{0.17em}10}^4`$ kg/m, and the density of the prominence is $`\sigma =1.35\times \mathrm{\hspace{0.17em}10}^4`$ kg/m. The total current is $`I_0=6.5\times \mathrm{\hspace{0.17em}10}^{10}`$ A. Such a prominence is expected to oscillate with a horizontal period of 2777 s and a vertical period of 1118 s. The sound speed was chosen $`c_\mathrm{s}=128.5`$ km/s, typical of the corona at a temperature $`T=10^6`$ K.

We now let the initial state as defined in the previous section relax to a numerical equilibrium. The artificial damping term ($`_t\rho 𝒗=\alpha \rho 𝒗`$) was used to prevent numerical instabilities. We choose $`\alpha =0.1`$ for the first $`10^4`$ seconds of the relaxation, and $`\alpha =0.01`$ during the latter $`10^4`$ seconds. Figure 3 shows the relaxation of the volume averaged horizontal and vertical momentum. At the end of the relaxation the flows in the larger part of the arcade are typically less than 5 m/s. Along the boundary of the prominence body a rather irregular flow field exists with velocities of 175 m/s at most. The prominence still moves downwards at a systematic speed of some 1 m/s (see Fig. 4). Considering both the magnitude of the velocity perturbation applied later and the total simulation time scale, we consider these residual flows to be unimportant. The agreement between simulations starting from different relaxations (with a duration of either $`\mathrm{9\hspace{0.17em}000}`$ or $`\mathrm{20\hspace{0.17em}000}`$ s) confirm this.

The numerical equilibrium does not deviate much from the initial state. The resulting field topology around the prominence is shown in Fig. 5. The prominence itself is a region of high plasma-$`\beta `$, surrounded by a region of magnetically dominated plasma (see Fig. 6). At large distances from the prominence the plasma-$`\beta `$ becomes larger than unity again, as the magnetic field strength decreases while the coronal density is constant (due to the absence of gravity). The Alfvén speed is also shown in Fig. 6. It increases as one approaches the prominence from the corona, once inside it falls off rapidly. The sound speed, of course, is the same everywhere (it is a free parameter of Eq. (10)).

The obtained numerical equilibrium is used to derive a series of numerical equilibria by adding (or subtracting) a constant value from the density in each cell. Since this does not create any additional forces (the force due to gas pressure is $`c_\mathrm{s}^2\rho `$) a new equilibrium is found. In the resulting equilibria we disturb the inner part of the prominence (all plasma within $`1.25\times \mathrm{\hspace{0.17em}10}^7`$ m from the center of mass) with a velocity perturbation of 10 km/s, at an angle of $`45^o`$ to the photosphere. The evolution of the system is studied for $`10^4`$ s, in some cases even for $`2\times \mathrm{\hspace{0.17em}10}^4`$ s.

In Fig. 7 the evolution of the volume averaged momentum is shown for all eight cases considered. Initially momentum is concentrated in the prominence, but it is redistributed throughout the corona in the subsequent evolution. A global oscillation is apparent, whose properties depend strongly on coronal density.

For one particular case, the flow field is shown in Figs. 8 and 9. Although only the velocity is shown, the prominence stands out clearly in most graphs. It seems to move through the coronal plasma as a rigid body. The largest velocities are usually found outside the prominence. The coronal plasma ‘washes around the filament’.

In particular, we studied the motion of the prominence. This was done by computing every so many time steps the center of longitudinal mass and current density. To rule out the contribution of the coronal part of the grid, only densities above a certain threshold where used. As the contrast between typical coronal and prominence current densities is larger than the contrast between typical coronal and prominence mass densities, the first provide a better estimate of the location of the prominence. Nevertheless, the results are always similar.

In essence, the motion of the filament can be described by two decoupled damped harmonic oscillators. The horizontal resp. vertical displacement of the center of longitudinal mass or current density of the prominence may be fitted to $`A\mathrm{e}^{t/T_{\mathrm{damp}}}\mathrm{sin}2\pi t/P`$. The resulting periods and damping times for the horizontal and vertical oscillations are listed in Table 1. Sometimes strong transient effects at the start of the simulation and noise at later times (when the velocities are smaller) cause deviations from a simple damped harmonic oscillator. The horizontal oscillations yield, in general, better fits. Typically the low density simulations lead to better fits than the high density simulations. For horizontal motions we used the cases 1–6, for vertical motions we used the cases 1–5 (see Table1). No reliable vertical time scales could be obtained for the last two cases due to strong transients and noise.

Comparing the results with simulations where a purely vertical or horizontal perturbation was applied shows that the horizontal and vertical motions of the filament are actually decoupled. This is understandable as the filament is located on the symmetry line of the coronal field (see also van den Oord et al. 1998). At the same time this suggests that non-linear effects are not important. Further evidence for the linearity of the numerical results is found in the good fits of the prominence oscillation curves to a damped harmonic oscillator. However, simulations with smaller or larger perturbations (3 km/s or 20 km/s instead of 10 km/s) yield different damping times for the horizontal motions (see Fig. 10). For all three simulations the horizontal periods differ by only $`0.1\%`$, while the horizontal damping times differ by $`3060\%`$! Apparently, the horizontal damping mechanism depends non-linearly on the flow speeds. We have no explanation for this result, but believe it to be genuine. Extensive tests with higher spatial and temporal resolution argue strongly against the possibility of a numerical effect. It is a strong indication that the damping mechanism for horizontal and vertical motion differ, at least in their dependencies on the flow field.

The horizontal periods show a strong dependence on the coronal density, while the vertical periods are almost constant (Fig. 11). The quality factors (Fig. 12) of the horizontal oscillation are typically larger than four, and hence damping contributes little to the oscillation frequency (see Eq. (6)). The dotted lines (in Fig. 11) represent the undamped quasi-stationary period for a prominence with a longitudinal mass density of $`\sigma =\sigma _0+\pi r_0^2\rho _{\mathrm{cor}}`$. Clearly this is a bad fit to the numerical results for horizontal oscillations. Let us assume that, due to the magnetic structure, the actual oscillating body has a radius $`r_{\mathrm{eff}}>r_0`$. The actual oscillating body now has an effective longitudinal mass density $`\sigma _{\mathrm{eff}}=\sigma _0+\pi r_{\mathrm{eff}}^2\rho _{\mathrm{cor}}`$. From Eq. (6) we obtain

$$\frac{1}{\rho _{\mathrm{cor}}}\left(\frac{2\mu _0M_\mathrm{d}^2}{\pi ^3}\frac{y_0}{(y_0+H_\mathrm{d})^5}P_x^2\sigma _0\right)=\pi r_{\mathrm{eff}}^2$$
(17)   
the left-hand-side of which can be fitted to a function of the form $`C\rho _{\mathrm{cor}}^\gamma `$. We find that $`r_{\mathrm{eff}}=5.0\times \mathrm{\hspace{0.17em}10}^7(\rho _{\mathrm{cor}}/10^{12})^{0.04}`$ m, weakly dependent on the coronal density. The typical radius $`r_{\mathrm{eff}}=5.0\times \mathrm{\hspace{0.17em}10}^7`$ m agrees well with the location of the magnetic surface across which the connectivity of the field lines changes from closed field lines in the corona to field lines anchored in the photosphere. An even better approximation to this surface is obtained by taking into account the effect of induced mass as described by Landau & Lifschitz (1989, p. 29; see also Lamb 1945). For an incompressible, potential hydrodynamical flow the effective mass is the sum of the mass of the actual oscillating body ($`\sigma _0+\pi r_0^2\rho _{\mathrm{cor}}`$) plus the mass of the fluid displaced by the body. In that case $`\sigma _{\mathrm{eff}}=\sigma _0+2\pi r_{\mathrm{eff}}^2\rho _{\mathrm{cor}}`$ (note the factor 2!) and we find $`r_{\mathrm{eff}}=3.5\times \mathrm{\hspace{0.17em}10}^7(\rho _{\mathrm{cor}}/10^{12})^{0.04}`$ m.

In Fig. 13, we have plotted the ratio of vertical damping time to horizontal damping time. From Eqs. (2), (7) and (9) this ratio is proportional to $`c_x/c_y`$. Here $`c_x`$ is the typical wave speed for perturbations travelling parallel to the photosphere and $`c_y`$ is the typical wave speed for perturbations travelling perpendicular to the photosphere. For simplicity, we assume that the cross sections for both directions have the same $`\rho _{\mathrm{cor}}`$ dependence, but given the azimuthal invariance of the prominence flux tube this is probably a fair approximation. Now $`c_x`$ and $`c_y`$ are equal to either the slow cusp speed $`c_\mathrm{T}`$, the Alfvén speed $`c_\mathrm{A}`$ or the fast speed $`c_\mathrm{f}`$. In a low plasma-$`\beta `$ environment, the cusp speed equals the sound speed $`c_\mathrm{s}`$, while the fast speed equals the Alfvén speed $`c_\mathrm{A}`$. Hence $`c_x/c_y\rho _{\mathrm{cor}}^{0.5},1`$ or $`\rho _{\mathrm{cor}}^{0.5}`$ depending on whether slow, Alfvén or fast wave emission prevails in a certain direction. From the data, we find $`c_x/c_y\rho _{\mathrm{cor}}^{0.53}`$ which strongly suggests that slow wave emission damps horizontal motions, while Alfvén or fast wave emission damps vertical motions. Since the prominence moves vertically to the field lines in the latter case, fast waves are more likely.

The conclusions regarding damping mechanisms can be substantiated further by comparing the relevant time scales for horizontal and vertical oscillations independently to the Eqs. (7) and (9) (see Fig. 14).

For the mass, we use the effective longitudinal mass density $`\sigma _{\mathrm{eff}}`$ as obtained previously. From Eq. (2) and either Eq. (7) or (9) we find

$$\frac{\sigma _{\mathrm{eff}}}{\rho _{\mathrm{cor}}T_{\mathrm{damp}}}=Ac,$$
(18)   
the left-hand-side of which can be fitted to $`C\rho _{\mathrm{cor}}^\gamma `$. Here $`A`$ is the cross-section perpendicular to the direction of motion. For horizontal oscillations we find $`c_x\rho _{\mathrm{cor}}^{0.06\pm 0.06}`$, and $`A2\times \mathrm{\hspace{0.17em}10}^7`$ m<sup>2</sup>. This suggests excitation of slow waves by horizontal motion of the prominence. For vertical oscillations $`c_y\rho _{\mathrm{cor}}^{0.54\pm 0.04}`$ which suggests excitation of fast waves by vertical motion of the prominence. Since the Alfvén speed varies in space, it is not possible to derive a value of $`A`$ from the product of $`Ac`$.

The cross-section $`A`$ for the horizontal damping mechanism yields a dimension for the oscillating body smaller than the extended size $`r_{\mathrm{eff}}`$ as obtained from the horizontal periods. This may be due to field line curvature. Since the size of the oscillating ‘solid’ body is determined by the field topology, Lorentz forces cause the waves that carry away momentum. However, Lorentz forces only act perpendicular to, not parallel to the field lines. Thus only perturbations of magnetic field with a strong vertical component can excite horizontally travelling waves.

## 5 Summary and conclusions

We have made a numerical investigation of prominence oscillations, by solving the isothermal MHD equations in two dimensions. First we computed a prominence equilibrium that is very similar to the Kuperus-Raadu (inverse polarity) topology. However, in our numerical equilibrium the prominence is not infinitely thin, but instead well resolved. From this equilibrium, we derived other equilibria with different coronal plasma densities. We then perturbed the system by instantaneously adding momentum to the prominence mass and followed the ensuing oscillations. The dependence of the characteristic time scales (periods and damping times) on the coronal plasma density was analyzed in terms of a solid body moving through a fluid. To our knowledge this is the first attempt at numerically simulating prominence oscillations.

In our numerical model, we ignored the effects of gravity and thermodynamics, for the sake of clarity and practicality. Also, up to now the numerical boundary conditions we are using do not seem to allow a stable gravitationally stratified corona.

We believe, however, that gravity and thermodynamics do not contribute significantly to the physics of the system. Gravity is of small consequence for the equilibrium of a Kuperus-Raadu prominence as detailed by van Tend & Kuperus (1978). Force balance is due to two Lorentz forces: one due to the coronal magnetic arcade, the other due to the photospheric flux conservation. The gravitational force can be ignored when describing global equilibrium. Furthermore, the scale height of a corona of $`T=10^6`$ K is $`3\times \mathrm{\hspace{0.17em}10}^8`$ m, which is larger than the typical vertical size of prominences. The inclusion of gravity would change the appearance of the prominence into a slab, with prominennce matter accumulating in the pre-existing (!) dips in the magnetic field lines. The magnetic field configuration would hardly change.

Likewise, the absence of a thermally structured corona and prominence seems of only minor influence. The most important effect would be that a cool prominence will be heavier than a prominence of coronal temperatures. In our model this is balanced by the absence of a longitudinal field. As a consequence, the pinching effect of the azimuthal prominence field has to be balanced by gas pressure only (in real prominences the longitudinal field pressure contributes significantly). The total longitudinal mass density of our prominence ($`1.3\times 10^4`$ kg/m<sup>3</sup>) agrees well with that of a slab of height $`5\times \mathrm{\hspace{0.17em}10}^7`$ m, width $`6\times \mathrm{\hspace{0.17em}10}^6`$ m and density $`5\times \mathrm{\hspace{0.17em}10}^{11}`$ kg/m<sup>3</sup>. Oliver & Ballester (1996) studied the influence of the prominence-corona transition region (characterized by strong temperature gradients) and found that it mainly influences the prominence internal oscillations, but not the global oscillations.

We point out that the inclusion of both gravity and thermodynamics would give the prominence proper the appearance of a cool slab, suspended in the dips of the field lines belonging to the prominence current. Presumably this will not change the global oscillation discussed in this paper, since that mode is determined by the overall field structure.

The results indicate that, for typical coronal densities ($`\rho _{\mathrm{cor}}=10^{13}10^{12}`$ kg/m<sup>3</sup>), the prominence structure can indeed be viewed as a solid body moving through a fluid. However, the mass of this solid body is determined by the magnetic topology, not the prominence proper. In particular, the mass of the body seems to be determined by coronal field lines that enclose the prominence proper. In a low plasma-$`\beta `$ environment this is to be expected. As a consequence, the total mass of the oscillating solid body is larger than the mass of the prominence proper.

Due to the symmetry of the system, horizontal and vertical prominence oscillations decouple. These oscillations can each be interpreted as the motion of a damped harmonic oscillator. The horizontal periods and the horizontal and vertical damping times can be explained by assuming that the actual oscillating structure is larger than the prominence proper, due to the magnetic field topology. Only the vertical periods do not agree with this model. They are nearly constant for different values of the coronal density and are best modelled by the oscillation of the prominence proper in a corona with vanishing plasma density ($`\rho _{\mathrm{cor}}0`$). For realistic coronal densities, the horizontal periods change with $`315\%`$ at most when taking the actual oscillating structure into account. However, for lower prominence longitudinal mass densities $`\sigma `$ in a stronger coronal background field (and hence larger prominence currents $`I_0`$), the effect will be much more pronounced.

Vertical oscillations lead to the emission of fast waves that carry momentum away from the prominence and damp the oscillation. Horizontal oscillations, on the other hand, lead to the emission of slow waves. These will damp the horizontal oscillation, but less effectively than fast waves ($`Q_x>Q_y`$).

The difference in wave emission between horizontal and vertical oscillations can be understood in terms of the coronal arcade in which the prominence is embedded. Due to the large scale height, the arcade field is close to horizontal near the prominence. Waves that travel in the vertical direction (up or down) therefor travel more or less perpendicular to the field lines and must be fast waves. Waves that travel in the horizontal direction travel along the field lines. They could be either Alfvén waves or magneto-acoustic waves. As regards excitation of the waves, it is obvious that the vertical motion of a prominence across field lines that are nearly horizontal will compress both gas and magnetic field and thus set off fast waves. The excitation of the magneto-acoustic slow waves (for our analysis strongly suggests they are slow waves) is not as well understood. But apparently the prominence acts as a piston during horizontal motions along magnetic field lines of the arcade.

We surmise that for smaller scale height of the arcade the slow waves might be replaced by fast waves.

* N.A.J. Schutgens was financially supported by the Netherlands Organisation for Scientific Research (NWO) under grant nr. 781-71-047. He gratefully acknowledges stimulating discussions with Max Kuperus and Bert van den Oord. The Versatile Advection Code (VAC) was developed by G. Tóth as part of the project on ‘Parallel Computational Magneto-Fluid Dynamics’, funded by the Netherlands Organisation for Scientific Research (NWO) Priority Program on Massively Parallel Computing, while he was working at the Astronomical Institute of Utrecht University. G. Tóth currently receives a post-doctoral fellowship (D25519) from the Hungarian Science Foundation (OTKA).