- Identification of gravity waves in hydrodynamical simulations
- 1 Introduction
- 2 The hydrodynamical model
- 3 Measuring IGWs excited by the oscillating bubble from classical methods
- 4 A new method using the anelastic subspace
- 5 Results
- 6 Conclusions
- References
- 7 Online Material
- Appendix A: Orthogonality of the anelastic eigenfunctions
- Appendix B: Anelastic oscillations of an isothermal atmosphere
- Appendix C: Damping rate of the modes using work integrals
- Appendix D: Parseval-Bessel relation in the anelastic subspace

A&A 421, 775-782 (2004)

DOI: 10.1051/0004-6361:20035606

**B. Dintrans ^{1} - A. Brandenburg^{2}**

1 - Observatoire Midi-Pyrénées, CNRS UMR5572, 14 avenue
Edouard Belin, 31400 Toulouse, France

2 - NORDITA,
Blegdamsvej 17, 2100 Copenhagen Ø, Denmark

Received 31 October 2003 / Accepted 2 April 2004

**Abstract**

The excitation of internal gravity waves by an entropy bubble
oscillating in an isothermal atmosphere is investigated using direct
two-dimensional
numerical simulations. The oscillation field is measured by a
projection of the simulated velocity field onto
the anelastic solutions of the linear eigenvalue problem for the
perturbations. This
facilitates a quantitative study of both the spectrum and the amplitudes
of excited *g*-modes.

**Key words: **hydrodynamics - waves - stars: oscillations - stellar dynamics

The problem of the excitation of internal gravity waves (hereafter IGWs)
in stably stratified media is often studied
in connection with the possible detection of solar *g*-modes
(see e.g. the latest attempt of *g*-mode detection in the GOLF data by
Gabriel et al. 2002).
As an example, two- and three-dimensional numerical simulations
of penetrative convection have shown that it is possible to excite
such waves below the convection zone from the penetration of
strong downward plumes into the stable radiative zone located below
(Hurlburt et al. 1986; Hurlburt et al. 1994; Brandenburg et al. 1996;
Brummell et al. 2002; Dintrans et al. 2003).

Furthermore, IGWs play a major role in stellar evolution as they can transport angular momentum and/or chemical elements over long distances in the stellar interior. This transport mechanism has been proposed to explain the quasi-solid rotation profile of the solar core as revealed by helioseismology (Kumar et al. 1999) as well as the lithium depletion observed in low mass stars (Talon & Charbonnel 1998). However, the correct amount of IGWs excited by, say, penetrative convection remains still unknown as numerical studies were essentially led from a qualitative point of view, so that the mode amplitudes were not deduced.

The main goal of this paper is to show that it is possible to determine
quantitatively the amplitudes of gravity waves propagating in a stable zone
of a numerical simulation using the anelastic subspace. This subspace is
built from the solutions of the associated anelastic linear eigenvalue
problem for the perturbations. The projection of the simulated velocity
field onto this basis yields the IGW amplitudes. We
will present here the application of this method to the simple problem
of *g*-mode oscillations in an isothermal atmosphere. The dynamics
of such an atmosphere is well known (Brandenburg 1988)
so that this problem is ideal to illustrate and test the validity of our method,
before applying it to the more difficult problem of
IGWs excited by penetrative convection, where the atmosphere is in
general non-isothermal.

After presenting the hydrodynamical model describing our isothermal
atmosphere
(Sect. 2), we apply two classical and widely used methods
to detect the *g*-modes excited by an oscillating entropy bubble and
show their limitations (Sect. 3). Next, we introduce our new
method based on the anelastic subspace and give the analytic forms we
found for both the eigenfrequencies and eigenvectors of the anelastic
set of linear equations for the perturbations (Sect. 4).
We then apply this formalism to the same simulation used to test the two
classical methods discussed in Sect. 3 and show that we now have access to
both the spectrum and amplitudes of the IGWs
(Sect. 5). In particular, we illustrate the
usefulness of our method by studying the influence of
the box geometry on the mode amplitudes. Finally, we conclude in
Sect. 6 by introducing the next step of this work, which
will consist of the application of the anelastic subspace method to the
detection of IGWs that are stochastically excited by penetrating convective plumes.

2 The hydrodynamical model

We consider the two-dimensional propagation of internal gravity waves
in an isothermal atmosphere consisting of a layer of depth *d*of a perfect gas at temperature *T*_{0} embedded between two horizontal
impenetrable boundaries. The fluid is non-rotating and stratified with
constant gravity and its properties like its kinematic viscosity, ,
and specific heats, *c*_{p} and *c*_{v}, are assumed to be constant (with
for a monatomic gas). Initially, the layer is in
hydrostatic equilibrium such that its density is given by

where

We choose the depth of the layer *d* as the length unit,
as
the time unit
and
and *T*_{0} as the units of density and temperature, respectively. The evolution
of this layer is governed by the equations of conservation of mass,
momentum and energy

where is the total derivative, the density, the velocity,

When imposing the boundary conditions, we assume the top and bottom
surfaces to be stress-free with the fixed temperature (or,
equivalently, internal energy), that is

In the horizontal direction, we shall impose periodic conditions for all fields.

The fully nonlinear set of Eqs. (2) is solved using the
hydrodynamical code described in Nordlund & Stein (1990) and
Brandenburg et al. (1996). The time
advance is explicit and uses a third order Hyman scheme (Hyman 1979).
All spatial derivatives are calculated using sixth order
compact finite differences (Lele 1992). The two-dimensional simulations presented
here were done with the same resolution
(i.e. 256 zones
in the horizontal direction and 300 in the vertical) and aspect
ratios,
,
between 2 and 6. Here, *L*_{x} and *L*_{z} denote
the domain size in the horizontal and vertical directions, respectively.
Moreover, the
dimensionless sound speed and the gravitational acceleration
are also the same for all
runs and set equal to
.
As a consequence, the dimensionless
pressure scale height is *H*=3/5 (i.e. the box extends over 1.66*H* in the
vertical direction so that the density contrast between bottom and top
is around
), while the Brunt-Väisälä frequency
is
.
Finally, the dimensionless kinematic
viscosity
has been set equal to
in all
simulations, except in the study of the influence of the viscosity
on the mode damping rates (Sect. 5.2) where we
also considered smaller values down to
.

As initial condition to excite IGWs in the layer, we choose the
perturbation to be an entropy bubble (Brandenburg 1988). For a given
entropy perturbation, it is indeed easy to find the approximate position
of the point around which the bubble will begin to oscillate
and thus generate IGWs. In dimensionless units we have

where is the dimensionless entropy. All simulations start with an initial (negative) entropy perturbation located at

Figure 1:
Velocity field superimposed on a grey scale representation
of the entropy perturbation for a two-dimensional simulation of an entropy
bubble oscillating in an isothermal atmosphere.
The asterisks shown in panel a) are used in
Fig. 2a. |

Figure 1 shows an example of such a two-dimensional hydrodynamical simulation
where
the velocity field has been superimposed onto a grey scale representation
of the entropy. At *t*=0, the negative entropy perturbation looks like a bubble at *x*=0and *z*=0.2 (Fig. 1a) and no velocity field is present. Under
the effect of gravity, the bubble begins to descend (Fig. 1b)
and stabilizes at the depth *z*=0.5. The bubble thus oscillates
around this equilibrium position and IGWs are generated in the whole domain
(Fig. 1c). The question is what is the spectrum and amplitude
of the internal waves excited by this oscillating bubble.

3 Measuring IGWs excited by the oscillating bubble from classical methods

Two techniques are commonly used to measure wave fields propagating in direct hydrodynamical simulations:

**Method 1 (hereafter M1):**- the simplest method consists of firstly,
recording the vertical velocity at a fixed point and, secondly, performing
a Fourier transform of the time series. This method was used for example
by Hurlburt et al. (1986) to detect IGWs in a simulation
of penetrative convection.
**Method 2 (hereafter M2):**- a more precise method consists of, firstly,
taking horizontal Fourier transforms of the vertical mass flux
for every time step and, secondly, computing power spectra for
the individual time series. Peaks corresponding to acoustic or gravity
oscillations thus appear at a given horizontal wavenumber
*k*_{x}in the -plane. This method was used for example by Stein & Nordlund (1990) to determine the acoustic modes that are excited in their numerical simulations of solar granulation.

Figure 2:
Detection of IGWs using method M1. a): time evolution of
the vertical velocity recorded at three different positions in the domain
(denoted by stars in Fig. 1a). b): corresponding temporal
power spectra. Peaks below the buoyancy frequency N=0.82 correspond to internal gravity waves whereas the peak around
is the fundamental acoustic mode. |

In Fig. 3 we summarize the application of the second method M2 to detect IGWs in the same simulation of Fig. 2. We now have
access to more informations than with method M1. At first, the diagram confirms: *(i)* that the mode around
is indeed
a radial one and that it corresponds to the fundamental acoustic mode
as no radial node is visible; *(ii)* that gravity modes are strictly
nonradial modes as no peaks are visible below the dotted line .

Second, gravity modes with similar frequency but different degrees are now well separated. As an example, the fundamental *g*-mode at and the first overtone at
have almost the same frequency,
.
Using the
-plane for different 's
allows us to separate modes and to show that one mode
corresponds to a fundamental one while the other one is a first
overtone (note the two "bumps'' around
in the
-diagram for ).

4 A new method using the anelastic subspace

The classical methods M1 and M2 presented above have some disadvantages which render them inappropriate for a careful detection of IGWs in hydrodynamical simulations:

- the main disadvantage of method M1 is its
*k*_{x}-degeneracy, that is, the horizontal dependence of the modes is not taken into account. As a consequence, we have only access to a rough spectrum where only very few peaks appear. - method M2 takes into account this horizontal dependence so
that more modes are detected. But what about modes amplitudes?
For example, this method does not permit to quantify the amount of kinetic energy
due to IGWs propagating in the box so that only qualitative information, such
as the knowledge of the spectrum of excited
*g*-modes, is possible.

Our new method, which eliminates these shortcomings, is inspired by
the work of Bogdan et al. (1993, hereafter BCM) who have
developed a tool to measure the *acoustic* emission generated in
their simulations of turbulent convection. They extracted the
acoustic field by projecting the convection field onto the acoustic
subspace build from the eigenfunctions of the associated linear
oscillations problem. Here we are interested in gravity waves rather
than sound waves and this allows us to adopt a simpler procedure than
theirs. Thus, we project our simulation data onto the anelastic
subspace, that is, we filter out the acoustic waves in the linear system
of oscillation equations. For a time dependence of normal modes of the
form
,
and in the ideal limit ,
the anelastic
set of equations reduces to (Dintrans & Rieutord 2001; Rieutord &
Dintrans 2002)

where denotes the Lagrangian displacement (the associated velocity field being ), the Eulerian perturbation in reduced pressure and the equilibrium density profile.

Because we adopt periodic boundary conditions in the *x*-direction,
we can seek solutions of the form

such that the anelastic system (3) reads after the elimination of

This last system can be formally rewritten as a generalized eigenvalue problem of the form

where is the eigenvector associated with eigenvalue and , denote two differential operators.

As shown in Appendix A, the eigenfunctions
which are
solutions of the eigenvalue problem (5) are orthogonal,
that is we have

for normalized eigenfunctions (here the symbol denotes the Hermitian conjugate). As a consequence, the set forms an orthogonal basis onto which the simulated velocity field can be projected, i.e.

where the -sum begins at as no gravity modes propagate radially, i.e. should be non-zero (see Sect. 3). In this decomposition, the first term corresponds to the IGW contribution while all other rates, such as the bubble displacement or the acoustic waves, were collected in a term which we identify hereafter as "rest". Finally, following BCM, we identify the amplitude of each gravity mode by the time-dependent coefficient

In the case of an isothermal atmosphere, for which and

where we have introduced the vertical wavenumber . Because of these analytic expressions, the construction of the anelastic basis used in the projection (7) is straightforward. However, in practice, we do not project directly onto the two-dimensional velocity field but rather its horizontal Fourier transform at a given

- 1.
- given a
*k*_{x}-value, we first deduce the anelastic eigenvalues and eigenvectors from their analytic expressions (9); - 2.
- we build the anelastic subspace at this horizontal wavenumber
*k*_{x}with, say, the first five anelastic modes*n*=[0,4]; - 3.
- we Fourier transform the simulated velocity field
to obtain the new complex-type field
;
- 4.
- we project
onto the anelastic basis, built
at step 2, and deduce the mode amplitudes (8)
for each order
*n*.

5 Results

Figure 4 shows four anelastic modes at
(a) and
(b) with *n*=0 (dark lines) and *n*=1 (grey lines). We remark that
anelastic eigenfunctions hardly change with the -values as does not explicitly depend on *k*_{x}, while only the -amplitude
involves *k*_{x}. As a consequence, it is the ratio
which makes sense so that
becomes more and more
important compared to
with increasing *k*_{x}, as observed in
Fig. 4 when comparing solid to dot-dashed lines for the
same order *n*. In the case of method M2,
the use of the vertical mass flux
(instead of, say,
)
is justified a posteriori
as a good proxy for detecting IGWs propagating in the
simulation (see Fig. 3).

Figure 5:
Upper: same as Fig. 3 at ,
except that amplitudes
have been restored using third direction. Lower: comparisons between
the normalized vertical mass flux
profiles deduced from
the upper power spectra (solid lines) with the ones obtained from the
anelastic modes shown in Fig. 4a (dotted lines). |

Once the anelastic eigenfunctions are known, one can compare the corresponding
vertical mass flux
with the ones deduced from the power
spectra used in method M2. Indeed, the peaks appearing at a given
frequency in Fig. 3 correspond to *g*-modes so that their vertical
profiles can be related to the eigenfunctions of the linear modes. This
was done by BCM for the case of acoustic modes propagating
in their convection simulations.

The upper panel of Fig. 5 shows the depth dependence of the
two peaks appearing at
in Fig. 3, which correspond to
the *n*=0 (
)
and *n*=1 (
)
modes, respectively. We thus recover what BCM called the "shark
fin'' profiles, that is, the power at any depth is located around
the theoretical anelastic frequency. Taking a mean profile around
eigenfrequencies
and
allows us to compare
the vertical mass flux observed in the simulation with those calculated
from the linear anelastic modes. The lower panel of Fig. 5
shows that the agreement between the two profiles is quite remarkable,
meaning that the anelastic modes reproduce well the long-period dynamics of the
oscillating isothermal atmosphere.

5.2 Time evolution of the mode amplitudes

Figure 6:
Left: real parts of the complex mode amplitude
for the three modes shown in Fig. 4. Right:
corresponding temporal power spectra. Dotted lines mark the location
of the anelastic eigenmodes
and
,
respectively. Amplitudes on each plot have
been multiplied by one thousand. |

The left hand panel of Fig. 6 shows the real parts of the complex
amplitudes
and *c*_{21}, that is,
the result of the projection of the simulation shown in Fig. 1
onto the three anelastic modes plotted in Fig. 4. By taking
a Fourier transform of these sequences (Fig. 6, right), one obtains
three peaks which agree remarkably well with the theoretical anelastic
frequency given by (9), i.e. the complex amplitudes
behave as

This relation means that when an eigenmode is excited in the simulation, its corresponding displacement is periodic with a period .

Figure 7:
Time-evolution of the amplitude modulus
for the
three modes shown in Fig. 4. Grey lines superimposed
correspond to an exponential decay of the form
. |

As we have access to the real and imaginary parts of , we now compute the amplitude modulus ; Fig. 7 shows its time evolution for the three modes in Fig. 6. A regression analysis of the curves allows us to determine the shape of the envelope and we found that follows a simple exponential decay law of the form . Because of the periodic behaviour (10), it means that the mode amplitude behaves like

where is a constant which is related to the damping process, here the viscosity. Using runs with different kinematic viscosities we have verified that , i.e. every mode obeys the same law as that of a linear damped harmonic oscillator. Such a relation between the damping of the mode and the viscosity is hardly surprising as it can easily be deduced from work integrals; see Appendix C.

We will now show that one of the main advantages of the anelastic subspace
resides in the fact that it allows to calculate the contribution of
every *g*-mode to
the time-averaged kinetic energy.
We demonstrate in Appendix D that a Parseval-Bessel relation
also exists in the anelastic subspace as

where, as before, the rest term contains all the contributions which are not due to IGWs.

This relation illustrates the usefulness of the anelastic subspace when
one wants to perform a detailed study of the IGW excitation. By using
the classical Parseval-Bessel relation applied in Fourier space,
it is already possible to find the amount of kinetic energy contained at
a given wavenumber *k*_{x}, since

For example, the dependence of the time-averaged kinetic energy on the degree , shown in Fig. 8, demonstrates that half of the kinetic energy in the simulation in Fig. 1 is contained in the degree (the solid line denoted by "

The main advantage of our anelastic subspace method is that it lifts this
degeneracy by isolating contributions that come from *g*-modes *only*.
Comparing
to its classical form (13), the anelastic Parseval-Bessel relation
(12) introduces the radial order *n* so that one can extract
the contribution of the *g*-mode
in the kinetic energy and then
deduce which modes contribute most to it, meaning the modes that are
the most excited by the oscillating bubble.
This is what we did in Fig. 8 where we also plotted, for each ,
the
contributions coming from
the *n*=0, *n*=1 and *n*=2 *g*-modes. It allows us to show that the contribution is in fact almost entirely composed of the (1,0) and (1,1) *g*-modes, as they respectively contain around 33% and 15%
of the time-averaged kinetic energy of this simulation. This kind of information
is particularly relevant when one deals with the problem of wave
excitation.

In order to assess the restriction imposed by the finite size of the box,
we now study the influence of the box geometry onto the wave excitation.
We therefore perform three runs with three different aspect ratios (2, 4 and 6)
and calculate the *g*-mode contributions to the time-averaged kinetic energy
in each case.

Figure 9 shows the spectral power in the
plane
for the three different aspect ratios.
One recovers clearly the well-known dispersion relation for *g*-modes
(e.g. Stein & Leibacher 1974); see also its anelastic counterpart
in Eq. (9).
One also verifies that
increasing the horizontal extent of the box allows the power to
be distributed among a larger number of horizontal modes. This phenomenon
has also been observed in the problem of *p*-modes excitation by turbulent
convection (Stein & Nordlund 2001; Nordlund & Stein 2001).
Unfortunately, such diagrams are generally not easy to interpret when the
system is not isothermal and the modes stochastically driven; see for
example the three-dimensional simulations of Brandenburg et al. (1996)
where the resulting
turned out to be very noisy.

Figure 9:
diagrams which three different aspect ratios (2, 4 and 6).
Crosses correspond to the anelastic frequencies given by (9). |

Contrasting this with the anelastic subspace method, it is straightforward
to extract at a given wavenumber *k*_{x} the energy from modes with different
radial order *n*.
As expected, the spectral distribution of mode energies integrated over
all *n* collapses onto an universal curve (Fig. 10).
It turns out that the mode energy scales as
for small
enough values of *k*_{x}, i.e. for
.
However, both
the exponent and the cut-off are not unique and depend on the size of
the initial entropy bubble, as verified by running several simulations
with bubble radii ranging from *R*=0.1 to *R*=0.4. We found that the
larger the bubble, the smaller is the cutoff value of *k*_{x}, meaning that big
bubbles preferentially generate *g*-modes at large scales. It would be
interesting to investigate in more detail the origin of such behaviour,
but this is clearly beyond the scope of this paper.

Figure 10:
Distribution of the energy contained in g-modes for
different aspect ratios (2, 4 and 6) in a log-log representation. |

6 Conclusions

We have presented a new and accurate method to measure the internal gravity waves propagating in a numerical simulation. This is of prime importance when one deals with the problem of the generation of IGWs in a stably stratified medium such as the radiative zone of a star. The knowledge of the vigour of the velocity field in such a zone is the main input of any theory of wave transport of angular momentum and/or chemicals in stellar radiative zones (e.g. Schatzman 1993).

However, we have shown that classical methods based only on Fourier
transforms in space and/or time of the simulated velocity field do not
allow a quantitative determination of the mode amplitudes. The reason is
mainly that the contribution of the propagating *g*-modes to the total
energy in the simulation is diluted among the other ones due to *p*-modes
and turbulent motions.

Our method is based on the projection of the simulated velocity field onto
the anelastic subspace built from the solutions of the linear problem
for the perturbations. Once the resulting time-dependent coefficients of
the velocity field are computed, one can deduce the amplitudes
of every *g*-mode and measure their kinetic energy
at every timestep of the simulation.

We have applied this formalism to the measurement of IGWs excited by an
entropy bubble oscillating in an isothermal atmosphere. We have shown that
the mode amplitudes follow the same law than those of a linear
damped harmonic oscillator, that is, the periodic amplitudes decrease
exponentially with time, with a time scale that is simply proportional to
the inverse of the viscosity; see Eq. (11). By using a
Parseval-Bessel
relation valid in the anelastic subspace, we have deduced the contribution
of every *g*-mode to the time-averaged kinetic energy and have shown
that large-scale motions are mainly composed of *n*=0 and *n*=1
modes.
Such a result is in general impossible to obtain in a more complicated
settings using a classical
diagram.

The fact that the mode amplitudes in the anelastic basis follow the same law than those of a linear damped oscillator is important in the case of a stochastic excitation of IGWs, such as the one encountered when one deals with numerical simulations of penetrative convection. Indeed, in this case the mode amplitude evolves either chaotically around zero when the mode is not excited, or in a periodic fashion following Eq. (11) when the mode is excited. By computing time-frequency diagrams of that time-dependent amplitude, one can extract the time intervals during which the mode has been really excited. In a forthcoming paper, we will present the application of the anelastic subspace method we tested here to the problem of IGWs stochastically excited by penetrative convection (Dintrans et al. 2004).

This work has been supported by the European Commission under Marie-Curie grant No. HPMF-CT-1999-00411. Calculations were carried out on the CalMip machine of the "Centre Interuniversitaire de Toulouse'' (CICT) which is gratefully acknowledged. We also thank Michel Rieutord and the referee for useful comments.

- Bogdan, T. J., Cattaneo, F., & Malagoli, A. 1993, ApJ, 407, 316 NASA ADS
- Brandenburg, A. 1988, A&A, 203, 154 NASA ADS
- Brandenburg, A., Jennings, R. L., Nordlund, Å., et al. 1996, J. Fluid Mech., 306, 325 NASA ADS
- Brummell, N. H., Clune, T. L., & Toomre, J. 2002, ApJ, 570, 825 NASA ADS
- Cowling, T. G. 1941, MNRAS, 101, 367 NASA ADS
- Dintrans, B., Brandenburg, A., Nordlund, Å., & Stein, R. F. 2003, Ap&SS, 284, 237 NASA ADS
- Dintrans, B., Brandenburg, A., Nordlund, Å., & Stein, R. F. 2004, in preparation
- Dintrans, B., & Rieutord, M. 2001, A&A, 324, 635 NASA ADS
- Gabriel, A. H., Baudin, F., Boumier, P., et al. 2002, A&A, 390, 1119 NASA ADS
- Hurlburt, N. E., Toomre, J., & Massaguer, J. M. 1986, ApJ, 311, 563 NASA ADS
- Hurlburt, N. E., Toomre, J., Massaguer, J. M., & Zahn, J. P. 1994, ApJ, 241, 245 NASA ADS
- Hyman, J. M. 1979, in Adv. in Computation Methods for Partial Differential Equations, Vol. III, ed. R. Vichnevetsky, & R. S. Stepleman, Publ. IMACS, 313
- Kumar, P., Talon, S., & Zahn, J.-P. 1999, ApJ, 520, 859 NASA ADS
- Lamb, H. 1909, Proc. Lond. Math. Soc., 7, 122
- Lamb, H. 1932, Hydrodynamics (Dover Publications)
- Lele, S. K. 1992, J. Comp. Phys., 103, 16 NASA ADS
- Nordlund, Å., & Stein, R. F. 1990, Comput. Phys. Commun., 59, 119 NASA ADS
- Nordlund, Å., & Stein, R. F. 2001, ApJ, 546, 576 NASA ADS
- Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1994, Numerical Recipes in Fortran. The Art of Scientific Computing. 2nd Edition (Cambridge University Press)
- Rieutord, M., & Dintrans, B. 2002, MNRAS, 337, 1087 NASA ADS
- Stein, R. F., & Leibacher, J. 1974, A&ARv, 12, 407 NASA ADS
- Stein, R. F., & Nordlund, Å. 1990, Progress in Seismology of the Sun and Stars. Lecture Notes in Physics (Berlin: Springer-Verlag), 367, 93
- Stein, R. F., & Nordlund, Å. 2001, ApJ, 546, 585 NASA ADS
- Schatzman, E. 1993, A&A, 279, 431 NASA ADS
- Talon, S., & Charbonnel, C. 1998, A&A, 335, 959 NASA ADS
- Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial oscillations of stars (University of Tokyo Press)

7 Online Material

Appendix A: Orthogonality of the anelastic eigenfunctions

We start from the complete set of equations written for two eigenmodes
and
(Dintrans & Rieutord 2001)

(A.1) |

where denotes the Eulerian fluctuation of the reduced pressure. Multiplying the -equation by and the -equation by , and subtracting the resulting two equations leads to

(A.2) |

where we have used the relation , as the anelastic approximation implies .

Finally, we integrate this last equation over the volume and use the Green-Ostrogradsky theorem to obtain

(A.3) |

The applying of periodic horizontal boundary conditions and the condition at

(A.4) |

Appendix B: Anelastic oscillations of an isothermal atmosphere

We start from the anelastic set of Eqs. (4)

where, for an isothermal atmosphere, and

All coefficients of this differential equation being constants, we can seek solutions of exponential type to obtain the characteristic equation

(B.3) |

whose general solution is of the form

leading to the following solution for

(B.5) |

To find the quantization rule for , we consider the boundary conditions on ,

(B.6) |

where is the radial order of the mode. (The presence of the

and their associated eigenfunctions

(B.8) |

The remarkable fact to be noted here is that the vertical lagrangian displacement does not depend on

B.2 Comparisons with the complete case

Order n |
rel. error | ||

0 | 0.572 054 | 0.567 455 | 0.804% |

1 | 0.363 084 | 0.362 606 | 0.132% |

2 | 0.257 381 | 0.257 295 | 0.033% |

3 | 0.197 644 | 0.197 621 | 0.012% |

4 | 0.159 920 | 0.159 912 | 0.005% |

The set of equations for the complete case is obtained by, first, linearizing the system (2) for infinitesimal adiabatic perturbations and, second, by assuming the time and horizontal dependences of the modes to be and or , respectively. We thus arrive at

Here and are the Eulerian perturbations in reduced pressure and density, respectively. This system has been obtained under the classical Cowling approximation which assumes that perturbations in the gravitational potential are negligible (Cowling 1941). We refer the reader to the book of Unno et al. (1989) for a complete derivation of this system in the case of a spherical geometry.

Lamb (1909, 1932) found analytical solutions of the system (B.9) in the case of an isothermal atmosphere. By eliminating and in the system (B.9), we obtain a second-order differential equation for only,

Comparison with Eq. (B.2) shows that the anelastic approximation applies provided that:

- the Lamb frequency term
is negligible, implying that

which is a condition naturally fulfilled by high-order*g*-modes with ; - the ratio
is very large.
^{}

while the same quantization rule applies due to the identical boundary conditions , that is,

We thus obtain the following quadratic dispersion relation for the complete eigenfrequencies

which is the same than that in Stein & Leibacher (1974). In the same way, the eigenvectors of the complete case can be deduced from the anelastic ones and we find

Despite its different -term (B.12), we note that in the complete case the vertical displacement is the same as in the anelastic one, which is simply a consequence of applying the same quantization rule (B.13).

Figure B.1:
Comparison of the complete (solid lines) and anelastic
(dashed lines) normalized eigenfunctions
for
and n=[0,1,2]. |

We summarize in Table B.1 the eigenfrequencies for the complete case deduced
from (B.14) for
the first five *g*-modes of degree
and compare them with their anelastic counterparts given by the
formulæ (B.7). The agreement is quite remarkable since the
relative error done by the anelastic approximation is less than from the fundamental *g*-mode. Concerning the associated eigenvectors,
we showed that the vertical displacements
are the same in both
cases so that only the horizontal ones differ.
However, Fig. B.1 shows that the agreement on the normalized
-eigenfunctions is also very good, except for the mode at *n*=0for which the condition (B.11) is not well fulfilled.

In conclusion, the filtering of acoustic waves from the dynamics of the
isothermal atmosphere does almost not change the *g*-mode eigenfrequencies
and eigenfunctions
.
The basis
built from these anelastic eigenfunctions can therefore be used with
good confidence to project out velocity fields from hydrodynamical simulations.

It should be pointed out that the application of the anelastic subspace method will only give good results when the acoustic and gravity spectra are well separated, that is, when the ratio is large. This is the case for the chosen isothermal setup, where the resulting anelastic eigenfrequencies/eigenvectors are very close to those of the complete problem (see Table B.1 and Fig. B.1). However, this ratio decreases with larger stratification, in which case one may be forced to solve numerically the oscillation Eqs. (B.9) of the complete problem instead of those of the anelastic subset (B.1), which does not however introduce any particular difficulty.

Appendix C: Damping rate of the modes using work integrals

We show in this appendix that it is possible to find the relation
governing the damping rate of a mode and the viscosity using the work
integrals formalism. We start from the anelastic equations written for
the velocity

(C.1) |

where has a non-zero real part due to the viscous dissipation. We note that the viscous term is normally more complicated as there may be a vertical variation of the dynamical viscosity (see the general form of this term in Eq. (2)). However, we take here a simpler form as our aim is just to formally prove that the damping rate is linearly related to the viscosity .

We multiply the momentum equation by to obtain

(C.2) |

where we have used . As in Appendix A, we finally integrate this last equation over the volume and use the Green-Ostrogradsky theorem to obtain

(C.3) |

where we also used . As we only have real integrals on both sides, we can deduce the following relation for the real part

(C.4) |

The damping rate of a mode is therefore linearly related to the viscosity, as found numerically in (11).

A similar relation also applies in the complete case, provided one includes the contribution of pressure fluctuations to the energy, that is,

where the total energy contains the kinetic ( ), potential ( ) and internal ( ) contributions. Such relations are useful from a numerical point of view as they allow to check the accuracy of the computed eigenvalues from their associated eigenvectors.

Appendix D: Parseval-Bessel relation in the anelastic subspace

In this appendix, we will show that it is possible to relate the total
kinetic energy embedded in the box to the one calculated in the anelastic
subspace. To do that, we need a similar relation than the Parseval-Bessel
one which states that the total energy is conserved when one works in
Fourier space, i.e.

(D.1) |

where denotes the Fourier transform of

In our case, we only perform a 1-D Fourier transform of the velocity field in the horizontal direction (corresponding to a wavenumber ), so that the Parseval-Bessel relation for the total kinetic energy at a given time is

where denotes the velocity field at a wavenumber

(D.3) |

The

= | |||

= | |||

(D.4) | |||

= | |||

= |

This last equation means that the total kinetic energy due to

(D.5) |

Copyright ESO 2004