The Eﬀects of Large-Scale Convection on Solar
Eigenfrequencies
byMichael Marchand Swisdak, III
B.S., Astronomy, University of Maryland, 1994
B.S., Mathematics, University of Maryland, 1994
B.S., Physics, University of Maryland, 1994
M.S., Astrophysics, University of Colorado, 1996

A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulﬁllmentof the requirements for the degree ofDoctor of Philosophy
Department of Astrophysical and Planetary Sciences
1999

This thesis entitled:
The Eﬀects of Large-Scale Convection on Solar Eigenfrequencies
written by Michael Marchand Swisdak, III
has been approved for the Department of Astrophysical and Planetary Sciences

Ellen G. Zweibel

Juri Toomre

Date

The ﬁnal copy of this thesis has been examined by the signatories, and we ﬁnd that both thecontent and the form meet acceptable presentation standards of scholarly work in the above
mentioned discipline.

iii
Swisdak, III, Michael Marchand (Ph.D., Astrophysics)
The Eﬀects of Large-Scale Convection on Solar Eigenfrequencies
Thesis directed by Professor Ellen G. Zweibel

We describe and implement an approach for determining the eigenfrequencies of solar
acoustic oscillations (p modes) in a convective envelope. By using the ray approximation, we
transform the problem into one in which we seek the eigenfrequencies of a Hamiltonian system. To ﬁnd these eigenfrequencies we have written a computer program which implements themethod of adiabatic switching. In this technique, we begin with a system with no convective
perturbations for which the eigenmodes and eigenfrequencies are known. The code adiabatically increases the strength of the convective structures, allowing the mode eigenfrequency toadjust from its initial value to the eigenfrequency of the perturbed state. The ray approximation
restricts our investigations to perturbations which are large compared to the mode wavelength.
For a simple class of structures we test our results against the predictions of semi-classical
EBK quantization and ﬁnd the two methods agree. We then examine more complicated perturbations, concentrating on the dependence of the frequency shifts on the radial and angularmode numbers as well as the perturbation strength. Among our results, we conclude that thefractional frequency shift is given by the weighted average of the perturbation over the resonantcavity. As a result, convective perturbations with horizontally anti-symmetric structures generate downward frequency shifts which are second-order in the perturbation strength. We also
examine more complex convective structures which we ﬁnd tend to produce downshifts whose
magnitude scales with the strength of the perturbation. These results may have implications forresolving the diﬀerences between eigenfrequencies derived from solar models and those deduced
from helioseismic observations.

Dedication

To MishaAcknowledgements

This work was supported by NASA GSRP 153-0972, NSF Grant AST-95-21779, theNASA Space Physics Theory Program, and the SOHO/MDI Investigation. Mark Rast and Nic
Brummell kindly provided me with the plume model and convective simulation used in §4.3 and§4.4, respectively. Deborah Haber, Jim Meiss, and Rex Skodje provided helpful discussions at
critical points along the way.
In addition to the professional acknowledgments, I would like to use this space to thank a
few of the people whose inﬂuences have been instrumental in seeing this dissertation to fruition.
Naturally this is a daunting task; fortunately the associated memories make it an enjoyable one.
Ellen Zweibel, my advisor, well deserves ﬁrst mention on this list for giving a nervousﬁrst-year graduate student a job and then nurturing me through the intervening years. Hopefully
I have managed to fulﬁll at least some small part of her expectations as a student, for she has
met and exceeded all of mine as an advisor.
My time at Colorado resonates with great memories, a happy circumstance in part attributable to the students with whom I have had the pleasure of working. Even if speciﬁcmemories — conversations, scientiﬁc and otherwise, social events, and so many Ultimate games— begin to fade, I will always associate a warm glow with my time here. In that spirit, and witha nod to my fear of inadvertently neglecting anyone, I extend a comprehensive acknowledgment
to my fellow students.
A few merit special notice. No thesis is written without journeys along roads which,while initially promising, lead to treacherous dead-ends. I am lucky enough to have found several

vi
people who were tolerant enough to listen to my oﬀ-track meanderings and, more often than not,gently nudge me onto the correct path. Marc DeRosa has gracefully consented to my, at times
seemingly never-ending, barrages of questions, theories, and explanations. Kelsey Johnson, theclosest thing to a sister I have ever had, provided uncountable words of encouragement. Rémy
Indebetouw showed me the power of forthrightness. And, my oﬃcemates Eli Michael and Rekha
Jain, who coped with frantic mumblings and assorted eccentricities, deserve special praise for
their patience in the midst of it all.
I have been further blessed with a wonderful family. My parents, Michael and Carole,
have been the single biggest inﬂuence on my life; despite not appearing on the title page, they
are as much the authors of this document as I. Stephen, my brother, must be commended forhis grace during the onerous task of growing up with me as a sibling. In light of this burden, his
accomplishments are even more impressive. And ﬁnally I am grateful to Misha, my wife, whose
love sustained me, not only through the low points, but along every step of the way.
Thank you all.

Contents

Chapter
1 Introduction1

1.1

Problem Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.2

Solution Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.3

Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92 Theoretical Justiﬁcation
2.1

12

Stellar Oscillations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.1.1

Fluid Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.1.2

Perturbations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

2.1.3

Separation of Angular Variables . . . . . . . . . . . . . . . . . . . . . . .

15

2.1.4

Derivation of a Wave Equation . . . . . . . . . . . . . . . . . . . . . . . .

16

2.2

Asymptotic Approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

2.3

Hamiltonian Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

2.4

Quantization Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

2.4.1

EBK Quantization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

2.4.2

A Simple Physical Model . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

2.4.3

P Mode Quantization Conditions . . . . . . . . . . . . . . . . . . . . . . .

30

2.4.4

Surfaces of Section . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

Adiabatic Switching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

2.5

viii
3 Implementation and Tests
3.1

39

The Reference State . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

39

3.1.1

Reference State Results . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

3.1.2

Quantization Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

3.2

Adiabatic Switching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

3.3

A Test of an Integrable System . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

3.4

Doppler Shifts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

4 Convective Structures

57

4.1

Cell-like Velocity Perturbation . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

4.2

Rayleigh-Bénard Boussinesq Convection . . . . . . . . . . . . . . . . . . . . . . .

64

4.3

Compressible Plume . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

4.4

Turbulent Compressible Convection . . . . . . . . . . . . . . . . . . . . . . . . .

81

5 Conclusions

88

5.1

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

88

5.2

Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

91

Bibliography

94

Appendix
A Table of Important Symbols

97

B Fluid Descriptions and Perturbations

100

B.1 Descriptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
B.2 Mathematical Operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
B.3 Perturbations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
B.4 Commutation Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

ix
B.4.1 Eulerian Perturbation, Commuting Operators . . . . . . . . . . . . . . . . 107
B.4.2 Eulerian Perturbation, Non-Commuting Operators . . . . . . . . . . . . . 107
B.4.3 Lagrangian Perturbation, Commuting Operators . . . . . . . . . . . . . . 108
B.4.4 Lagrangian Perturbation, Non-Commuting Operators . . . . . . . . . . . 108
C Alternate Forms of the Wave Equation

110

D The Computer Code114

Figures

Figure1.1

The diﬀerence between observationally and theoretically determined values of therelative squared sound speed versus solar radius. . . . . . . . . . . . . . . . . . .

4

2.1

Acoustic raypaths in a two-dimensional Sun.

29

3.1

A sample raypath in an adiabatically stratiﬁed, two-dimensional, polytropic en-. . . . . . . . . . . . . . . . . . . .

velope with µ = 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

3.2

A raypath and the evolution of the eigenfrequency for the perturbation of §3.3. .

52

3.3

A raypath and the evolution of the eigenfrequency for the perturbation of §3.4. .

56

4.1

Streamlines for the velocity perturbation of §4.1. . . . . . . . . . . . . . . . . . .

59

4.2

A raypath and the evolution of the eigenfrequency for the perturbation of §4.1. .

65

4.3

Streamlines and isotherms for the convective perturbation of §4.2. . . . . . . . .

67

4.4

A raypath and the evolution of the eigenfrequency for the perturbation of §4.2. .

73

4.5

Sound speed perturbations for the plume of §4.3. . . . . . . . . . . . . . . . . . .

75

4.6

A raypath and the evolution of the eigenfrequency for the perturbation of §4.3. .

78

4.7

Fractional frequency shift versus strength of the plume of §4.3. . . . . . . . . . .

79

4.8

Horizontal slices of the plume of §4.3. . . . . . . . . . . . . . . . . . . . . . . . .

80

4.9

Sound speed and velocity perturbations of the convective model of §4.4. . . . . .

82

4.10 A raypath and the evolution of the eigenfrequency for the perturbation of §4.4. .

84

4.11 Horizontal slices of the convective model of §4.4. . . . . . . . . . . . . . . . . . .

85

xi
4.12 Fractional frequency shift versus strength of the convective model of §4.4. . . . .

87

B.1 A pictorial demonstration of the Eulerian and Lagrangian perturbations. . . . . . 104

xii

TablesTable
3.1

Frequency shifts for the perturbation of §3.3. . . . . . . . . . . . . . . . . . . . .

50

3.2

Frequency shifts for the perturbation of §3.4. . . . . . . . . . . . . . . . . . . . .

54

4.1

Frequency shifts for the convective perturbation of §4.1. . . . . . . . . . . . . . .

60

4.2

Frequency shifts for the convection cells of §4.2 with dominant horizontal velocity
perturbations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.3

70

Frequency shifts for the convection cells of §4.2 with equal strength horizontal
velocity and sound speed perturbations. . . . . . . . . . . . . . . . . . . . . . . .

72

A.1 Important symbols. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97

D.1 Modules used in our adiabatic switching code. . . . . . . . . . . . . . . . . . . . . 115

Chapter 1

Introduction

Plasma in a thin exterior layer, the photosphere, produces all of the observable light fromthe Sun; photons generated interior to that layer are unable to escape the Sun without being
absorbed and re-emitted. This limitation restricted early observations, but astrophysicists could,
and did, construct models of the solar interior based on the available data (primarily the Sun’smass, radius, luminosity, photospheric composition, and age). Unfortunately, their predictionsof the interior structure could not be directly conﬁrmed. Only in the last 25 years has it beenpossible to peer beneath the surface of the Sun and examine its interior structure. The newsource of information is the discipline of helioseismology.
Leighton et al. (1962) used observations of the Doppler shifts of photospheric spectral lines
to detect the ﬁrst evidence of solar oscillations. (Modern observations use either a modiﬁcationof their technique or another method which searches for variations in the total intensity ofemitted radiation.) Although other mechanisms were considered, it was eventually proposed(Ulrich 1970; Leibacher & Stein 1971) that acoustic waves forming resonant oscillations, ormodes, within the Sun were the source of the signal. Since pressure provides the restoring force
for these oscillations, they are now known as p modes. The Sun, like an organ pipe, is a resonant
structure capable of sounding only particular notes. In general, an object’s structure dictates
which notes, or eigenfrequencies as they are called on the Sun, are generated. For an organ
pipe the key structural element is its length: A pipe with a C fundamental can also soundthe C an octave higher, but will never produce the intervening E  . Similarly, the composition

2and structure (particularly the speed of sound) of the Sun determine its eigenfrequencies. Oncecorrectly identiﬁed, it was quickly realized that p modes provide an excellent mechanism fordirectly exploring the solar interior since p-mode frequencies depend on the solar structure muchas organ notes depend on the length of the pipe. Thus, measurements of photospheric quantities
permit the probing of the solar interior. Since terrestrial earthquakes produce acoustic waves
which are subjected to a similar analysis, the new discipline was termed helioseismology.
Although recent helioseismic investigations conclusively demonstrate that modern solarmodels closely describe the Sun, there are still statistically signiﬁcant discrepancies between thepredictions and the observations. In this work, we investigate a probable cause of some of thediscrepancy: large-scale convection as it aﬀects p-mode frequencies. This ﬁrst chapter is devoted
to the presentation of introductory material. A more detailed outline of the problem is given in§1.1, while §1.2 describes our attempt to resolve it. In the context of helioseismology, we believe
our methodology is a novel approach to the problem. However, several antecedent works haveattacked the issue from other directions; a review of the relevant literature may be found in §1.3.
Subsequent chapters develop the solution method and its applications. Chapter 2, in
which we discuss the mathematical underpinnings of our work, is a visually imposing but crucialcomponent, as it delineates the assumptions on which the following chapters are based. In
Chapter 3 we use two unrealistic, but simple, models to illustrate both the implementation of
the method and some general results. The more detailed convective models of Chapter 4 aresimultaneously more interesting and more diﬃcult to interpret. In particular we consider a cellcomposed solely of velocity perturbations, a Rayleigh-Bénard model of convective cells, a cold
plume descending in a compressible atmosphere, and a slice of a model of turbulent compressible
convection. Chapter 5 contains a summary of the results as well as our thoughts concerning
future applications of the method.

3

1.1

Problem Outline
In broad terms, astrophysicists understand the structure of the Sun. The inner 20%, byradius, is known as the core and produces essentially the entirety of the Sun’s energy throughthermonuclear fusion. Photons carrying this energy slowly diﬀuse through the radiative interior
to approximately 0.7R ,1 where a sharp transition occurs. At this radius the Sun becomes
convective, meaning that energy transport occurs primarily through bulk ﬂuid motions rather
than photon diﬀusion. The convection zone extends to the solar surface, the photosphere,
whereupon photons free-stream outward carrying the Sun’s energy through the solar atmosphere— the chromosphere and corona — and into space.
Although this basic description of solar structure is almost universally accepted, certain
details remain contentious. After creating a solar model astrophysicists can compare its structurewith the solar structure as inferred from p-mode eigenfrequencies. As can be seen in Figure 1.1,
the diﬀerences between the two are small but signiﬁcant. The largest discrepancies occur in three
distinct regions: near the solar core, at the tachocline (the boundary between the radiativeinterior and the convective envelope), and near the solar surface. This work focuses on thediscrepancy near the surface which is thought to be partially a consequence of the incomplete
treatment of convective motions, both in computing the solar models and in analyzing thehelioseismic data. On one hand, solar models based on mixing-length theories of convectionslightly miscalculate the interior solar structure and thus disagree with the true Sun. On theother hand, inversions of p-mode eigenfrequencies to determine the solar structure ignore theeﬀects of convective perturbations. Although both eﬀects contribute to the discrepancies seenin Figure 1.1, we will primarily concern ourselves with the latter in this work.
Solar convection is notoriously diﬃcult to model. Important physical processes occur
on lengthscales varying by at least 6 orders of magnitude and current simulations (Brummell
et al. 1995) fall a factor of more than a billion short of resolving these diﬀerent scales. Solarmodelers, long recognizing the inadequacy of numerical simulations, developed and embraced a
1 The symbol R stands for radius and the subscript  symbolizes the Sun.

4

Figure 1.1: The diﬀerence between observationally and theoretically determined values of therelative squared sound speed versus solar radius. The observational values come from inversionsof p-mode eigenfrequencies while the theoretical values come from a solar model. Courtesy of
SOHO/MDI consortium. SOHO is a project of international cooperation between ESA and
NASA.

5
parameterization of the convection zone known as mixing-length theory (see Spiegel 1971, for areview of the ﬁeld’s development), in which all of the complex convective motions are collapsed
into one parameter, the mixing length. Mixing-length theory enjoys the dual beneﬁts of being
conceptually simple as well as generally correct — solar models which employ it do a remarkably
good job of matching observations. In reality however, the outer portion of the convection
zone harbors turbulent convection on scales at least as small as granules (1 Mm) and as big as
supergranules (30–50 Mm). In addition, giant cells spanning the entire convection zone (200 Mm)
may exist (Beck et al. 1998). The existence of diﬀerent scales of motion implies mixing-length
theory is, at some level, inadequate. Since it is reasonable to assume that the solar convection
zone is not precisely described by mixing-length theory, one naturally anticipates disagreements
between observational data and theoretical models of the solar sound speed structure.
The inversion methods which deduce the true solar structure of Figure 1.1 from p-mode
eigenfrequencies contribute a separate set of errors. Although in principle it is possible to use
only the information contained in the eigenfrequencies to reconstruct the solar interior,2 theprocess is computationally diﬃcult. In practice, inversions begin with a reference model suchas one of the mixing-length models described above. Eigenfrequencies are computed for thereference model and compared to the solar data; discrepancies are resolved by perturbing thestructure of the model so that its eigenfrequencies agree more closely with the observations.
Iteration produces the ﬁnal result, the inverted solar structure (which is hopefully independentof the initial reference model). Any errors in the model stratiﬁcation will be corrected in theiteration and will not aﬀect the inverted structure. However, current inversion techniques ignore
the eﬀects of convective perturbations. Since, as we show in this work, convective structures do
aﬀect p-mode eigenfrequencies, not treating their eﬀects will produce errors in the inversion for
the true solar structure.
Thus, two possible errors can produce the discrepancies of Figure 1.1. Models of the solar
convection zone created with mixing-length theories do not completely describe the average
2 Sound speed and density are the most common variables sought in inversions, but other choices, such as thehelium abundance, are possible.

6
interior structure. However, even if models employing mixing-length theory were completely
accurate, the true solar structure to which they should be compared is unknown since helioseismic
inversions do not account for the eﬀects of convective structures on p-mode eigenfrequencies. In
this work we concentrate on the second of these eﬀects and demonstrate that convective motions
change p-mode eigenfrequencies by a helioseismically signiﬁcant amount.

1.2

Solution Method
P modes are three-dimensional analogues of the standing waves (or modes) formed in an

organ pipe sounding at one of its natural frequencies. To produce a note in an organ, a soundwave propagates the length of the pipe, reﬂects from the opposite end and returns down the
tube. If the frequency of the sound wave corresponds to a mode of the pipe, the superpositionof the counter-propagating waves will resonate, producing a standing wave with the proper
frequency. Waves with other frequencies suﬀer destructive interference and will not produce atone. A similar mechanism occurs in the Sun, although complicated by spherical geometry, when
acoustic waves of the proper frequency and structure resonate to form standing waves — theaforementioned p modes. Unfortunately, this description of p modes, although mathematicallycorrect, causes diﬃculties with both visualizations and calculations. Ray theory is an alternative
approach which simpliﬁes each of these concerns at the expense of being only an approximation
to the true, modal, description.
Standing waves are the superposition of traveling waves which can mathematically, toa good approximation, be represented as plane waves — waves with the same amplitude andpropagation direction in all of space. Since plane waves travel in only one direction, they areconveniently described by a ray which points parallel to the direction of motion. Interaction ofwaves with the sides of any enclosure, such as the organ pipe considered above, will cause thisapproximation to break down. However, if the wavelength of the sound wave is small comparedto the dimensions of the enclosure, the eﬀects of the boundary can usually be ignored. Returning
to the solar case, if we consider a Sun where any inhomogeneities are large compared to the7wavelength of the p-mode (just as the dimensions of the organ pipe are large compared to thewavelength of the sound wave), then, by analogy, p modes are well approximated by rays.
In the spirit of this approximation, consider a beam of light incident upon a screen witha hole. When the hole is large compared to the wavelength of light, the ray approximation is avalid one and the light travels through the hole in a straight line. However, if the aperture issmall compared to the light’s wavelength, the light will diﬀract as it travels through the hole —a phenomenon which can only be explained by considering light as a wave. Within the domainof classical (non-quantum) mechanics, the wave picture will provide the correct description forboth situations. Nevertheless, it can be both conceptually and mathematically simpler to use
the ray approximation when it is valid. One of the reasons we employ ray theory is precisely
this simplicity. Recall that we seek the eigenfrequencies of p modes in a solar model containing
turbulent convection. The governing equations for the complete problem are quite diﬃcult —as is, by analogy, solving the problem of light traveling through a large hole using wave theory— particularly for the complex convective motions of the solar interior. By considering p modes
whose wavelengths are smaller than the dimensions of the convective motions, we can utilize thesimplifying ray approximation.
Until this point, we have been deliberatively vague about the nature of the computational
advantages provided by using the ray approximation. Mathematical analysis demonstrates thatthe equations describing the ray-like propagation of a p mode are cast in a special form called
Hamilton’s equations. Hamilton’s equations, which are equivalent to Newton’s laws of motion,
describe the time evolution of an object’s position and momentum and can be used, for instance,
to calculate particle trajectories. The correspondence between ray theory and particle propagation may not be surprising after considering the analogy presented above. When light passes
through a large aperture, its path can be treated with the ray approximation and is found to
be straight lines, just as a beam of particles would behave in a similar situation. Extensive
work has focused on systems obeying Hamilton’s equations, in particular on determining the
eigenfrequencies and modes of oscillation of such systems. By connecting the problem under8
consideration (p-mode oscillations in a convective Sun) with Hamiltonian dynamics, we candirectly apply these results to the system at hand.
Certain Hamiltonian systems possess quantities known as adiabatic invariants. In somesystems — including those with which we are concerned — the adiabatic invariants are foundto have special values when the system occupies a mode of oscillation. This property is known
as EBK quantization (Einstein 1917; Brillouin 1926; Keller 1958) and suggests that a proﬁtable
method for ﬁnding eigenmodes of complex systems may be a direct search for those rays whose
adiabatic invariants satisfy the quantization conditions. As we will show, this method, althoughconceptually useful, is unreliable for calculating the eigenfrequencies of suﬃciently complex
convective models.
However, adiabatic invariants possess another special property; as their name suggests,
they are unchanged when the system undergoes slow alterations. We use a technique known
as adiabatic switching (reviewed in Skodje & Cary 1988), which relies on both the invarianceand the quantization of the adiabatic invariants, to determine the eigenfrequencies of p modes
in convective systems. Beginning with a simple model for which theorists can calculate the pmode eigenfrequencies, we choose a ray corresponding to an eigenmode. Doing so ensures theray possesses correctly quantized adiabatic invariants. Next, we impose a slow alteration onto
the envelope in which the ray travels by, for example, gradually increasing the amplitudes of acomplicated set of convective motions. Since the adiabatic invariants remain constant duringslow changes, they continue to satisfy the correct quantization conditions during the alterationand the ray continues to correspond to an eigenmode. Although the invariants are constant,the eigenfrequency of the ray is not and it slowly changes from its initial value (known fromthe initial calculation) to a ﬁnal value, which is not known beforehand. Thus, by comparing
the two frequencies, we calculate frequency shifts caused by the convective motions. We willshow that this method successfully generates eigenfrequency shifts for complicated convectivemodels. We also demonstrate that these shifts oﬀer a possible solution to the discrepancy
between theoretically and observationally based solar models of the structure near the solar

9
surface.

1.3

Literature Review
The deleterious eﬀects of neglecting ﬂuid motions has been understood for some time andseveral attempts have been made to untangle the intertwined eﬀects of convection and stellar
pulsation (see Baker & Gough 1979, and references therein). Among the ﬁrst to consider the helioseismic implications was Brown (1984). While considering only vertical velocity perturbations
with a zero horizontal average, his work demonstrated that such perturbations always produce
eigenfrequency downshifts and suggested that these shifts could be of suﬃcient magnitude to
produce an observable eﬀect. These results are conﬁned to a special subset of photospheric models (those without perturbations in the horizontal direction or to thermodynamic quantities).
We conﬁrm them in this work and extend our results to more complicated ﬂows.
Others (Balmforth 1992; Rosenthal 1997; Rosenthal et al. 1998) have examined whether
the eigenfrequency discrepancies could be due to other eﬀects, notably turbulent pressure arisingfrom convective motions. Their analysis shows the discrepancy between theory and observation
can be reduced, although not enough to bring the two data sets into complete agreement. Furthermore, since the f mode (a surface gravity wave which propagates on the density interface
between the photosphere and the chromosphere) is a surface oscillation, it is not aﬀected by
the mechanism they consider even though its observed frequencies do not agree with theoretical
calculations. From this evidence, they acknowledge that the convective eﬀects which we consider
may also be important. They do conclude, however, that the resolution will come from consideration of modal physics, arising from the interaction between p modes and the solar structure,
rather than model eﬀects which arise from inaccuracies in the mean solar structure.
A series of papers (Murawski & Roberts 1993a,b; Murawski & Goossens 1993) addresses
the eﬀects of photospheric ﬂows and a chromospheric magnetic ﬁeld on the f mode. Although
the oscillation physics of f and p modes are quite diﬀerent, it is intriguing to note that they ﬁndthat random photospheric ﬂow ﬁelds lower the f -mode frequencies. Although they work in the10
incompressible ﬂuid limit, they are able to ﬁt the observational data (which show that the actual
f -mode eigenfrequencies are lower than theoretical predictions) with reasonable values for the
ﬂow speeds, and magnetic ﬁelds. Both Ghosh et al. (1995) and Duvall et al. (1998) extend thiswork, each showing that several features of the f -mode spectrum can be accounted for by theeﬀects of shearing velocity ﬁelds.
More applicable to this work are two papers (Lavely & Ritzwoller 1992, 1993) which
use quasi-degenerate perturbation theory to investigate the eﬀects of steady-state, large-scaleconvection on helioseismic linewidths and frequencies. In particular, they ﬁnd that convection
on the scale of giant cells (at least in one model) has a systematic eﬀect on the linewidths ofp modes and hence that such a signal can be used as a diagnostic tool. Although they also
ﬁnd a small eﬀect on p-mode frequencies, they consider only small to medium harmonic degrees( . 100, where  is deﬁned in §2.1.3) and hence their results do not apply to the same parameter
range as ours.
Gruzinov (1998) also uses an analysis based on perturbation theory, but with an emphasis
on ﬁnding analytic expressions for the eigenfrequency shifts. He derives a general result, in termsof integrals over the ﬂow ﬁelds, which he uses to examine turbulent shifts in the f mode. The
results include an unobservably small positive frequency shift for low-degree modes, a negative
frequency shift at high degrees, and a credible ﬁt with observational data. It is not clear,
however, if this method is computationally practical for complicated convective motions.
Also of note is a series of papers investigating acoustic waves in a structured medium(Zhugzhda & Stix 1994; Zhugzhda 1998; Stix & Zhugzhda 1998). They consider corrections
arising from a single sinusoidal perturbation to the sound speed and the vertical velocity, ignoring
any horizontal ﬂows. Again, they ﬁnd that the result is a downshift in the eigenfrequencies which
increases with both frequency and mode degree . Furthermore, for solar-like conditions theyﬁnd that velocity, and not sound speed, perturbations have a stronger eﬀect. As will be shown,
our results agree with these although we can extend our analysis to somewhat more complicated
ﬂows.

11
Our ﬁrst attempts at this problem were motivated by Gough (1993), who discusses EBK
quantization (see §2.4) as a means for determining eigenfrequency shifts. He assumes the convective eﬀects are weak, allowing him to use a perturbative expansion and express the frequency
shifts in terms of certain integral constraints. As will be shown below, we ultimately abandon
this technique in favor of adiabatic switching due to the ﬂexibility of the latter in handlingcomplicated convective motions. Our early work on this subject, portions of Chapters 2–4, hasbeen published elsewhere (Swisdak & Zweibel 1999).

Chapter 2Theoretical Justiﬁcation

This chapter provides the theoretical justiﬁcation for our approach to the problem. In§2.1 we discuss the theory of stellar oscillations and derive the governing diﬀerential equation.
In §2.2, we describe the WKB solution to this equation and show how it can be cast in terms of
the ray approximation. Section 2.3 outlines the general properties of Hamiltonian systems. The
use of EBK quantization to ﬁnd Hamiltonian eigenvalues is discussed in §2.4. Finally, in §2.5, wegive a synopsis of the method of adiabatic switching, treating both its theoretical underpinnings
and its application to the speciﬁc problem at hand.

2.1

Stellar Oscillations
The theory of stellar oscillations is a well explored ﬁeld with results ranging far wider

than is necessary for this work. The basic mathematics is well known and we present only a briefsummary of the germane topics; for more details, consult one of the many available references(Cox 1980; Gough 1993; Christensen-Dalsgaard 1997).

2.1.1

Fluid Equations
To derive the equations governing p-mode oscillations, we begin with the equations of ﬂuid

dynamics; for a derivation of these equations starting from Newtonian mechanics see Fetter &Walecka (1980). The mass conservation, or continuity, equation is∂ρ+ ∇ · (ρv) = 0,
∂t

(2.1)

13
where ρ and v are the ﬂuid density and velocity, respectively.1 Solar Reynolds numbers are &1012 (Brummell et al. 1995), implying that the viscous terms in the Navier-Stokes (conservationof momentum) equation are unimportant. We note that, even though they are numericallysmall, the dissipative eﬀects of viscosity can have signiﬁcant eﬀects on p modes, particularly inthe outer convection zone. We ignore them here, however, and take1∂v
+ (v · ∇)v = − ∇p + f ,
∂t
ρ

(2.2)

where p represents pressure and f is the sum of all body forces per unit mass. In this work, theonly body force we consider is gravity, so f = g where g is the acceleration due to gravity.2 The
gravitational acceleration also satisﬁes Poisson’s equation,∇ · g = ∇2 Φ = −4πGρ,(2.3)

for the gravitational potential Φ, where G is the Newtonian gravitational constant. A ﬁnal
relation is necessary to complete the set of equations. Beginning with the ﬁrst law of thermodynamics and assuming that the ﬂuid motions are adiabatic, it may be shown that the pressureand density are related by:dp
Γ1 p dρdρ== c2 ,
dt
ρ dt
dt

(2.4)

where the adiabatic exponent isΓ1 =∂ ln p∂ ln ρ


,

(2.5)

s

with the subscript s indicating that the derivative is taken at constant speciﬁc entropy. Equation(2.4) also deﬁnes the adiabatic sound speed c.
Justiﬁcation for the assumption of adiabatic ﬂuid motions in the context of p-mode oscillations comes from an examination of the relevant timescales. The motion of a ﬂuid element
may be approximated as adiabatic if the element undergoes several oscillations before heat is
1 Appendix A contains a table of important symbols used in this work.
2 Although Lorentz forces arising from magnetic ﬁelds could also be included, their eﬀects in the solar convection zone are usually signiﬁcant only in active regions. We disregard them here.

14transferred between it and the surroundings. In other words, the oscillation timescale must be
shorter than the timescales for radiation and convection, the relevant heat transfer mechanisms
in the Sun. The canonical period of a solar p-mode is 5 minutes. In the solar convection zone
the radiation timescale is ≈ 102 years, which is signiﬁcantly longer than an oscillation period.3
Convective timescales, however, are much shorter. Observations show that granules have lifetimes of ≈ 10 minutes. Larger convective motions such as mesogranules (≈ 1 hour lifetimes)
and supergranules (≈ 1 day lifetimes) have been detected and giant cells persisting for a solar
rotation period of 30 days may also exist. Granules, with lifetimes of ≈ 1 − 2 oscillation periods,
clearly violate the assumption of adiabatic ﬂuid motions. But, as will be shown below, p modes
rarely occupy the same regions as granules, hence justifying our dismissal of their non-adiabaticeﬀects.

2.1.2

Perturbations
The amplitudes of solar oscillations are small in comparison to other ﬂuid motions in

the Sun. The ﬁrst discovery of p modes (Leighton et al. 1962) showed the solar surface tobe tiled with oscillating ﬂuid elements having a characteristic velocity amplitude of 1 km s−1 .
Later observations found that these oscillations are the superposition of ≈ 107 p modes withindividual amplitudes of ≈ 1 m s−1 . For comparison, the adiabatic sound speed near thephotospheric surface is c ≈ 10 km s−1 ; in the interior the sound speed rapidly rises while theoscillation amplitude falls. The resulting Mach numbers, ≈ 10−4 at the surface maximum, imply
p modes are linear oscillations, an assumption we carry throughout this work.
We seek oscillations about a reference state satisfying equations (2.1)–(2.4) and certain(unspeciﬁed) boundary conditions. As p-mode oscillations are small, we use linear perturbation
theory. Assuming the existence of a static,4 spherically symmetric reference model satisfying
equations (2.1)–(2.4), we expand each variable around its equilibrium state. Any variable Q
3 The radiation timescale drastically decreases at the photosphere. However, processes in this region cause
diﬃculties for several reasons, as will be seen below, and are ignored in this work.
4 Direct observation supports this assumption as the mean solar structure does not vary on the timescales inquestion.

15
becomesQ(x, t) = Q0 (x) + Q (x, t),(2.6)

where x is the position vector, Q0 is the variable evaluated in the reference model and Q is
an Eulerian perturbation. (For a full description of Eulerian and Lagrangian perturbations, see
Appendix B.)
To continue, we replace the variables in equations (2.1)–(2.4) with their linearized counterparts, dropping the 0 subscripts from reference quantities. The continuity equation becomesρ + ∇ · (ρδx) = δρ + ρ∇ · (δx) = 0,

(2.7)

where δ represents a Lagrangian perturbation, while the momentum equation isρ∂ 2 δx= −∇p + ρg + ρ g.
∂t2

(2.8)

From the adiabatic approximation of equation (2.4), the density and pressure perturbations areconnected throughδp = c2 δρ.

(2.9)∇2 Φ = ∇ · g = −4πGρ .

(2.10)

Finally, Poisson’s equation becomes

2.1.3

Separation of Angular Variables
Since the reference state is spherically symmetric, it proves helpful to consider the angularbehavior of the independent variables in spherical coordinates (r, θ, φ). After rewriting the
perturbed form of Poisson’s equation (2.10) as
1 ∂
r2 ∂r



∂Φ
r2
+ ∇2h Φ = −4πGρ ,∂r(2.11)

we see that derivatives with respect to the angular variables θ and φ appear solely as part of
the horizontal Laplacian
∇2h =

1
∂
2
r sin θ ∂θ



∂
1
∂2
sin θ
+ 2 2
.
∂θ
r sin θ ∂φ2(2.12)

16
Assuming that ρ and Φ have the same angular dependence, separation of the angular coordinates, leads to the well known result that the angular eigenfunctions are the spherical harmonics,Ym (θ, φ). To be precise, we haveρ (x, t) = ρ (r, t)Ym (θ, φ)(2.13)

Φ (x, t) = Φ (r, t)Ym (θ, φ),(2.14)

where , a non-negative integer, is the degree of the harmonic, while m is restricted to the range− ≤ m ≤  and is called the azimuthal order. Given equations (2.13)–(2.14), we can show thatthe spherical harmonics are also the angular eigenfunctions of the other scalar perturbations.
Although the spherical harmonics are familiar from other problems of mathematicalphysics, it is useful to discuss a few of their properties as they relate to p-mode structure.
A full solution of the eigenvalue problem produces radial eigenfunctions ρ (r, t) demarcated byan integer n; the full eigenfunctions ρ (x, t) are then denoted by three integers: n,  and m. The
radial eigenvalue roughly measures the number of oscillations in the radial direction,  and mcharacterize the oscillations on a spherical surface. Hence as n and  increase, the radial andazimuthal wavelengths decrease.
The solar acoustic structure traps p modes in a radial cavity. Within the cavity the mode
resonates and produces a standing wave, while outside of the cavity the mode is evanescent. The
radial location of the cavity depends on the mode numbers, roughly increasing linearly with and decreasing linearly with n. For n = 1,  > 15 the cavity of a p mode is conﬁned within thesolar convection zone (r & 0.7R ) and oscillations with n = 1,  > 80 have a minimum radiusof 0.95R.

2.1.4

Derivation of a Wave Equation
Equations (2.7)–(2.10) may be recast in a suggestive form through a derivation of Gough(1993) based on the earlier work of Lamb (1932). Deﬁne the variableχ ≡ ∇ · δx = −δρδp=− 2 ,
ρc ρ(2.15)

17
where the last two equalities follow from equations (2.7) and (2.9). Substituting equation (2.15)
into the perturbed momentum equation (2.8) yields
ρ∂ 2 δx= −∇(ρδx · g − ρc2 χ) + ρg − ρχg − (δx · ∇ρ)g,∂t2(2.16)

where we have made use of the equation of hydrostatic equilibrium, ∇p = ρg, and equation(B.10). Dividing through by ρ and rearranging, we ﬁnd that∇ρ1∂ 2 δx∇ρ= ∇(c2 χ + δx · g) + c2 χ+ δx · g+ g − χg − (δx · ∇ρ)g.
2
∂t
ρ
ρ
ρ

(2.17)

The curl of the equation of hydrostatic equilibrium shows that g and ∇ρ are parallel and allows
us to cancel the third and sixth terms on the right-hand side. After a small rearrangement, wehave∂ 2 δx∇ρ= ∇(c2 χ + δx · g) + (c2− g)χ + g .
∂t2
ρ(2.18)

Finally, recall that the reference state was assumed to be spherically symmetric. Denoting
the radial unit vector as r̂ gives g = −g(r)r̂ and deﬁning the radial component of the Lagrangian
displacement asξ ≡ r̂ · δx

(2.19)∂ 2 δx∇ρ= ∇(c2 χ − gξ) + (c2+ gr̂)χ + g .
∂t2
ρ(2.20)

produces

We can proceed further, using equations (2.10), (2.15), and (B.10) to substitute for gand obtain an equation for δx only in terms of the parameters of the reference model. A more
fruitful approach comes from making several well motivated and simplifying assumptions which
we discuss below.

2.1.4.1

Cowling’s Approximation

In making Cowling’s approximation we assume that Φ , the perturbation to the gravitational potential due to the oscillations, is negligible compared to the density variation ρ and18thus that we can neglect g in equation (2.20). First demonstrated by Cowling (1941), it arises
from a consideration of the solution to the perturbed Poisson’s equation (2.10) in integral form:Φ (r) =

  ∞

 r
4πG
r
ρ (x)x1− dx + r−1−ρ (x)x+2 dx .
2 + 1r
0(2.21)The perturbed potential Φ may be safely neglected when the spherical harmonic degree is large,
as both integrals contain terms which rapidly decay with increasing . In addition, Φ can safely
be neglected in modes where ρ rapidly changes sign, which happens when n is large. Since
factors we discuss in §2.2 restrict our investigation to modes with  & 100, we hereby set g = 0
in what follows.

2.1.4.2

Further Approximations

Since we desire oscillatory solutions, we expect the governing equation to be similar inform to the standard wave (Helmholtz) equation for a scalar variable. After making Cowling’s
approximation, equation (2.20) becomes∂ 2 δx= ∇(c2 χ − gξ) + βχr̂.
∂t2(2.22)

where
β≡

c2 dρ
c2 N 2
+g =−
ρ dr
g

(2.23)and N is the Brunt-Väisälä frequency. Modes of high spherical degree occupy a cavity with asmall radial extent near the top of the convection zone. As a simple, but useful, approximation
we take gravity to be constant throughout this region and simplify the governing equations. To
be precise, we assume the variation of gravity on the vertical scale of the oscillations is smalland hence that all derivatives of g with respect to r vanish.5
After applying the divergence operator ∇· to equation (2.22), we have∂2χ2βχ
= ∇2 (c2 χ − gξ) + r̂ · ∇(βχ) +.
∂t2
r(2.24)

5 To be consistent we should work in the plane-parallel limit r → ∞. However for now we keep terms
depending inversely on r noting that they, by chance, disappear from the ﬁnal result of equation (2.30). For aconsistent treatment see Appendix C.

19
To eliminate the variable ∇2 ξ, apply the operator r̂ · (∇ ×∇× to equation (2.22). Using various
vector identities we obtain2βχ∂2(r̂ · ∇χ − ∇2 ξ) = −∇2h (βχ) − 2 ,2∂tr(2.25)

where ∇2h is the horizontal Laplacian operator. Combining equations (2.24) and (2.25) gives
∂2
∂4χ
−
∂t4
∂t2



2βχ
2βgχ
∇2 (c2 χ) + r̂ · ∇ [(β − g)χ] ++ 2 + g∇2h (βχ) = 0.
rr

(2.26)

Since this equation contains odd spatial derivatives of χ, it is not quite in the desired Helmholtzlike form. The substitution Ψ = c2 ρ1/2 χ = −ρ−1/2 δp eliminates the odd-order derivatives and
gives




2βχ ∂ 2 Ψ
2βgc2
∂4Ψ
2 2
2
2 2 2
−
c
∇
−
ω
−
−
c
N
∇
−
Ψ = 0.
c
h
∂t4
r
∂t2
r2

(2.27)

The acoustic cutoﬀ frequency ωc is deﬁned by
ωc2 ≡

c2
4Hρ2



4Hρ
1 − 2r̂ · ∇Hρ −r(2.28)whereHρ ≡ −

d ln ρ
dr−1(2.29)

is the density scale height.
Solar convection is so eﬃcient that the underlying reference state is nearly adiabaticallystratiﬁed. Non-adiabatic stratiﬁcation is a signiﬁcant eﬀect only in the outer portions of theconvection zone where other aspects of our analysis (primarily the assumption of adiabatic ﬂuid
motions) also break down. Thus N 2 (and therefore β) are ≈ 0 in the bulk of the convection
zone, the region of interest in this analysis. Neglecting the appropriate terms in equation (2.27),
we arrive at a modiﬁed wave equation



∂2
2
2 2
+
ω
c Ψ − c ∇ Ψ = 0.
∂t2(2.30)

The analysis leading to equation (2.30) required several assumptions. In roughly theorder of occurrence, they are

20(a) Non-dissipative, adiabatic ﬂuid motions(b) Linear perturbations(c) Spherically symmetric, non-magnetic, static reference state(d) Cowling’s approximation(e) Constant gravity(f) Adiabatic stratiﬁcation.
Appendix C discusses the ramiﬁcations of relaxing restrictions (e) and (f). Note that the derivation of equation (2.30) did not assume a plane-parallel geometry, only that variations in thegravitational acceleration are small over the region of interest. Below we will adopt a planeparallel geometry, the only eﬀect of which is to eliminate the ﬁnal term in equation (2.28).

2.2

Asymptotic Approximations
The ultimate goal is to determine the eigenfrequencies of a solar model which includes

convective motions. In principle this means solving equation (2.30) in a domain where c andωc are complicated functions of position — a diﬃcult computational problem which has thepotential to obfuscate the underlying physics within a morass of numerical work. Instead, weseek the eigenvalues for a particular limit of the general problem, the short-wavelength limit
where the lengthscales of the oscillations are much shorter than those of the other varying
quantities. Since, for example, a p mode with an  = 200 has a typical horizontal length scale ofR / = 3.5 Mm, our analysis will only be valid for convective motions of mesogranular scales
and larger. Of course, as granules already violate the assumption of adiabatic ﬂuid motions usedin the derivation of equation (2.30), our mathematical formalism has already implicitly excluded
these motions. On the other hand, granulation near the solar surface is the leading candidate
for the source of p modes and hence would be expected to leave an imprint on the oscillation
frequencies. Our analysis is insensitive to any such signal. The short-wavelength approximation,

21
also known as geometrical acoustics or optics, is equivalent to the WKB method used in theasymptotic analysis of diﬀerential equations.
Our assumption of small-wavelength modes allows us to use the well known analogy
between geometrical optics and particle motion to cast the problem in a Hamiltonian formalism.
A method analogous to EBK semi-classical quantization then provides a relatively simple routeto the desired eigenvalues. However, we will demonstrate that this method is diﬃcult, if not
impossible, to implement for the complicated domains under consideration. However, due to the
Hamiltonian nature of the system, a related technique known as adiabatic switching produces
eigenvalues even for complex systems.
Consider, for the moment, an unrealistic solar model where ωc = 0 and c is a constant.
Equation (2.30) reduces to the Helmholtz equation, which has plane-wave solutionsΨ = Ψ0 e±i(k x−ωt) ,(2.31)

where Ψ0 is the constant amplitude of the wave and the wavenumber k and frequency ω arerelated by the dispersion relation ω 2 = c2 k 2 . These plane waves have the distinctive property
that their propagation directions and amplitudes are the same for all space. In this case, we canactually ignore the wave nature of equation (2.31) and treat the solution as a one-dimensionalray, propagating perpendicular to the wavefronts.
Now, consider instead a system where the lengthscale over which equilibrium quantities(such as c, p and ρ) vary is much longer than the wavelength of a mode. Letting the ratio of theselengthscales be denoted by Λ (∼ Hρ |k|) — a constant, large parameter used for bookkeeping
purposes — we write the solution to equation (2.30) as a plane wave with varying amplitude,A, and phase, ϕ (assumed to be real):Ψ = A(x, t)eiΛϕ(x,t) .

(2.32)

Although we have not made any approximations in assuming this form for the solution, thepresence of Λ suggests that Ψ rapidly oscillates in comparison to the background state. Thus, it
is reasonable to think of equation (2.32) as describing waves which are locally planar and hence

22
amenable to a ray description, albeit one where the rays undergo changes in both direction andamplitude during propagation. If, as is the case for p modes, an eigenmode is trapped withina cavity the wavenumber of the solution disappears at the boundaries. There, Λ → 0 and it
is impossible to satisfy the locally planar condition. This diﬃculty is merely the well known
breakdown of the WKB approximation near the classical turning points of a trajectory.
Using the ansatz of equation (2.32) and following the work of Gough (1993), we expand
equation (2.30) and equate powers of Λ. This process, equivalent to the WKB approximation,
gives the leading equation


∂ϕ
∂t

2
−

 ω 2cΛ

2− c2 |∇ϕ| = 0.

(2.33)In a slightly diﬀerent form, this is the eikonal equation of geometric optics. Next, making theanalogy between equations (2.31) and (2.32), we deﬁne the local frequency ω̃ and wavenumberk asω̃(x, t) ≡ −Λ

∂ϕ
∂t

and

k(x, t) ≡ Λ∇ϕ.

(2.34)

These identiﬁcations allow us to write equation (2.33) in the form of a dispersion relation:1

ω(x, k, t) = (c2 k 2 + ωc2 ) 2 ,(2.35)

where ω diﬀers from ω̃ only in the inclusion of k as an independent variable. Furthermore,equation (2.34) also implies that∂k+ ∇ω̃ = 0,∂t(2.36)

 

dk
∂ω dx
+ ∇ω +−· ∇k = 0,dt∂kdt(2.37)

from which we ﬁnd thatdwhere dt
is the total derivative with respect to time along a ray path. As the grouping of theterms suggests, this equation is satisﬁed when each expression is identically zero. The result isa set of ﬁrst-order diﬀerential equations which describe the trajectory of a ray:
dk
∂ω
=−
dt
∂x

and

dx
∂ω
=
.
dt
∂k

(2.38)

23
Equations (2.38) are Hamilton’s equations for a Hamiltonian ω in terms of the canonicalpositions qi = xi and momenta pi = ki . The dispersion relation of the ray, ω(x, k, t), is identiﬁed
as the functional form of the Hamiltonian. Using the WKB approximation we have recast theoriginal problem; instead of solving the eigenvalue problem presented by the partial diﬀerential
equation (2.30), we seek the eigenfrequencies of the Hamiltonian given by equation (2.35). This
reworking comes at a price since once we pass to a ray description we lose information concerning
the structure of the eigenmode. However, we are primarily concerned with the eigenfrequencies
which, as we will show, can be determined in a straightforward manner for the Hamiltoniansystems under consideration.
While the appearance of Hamilton’s equations may appear serendipitous, it has been
known since the time of Hamilton himself that the eikonal equation is equivalent to the HamiltonJacobi equation of classical mechanics. As an aside, it is well known in quantum mechanics(Goldstein 1965) that high-frequency solutions to the wave equation (in other words, solutions
which oscillate rapidly compared to any background variation) follow ray-like trajectories similar
to classical particles. In that case, the semi-classical limit, the inverse of Planck’s constant, h−1 ,plays an analogous role to our parameter Λ.
Finally, the dispersion relation of equation (2.35) is not a fully satisfactory descriptionof ray motion in a convective domain. While it does describe the trajectory of a ray in a Sun
where c and ωc are arbitrary functions of position, any method which wishes to treat realistic
convective structures must also consider the possibility of advective motions. This oversight isa result of our initial description of the problem: equations (2.7)–(2.10) do not allow for ﬂuid
motions in the reference state and hence they do not appear in the ﬁnal dispersion relation. To
account for this discrepancy, we will add a term to equation (2.35) in §3.2 which takes the formof a Doppler shift and accounts for the eﬀects of advective motions.

24

2.3

Hamiltonian Systems
The identiﬁcation with Hamiltonian mechanics allows us to take advantage of the exten-

sive results available for Hamiltonian systems. Before discussing how to ﬁnd the eigenfrequencies
for such systems, we review several results from Hamiltonian theory which will be of use in later
sections. More details on these subjects can be found in a number of texts on classical mechanics
and Hamiltonian systems (for example Goldstein 1965; Arnold 1978; Lichtenberg & Lieberman
1983; Tabor 1989).
A general Hamiltonian H(q, p, t), satisﬁes Hamilton’s equations
dq
∂H
=
dt
∂p

and

dp
∂H
=−
dt
∂q

(2.39)

for the canonical position vector q and canonical momentum vector p. If the Hamiltonian
is cyclic in one of the canonical positions qi (i.e., the Hamiltonian does not explicitly dependon that coordinate), then the corresponding canonical momentum pi is also a constant of the
motion:∂H=0∂qi

implies

dpi
= 0.
dt

(2.40)

Furthermore, the time derivative of H is given bydH∂H∂H dq ∂H dp
∂H
=
+
+
=
,
dt
∂t
∂q dt
∂p dt
∂t

(2.41)

where the ﬁnal equality follows from equation (2.39). Hence, for systems where the Hamiltonian
is not explicitly time-dependent, the Hamiltonian itself is a constant of the motion. Finally, we
note that if a Hamiltonian of degree n does contain an explicit time-dependence, it is possible
to rewrite it as a time-independent Hamiltonian of degree n + 1 in which the time and the old
Hamiltonian become elements of the canonical position and momentum vectors, respectively.
Certain special Hamiltonians, named integrable systems, have equal numbers of independent constants of the motion and degrees of freedom. Such systems permit a canonicaltransformation (one which preserves the Hamiltonian nature of the system) in which all of the

25
new canonical positions, usually termed angle variables and denoted Θk , are cyclic. The new
canonical momenta are the action variables Ik , and the transformed Hamiltonian H(Ik ) satisﬁes
a certain set of Hamilton’s equations:
∂H
dIk
=−
=0
dt
∂Θk

dΘk
∂H
=
= νk (I),
dt
∂Ik

and

(2.42)

where νk is the characteristic frequency for motion in the kth coordinate.6 Since the angles are
cyclic, each action variable is one of the n independent constants of the motion.
Thus, a Hamiltonian system with n degrees of freedom deﬁnes a 2n-dimensional phasespace, but integrability ensures that the phase-space trajectories are conﬁned to (2n − n) = ndimensional surfaces. It can be shown both that this surface is topologically equivalent to ann-dimensional torus and that the action variables may be written as
Ik =

1
2π


p · dq,(2.43)

Ck

for k = 1, . . . , n. The curves Ck are topologically independent closed paths on the n-dimensional
torus. An integrable system is conﬁned to a torus by its initial conditions (which set the values of
the invariant action variables) and does not stray from this torus as long as the system remains
integrable.
A concrete example is provided by a two degree of freedom Hamiltonian; in particular,
we shall consider equation (2.35). In the appropriate notation, the Hamiltonian (ω) is writtenasH(q1 , q2 , p1 , p2 ) =c2 (p21 + p22 ) + ωc2(2.44)

where the sound speed c and the acoustic cut-oﬀ frequency ωc are functions of the canonicalpositions q1 and q2 . The system has two degrees of freedom and thus a four-dimensional phasespace. As the Hamiltonian is independent of time, H itself is a constant of the motion. Furthermore, if we specialize to a system where c and ωc are both independent of one of the position
coordinates, say q1 , then the corresponding wavenumber p1 is also a constant of the motion. The
6

These are usually represented as ωk . To avoid confusion with the acoustic cutoﬀ frequency ωc and the
Hamiltonian ω, we adopt this non-standard notation.

26
existence of two constants ensures the integrability of the system and restricts phase-space motions to a two-dimensional surface, topologically equivalent to a familiar two-dimensional torus.
In fact, a family of invariant tori exist, each individual torus corresponding to a particular pair(H, q1 ). As the tori do not intersect, initial conditions constrain the system to a particular torus
for all time. To write the Hamiltonian of equation (2.44) in action-angle coordinates requires
the speciﬁcation of the functional form of c and ωc and would lead us astray at this time. We
return to the subject in Chapter 3.
General Hamiltonian systems, however, are usually non-integrable and thus their trajectories are not conﬁned to n-dimensional phase-space tori. However, special attention has been
given to systems which are nearly-integrable in the sense that they may be written as
H = H0 (I) + H1 (I, Θ),(2.45)

where I and Θ are vectors of the action and angle coordinates, respectively, and  is a small
parameter. The fate of invariant tori in non-integrable systems is addressed by the KAM(Kolmogorov-Arnold-Moser) theorem. We will return to the implications of this theorem inlater sections, but roughly it states that for suﬃciently small values of , most invariant tori
are preserved. In other words, for small perturbations from integrability, most sets of initialconditions remain on invariant tori. However, the destroyed and invariant tori are intermingled
in phase space and as the strength of the perturbation increases, so does the density of thedestroyed tori. For a strong enough perturbation, no invariant tori remain.

2.4

Quantization Conditions
We now return to the question of determining the eigenfrequencies of a convective solarmodel. Before deriving the p-mode quantization conditions, we brieﬂy digress with a discussionof semi-classical EBK quantization as it is formulated in quantum mechanics. Then, we examine
a simple solar model to derive insight into the behavior of the ray solutions. Finally, we determine
the p-mode quantization conditions and discuss their implications.

27
2.4.1

EBK Quantization
Prior to the development of wave and matrix mechanics, quantum theory postulated

that the correct description of the quantum world arose from quantizing the classical action
variables. Perhaps the most familiar example is Bohr’s derivation of the hydrogen spectra fromthe quantization of orbital angular momentum. The general result is displayed in the BohrSommerfeld-Wilson quantization condition:Ij =pj dqj = nj h,(2.46)

Cj

where h is Planck’s constant, nj is an integer, and Cj is the closed contour associated withthe motion in the jth degree of freedom. Disregarding the factor of 2π, this action is slightly
diﬀerent than that of equation (2.43) since it assumes that motions in diﬀerent coordinates are
separable. But, it was quickly realized that the separation procedure is not unique: evaluation in
diﬀerent coordinate systems produces diﬀerent quantization conditions. As an additional defect,
the zero-point energy (nj → nj + 1/2) found in, for example, the quantum harmonic oscillator,
could not be accounted for with equation (2.46) and was added as an empirical correction.
Einstein (1917) proposed a solution to the ﬁrst of these problems by noting that since thequantity p · dq, unlike the individual terms pj dqj , is invariant under coordinate transformations,the action variable to be quantized should be that given in equation (2.43). Soon thereafter,
modern quantum theory arose and it was discovered that Schrödinger’s wave equation reduced
to the Hamilton-Jacobi equation (and hence classical mechanics) in the semi-classical limit of
h → 0. Brillouin (1926), working in this limit, conﬁrmed Einstein’s supposition concerning theform of the quantized action variables by requiring that the wavefunction solutions be singlevalued. Finally, Keller (1958) added a term to the quantization condition which accounts forphase loss at caustic surfaces.7 The result is EBK semi-classical quantization which states that
7 Caustics occur when rays converge onto a lower-dimensional surface. A telescope focus, where a threedimensional bundle of rays coalesces to a point, is one example.

28
for quantum phenomena the classical action variable Ik of equation (2.43) is quantized:
Ik =

1
2π




sk 
p · dq = ~ nk +,4
Ck

(2.47)

where nk is an integer and sk , Keller’s contribution, is the number of encountered caustics. It
is interesting to note that the ﬁnal term can naturally account for the zero-point energy absent
from the Bohr-Sommerfeld-Wilson quantization condition.

2.4.2

A Simple Physical Model
Before deriving the appropriate quantization conditions for p modes, we ﬁrst examine a

simple system for some physical insight into the procedure. In §2.2 we demonstrated that, inthe small-wavelength limit, the solutions to equation (2.30) are rays. For now we work in plane
polar coordinates (r, φ), use the simpliﬁed dispersion relation ω 2 = c2 k 2 , and assume that c isazimuthally symmetric, c = c(r). The equations governing ray motion are found from equation(2.38):
dr
c2 kr
=
dt
ω

and

r

dφ
c2 kφ
=
,
dt
ω

(2.48)and the raypaths are solutions to
rkr
dr
=
=r
dφ
kφ

ω2
− 1.
c2 kφ2(2.49)In order to proceed, we need to specify several quantities: the form of c(r), a frequency ω (aconstant of the motion as c is independent of time), an azimuthal wave number kφ (anotherconstant as c is independent of φ), and initial conditions. For a given c, certain combinationsof ω and kφ will yield eigenmodes (we discuss below how to ﬁnd these combinations), however
raypaths still exist even for non-modal values of those parameters. Choosing a reasonable form
for the sound speed,8 as well as values for ω and kφ , results in the raypath shown in Figure 2.1.
Some general features of the raypath, independent of the initial conditions or speciﬁcvalues of ω and kφ , are noteworthy. A raypath propagates within a well deﬁned radial cavity
8 Although the actual form is not important, we use a polytrope with c2 ∝ r. Polytropes are discussed furtherin §3.1

29

Figure 2.1: Acoustic raypaths in a two-dimensional Sun. The caustic at the lower turning pointis clearly visible. Since the dispersion relation does not include an acoustic cutoﬀ frequency, theupper caustic coincides with the surface. As will be seen later, inclusion of ωc eliminates thecusps at the upper caustic, producing raypaths which are everywhere diﬀerentiable.

30bordered by two caustic surfaces at which kr changes sign. The quadratic dependence of thedispersion relation on the wavenumber implies that at any point within the cavity two rays, withkr equal and opposite in magnitude, can exist. However, it will later prove useful to think of therays as propagating in a slightly more complicated geometry in which two sheets are connected
at the caustics, one in which kr is restricted to positive values, the other in which kr is negative.
This domain is periodic in both the r and φ directions and is equivalent to the surface of a torus(unrelated, however, to the phase-space tori of §2.3). The construction of this new domain canbe extended to higher-dimension systems, although visualization is more diﬃcult.
We have not discussed how the values of ω and kφ in Figure 2.1 were picked. In fact, they
were chosen so that the ray corresponds to a p mode although, again, we have not yet indicated
how this was done — indeed, the exact point of this work is to determine ω for arbitrary solar
structures. Note, however, that even though it corresponds to an eigenmode, the ray does not
close on itself. Instead it will eventually ﬁll the entire cavity. Although the raypath itself gives
no visual indication of whether it corresponds to an eigenmode, there is a technique similar to
EBK quantization which allows one to make such determinations.

2.4.3

P Mode Quantization Conditions
We now consider how to determine the eigenvalues of equation (2.30) in the short-wavelength limit. The result will be almost identical to the EBK condition of equation (2.47),
although the application to non-quantum systems was probably ﬁrst discussed in Keller & Rubinow (1960). Our starting point is the wavefunction form of the solution as written in equation(2.32)
Ψ = A(x, t)eiΛϕ(x,t) ,(2.50)

with ϕ real. As seen in the discussion of §2.4.2, the solution is a sum of multiple waves with
opposing wavenumbers, however the forthcoming argument applies to each term individuallyand hence we work with a single term. Although this wavefunction is a classical solution, it is31
analogous to a short-wavelength solution to the Schrödinger equation Ψ = A exp(ih−1 S) where
S is the action and h−1 is equivalent to Λ. This is the basis for the strong similarity between
the EBK result and the result of this section.
Since the domain is periodic, the wavefunction Ψ is constrained to be single-valued. No
such restriction is placed on either the amplitude A or the phase ϕ. If we imagine traversing aclosed circuit in the domain, the requirement that the solution be single-valued impliesΛ(∆ϕ) = 2πn + i(∆ ln A)(2.51)

for any integer n, where ∆ represents the diﬀerence accumulated along the circuit. However, we
can quickly rewrite this condition by noting thatΛ(∆ϕ) = Λ

∇ϕ · dx(2.52)and thus, from equation (2.34),
k · dx = 2πn + i(∆ ln A).

(2.53)

Keller (1958), in the context of EBK quantization, ﬁrst evaluated the ﬁnal term of this equation,
realizing that, as in optics, the phase of the amplitude is retarded by π/2 whenever a ray
encounters a caustic. With this contribution, the condition becomessk · dx = 2π(n + ),4

(2.54)

where s is the number of caustics encountered along the path.
Although the quantization condition of equation (2.54) is valid around any closed curvein the domain, this does not lead to an inﬁnite number of quantum conditions. Instead, there
are as many independent curves, and hence quantization conditions, as there are degrees of
freedom in the system.9 In the example considered in §2.4.2, evaluation of equation (2.54) along
any two independent curves will yield two quantization conditions. Of course, some curves are
more equal than others: wise choices in that case would be a curve at constant r and another
9 Topologically, two curves are independent when they cannot be continuously deformed into each other.

32
at constant φ. So, the complete quantization conditions arek · dx = 2π(nk +
Ck

sk
)
4

(2.55)

where the Ck are independent curves and nk and sk are integers.
From our previous work relating k and x to the canonical momenta and positions, we
can rewrite (2.55) in the suggestive form
1
Ik =
2π


k · dx = nk +
Ck

sk
,
4

(2.56)

which both connects the quantization condition to the realm of Hamiltonian mechanics discussedin §2.3 and shows that it is analogous to the EBK result of equation (2.47).
Unfortunately, while (2.56) appears to give a general method for determining eigenfrequencies, its use is limited to integrable systems. As Einstein (1917) recognized, in the context
of EBK quantization, for non-integrable systems the curves Ck of equation (2.56) no longer exist
and evaluation of the action variables becomes, in a strict mathematical sense, impossible. This
restriction severely limits the convective structures which can be studied with this technique.
Recall that Hamiltonian systems with two degrees of freedom, such as that of §2.4.2, require twoconstants of the motion to be integrable — in that case, ω and kφ . Furthermore, each of these
constants reﬂects a symmetry in the governing Hamiltonian, time-independence in the ﬁrst caseand azimuthal symmetry in the second. Although most convective ﬂows can be assumed to bestationary on the timescale of p-mode oscillations, they do not have azimuthally-independent
structures. Hence, kφ is not a constant of the motion in such ﬂows, the system is non-integrable,
and equation (2.56) is an unsuitable method for determining the system’s eigenfrequencies.

2.4.4

Surfaces of Section
As we will see, adiabatic switching addresses these diﬃculties. Before its introduction in§2.5, we brieﬂy consider an aesthetically pleasing interpretation of the quantization conditionsof the previous section. Although Noid & Marcus (1975) used this approach with some successto ﬁnd the energy levels of an anharmonically coupled pair of oscillators, it does not extend to33
non-integrable systems and is hence unsuitable for our needs. However, it gives the problem ageometric interpretation and, in addition, was the path we took on our ﬁrst explorations of thesubject. Although it can be applied to higher order systems, this method is most useful for asystem with two degrees of freedom. We will avoid a formal treatment and instead concentrateon the system of §2.4.2; Tabor (1989) has a more general discussion.
Instead of treating equation (2.55) as an integral equation which must be solved to ﬁnd theeigenfrequencies, we interpret it as a geometric statement. The Hamiltonian of §2.4.2 possesses
two degrees of freedom and hence a ray propagates in a four-dimensional (r, φ, kr , kφ ) phasespace. However the constancy of ω and kφ restrict the trajectory to a two-dimensional surface,
topologically equivalent to a torus. Imagine the trajectory intersecting a plane at an arbitraryvalue of φ and making a record of the values of r and kr at each intersection. After suﬃcienttime, a surface of section will be traced: a closed curve which for our case can be found fromsolving the dispersion relation
kr2 =

ω2
− kφ2 .
c2

(2.57)

A similar curve arises from considering the values of kφ and φ at the intersection of the torus
with a plane at constant r. Equation (2.55) is thus a statement concerning the area traced out
by these curves: for eigenmodes the areas have certain quantized values. Rays not corresponding
to modes trace out curves enclosing other areas not given by equation (2.55). Hence, one candetermine the eigenmodes of a system by ﬁnding those rays which produce surfaces of section
with the proper areas.
For integrable systems with two degrees of freedom, the trajectories are restricted to twodimensional surfaces in phase space and hence the surfaces of section consist of smooth curves
enclosing well deﬁned areas. As a system is perturbed from integrability, the KAM theorem
ensures that increasing numbers of invariant tori are destroyed. Rays whose initial conditions
correspond to the destroyed tori are then free to wander in a three-dimensional region of phasespace. A two-dimensional slice through a three-dimensional region of phase space does not34
produce curves which enclose well-deﬁned areas and hence makes the geometric interpretationof equation (2.55) unclear. The surfaces of section acquire a characteristic “fuzzy” appearance,
an early example of which is seen in the work of Hénon & Heiles (1964). This is a consequenceof the assertion made in §2.4 that the quantization conditions are undeﬁned for non-integrable
systems.

2.5

Adiabatic Switching
Ehrenfest (1917) was among the ﬁrst to propose the concept of adiabatic switching,although his work has since been superseded by modern treatments. In essence, his hypothesis
was that if a system is changed in a reversible, adiabatic way, then allowed motions will smoothlytransform to allowed motions. Perhaps the most familiar example is provided by a smallamplitude pendulum. If the length of a pendulum is slowly increased, the energy, frequency,
and oscillation amplitude will each decrease. However, the ratio of the energy to the frequency
remains constant. In fact, this ratio remains constant for any perturbation, such as changing
the gravitational acceleration or the mass of the bob, as long as it is performed on a timescale
longer than the dynamical period of the system.
For a simple harmonic oscillator, of which the small-amplitude pendulum is an example,
the ratio of the energy to the frequency is the classical action variable of §2.3. In a rough sense,
an adiabatic change smoothly deforms the invariant tori between those of the initial and thoseof the ﬁnal state. This property is of particular interest if a system has eigenmodes, for in thiscase the action variables — at least in the semi-classical limit — are quantized. If we start thesystem in an eigenmode for which the eigenfrequencies are known and allow it to evolve as thestrength of a perturbation is increased, eventually we are left with a new system for which the
eigenfrequencies were not known beforehand. But, due to the adiabatic nature of the transition,
the action variables remain quantized and hence the original eigenmode has slowly relaxed into
an eigenmode of the new system. Although we lend a degree of formalism to this qualitativedescription below, we do not present a complete derivation. A more thorough review of the35modern implementation of the method of adiabatic switching and its applications can be foundin Skodje & Cary (1988).
Implementing adiabatic switching is fairly straightforward. The non-integrable Hamiltonian for which the eigenfrequencies are desired is written in action-angle variables as a sum of
two terms
H(I, Θ) = H0 (I) + H1 (I, Θ),(2.58)

where H0 is an integrable Hamiltonian for which the eigenfrequencies are known (for instance,from the quantization procedure of §2.4), and H1 includes the non-integrable terms. We introduce into equation (2.58) a time-dependent switching function λ(t):H = H0 + λ(t)H1(2.59)

which satisﬁes
λ(t ≤ 0) = 0andλ(t ≥ T ) = 1(2.60)

for some time T which, in order to ensure that the transition from H0 to H is adiabatic, istaken to be much longer than other timescales associated with the system (for instance, thecharacteristic propagation times of equation (2.42), νk−1 = [∂H0 /∂Ik ]−1 ). In addition the ﬁrst
few derivatives with respect to t of λ are, to ensure the smoothness of the transition, chosen to
be continuous at t = 0 and t = T . In this work, λ is taken to bet1
λ(0 ≤ t ≤ T ) = −
sin
T
2π



2πt
T


.

(2.61)

Johnson (1985) discusses the eﬀects of diﬀerent switching functions.
At t = 0 we begin with an eigenmode. That is, our initial conditions are such that
we satisfy the quantization conditions for the Hamiltonian H0 . We then numerically integrate
Hamilton’s equations of motion (2.38) under the inﬂuence of the Hamiltonian H given in equation(2.59). As t increases so does λ(t), adiabatically switching on the inﬂuence of the non-integrable
term H1 , and allowing the eigenfrequency of the system to slowly adjust from its original value.
At t = T we arrive at the eigenfrequency of the non-integrable Hamiltonian H.

36
Unfortunately, the mathematical justiﬁcation of adiabatic switching is complete only forsystems with one degree of freedom. For these systems, which include the pendulum discussed
above, it can be proven (for example Arnold 1978, Section 52F) that the action variable is astrict adiabatic invariant if the characteristic frequency is never equal to 0. Qualitatively, thisis a reasonable statement. Since the switching is slow, at any instant the system is governedby a one-dimensional, time-independent Hamiltonian. Such systems are always integrable since
they possess one constant of the motion, the Hamiltonian itself, and hence the phase-space
trajectories are restricted to a one-dimensional closed curve (that is, a one-dimensional invariant
torus). Because of this conﬁnement, a trajectory samples all of the available phase space during
one period of motion and the system can, in a sense, average over the entire disturbance andkeep the action variable constant. The restriction to νk = 0 is reasonable when one realizes
that if νk = 0 at any point, the period of the system is inﬁnite. It is then impossible for theswitching to occur on a timescale which is long compared to the system’s motion and a complete
exploration of phase space cannot be made.
For systems with n > 1 degrees of freedom, the situation is less certain. A primaryconcern is the existence of the invariant tori, as time-independent Hamiltonians are always
integrable only when n = 1. As we brieﬂy discussed in §2.3, the KAM theorem concerns theexistence of these tori. A more formal, although still not rigorous, statement of the theorem is:10
If an unperturbed, integrable system is non-degenerate, then for suﬃciently small perturbations,most non-resonant invariant tori do not vanish. These invariant tori form a majority when theperturbation is small.
The degeneracy condition is a generalization of the requirement we encountered in the
one degree of freedom case that νk = 0. It accounts for the possibility that if the characteristicfrequencies are commensurable, a trajectory will not fully explore an invariant torus and hence
10 A precise statement may be found in Arnold (1978), Appendix 8B.

37will not sample the available phase space. As usually stated, it is written asdet∂ 2 νi (I)∂ 2 H0 (I)
= det
= 0.
∂Ij
∂Ii Ij(2.62)

This condition is suﬃcient but not necessary, and the unperturbed system which we consider
in Chapter 3 violates it. For a system with two degrees of freedom, another less stringent
requirement (see Lichtenberg & Lieberman 1983, equation (3.2.12)) states that if the ratio of
the characteristic frequencies is ν2 /ν1 = r/s with r and s integers, thenr22
∂ 2 H0∂ 2 H02 ∂ H0
−
2rs
+
s
= 0
∂I12
∂I1 ∂I2∂I22(2.63)

is suﬃcient to ensure the existence of invariant tori under small perturbations. Our unperturbed
state does satisfy this condition.
Although the KAM theorem assures us of the existence of invariant tori (at least for small
perturbations), during the switching we are essentially guaranteed to encounter regions where
the degeneracy condition is violated and hence no tori exist. For a system with two degrees
of freedom, the two-dimensional tori partition the three-dimensional energy shell into an inside
and an outside, eﬀectively blocking the trajectories on destroyed tori from wandering throughall of phase space. However, if our original time-independent system has two degrees of freedom,
the time-dependent switching Hamiltonian of (2.59) is equivalent to a time-independent systemwith three degrees of freedom. The complete phase space is six-dimensional with the timeindependent Hamiltonian conﬁning the system to a ﬁve-dimensional energy shell. The invariant
tori are three-dimensional and do not partition the ﬁve-dimensional space. Therefore, when we
switch through a region violating the degeneracy condition, the trajectory is free to wander inphase space away from our original family of invariant tori. This process is known as Arnold
diﬀusion and is ubiquitous in non-integrable systems with more then two degrees of freedom.
Arnold diﬀusion is of particular importance because it places an eﬀective limit on the size of T
in equation (2.61). As T → ∞ a switched system remains near destroyed tori for longer periodsof time, undergoes more Arnold diﬀusion, and eventually wanders randomly in phase space. We
discuss our choice for T in §3.2.

38
There is the ﬁnal possibility that the end state to which we switch does not contain any
invariant tori. The KAM theorem merely asserts that most invariant tori exist for small perturbations; as the perturbation strength grows, more and more tori are destroyed until eventuallynone remain. This diﬃculty explains why we were forced to abandon the quantization procedureof §2.4 and adiabatic switching cannot be completely justiﬁed when this occurs.
Although the formal justiﬁcation for adiabatic switching is incomplete for the systems
under consideration, the empirical evidence for its success is strong. It has been used to investigate molecular spectra (Patterson 1985; Smith et al. 1987), the energy of nucleon systems(Brut & Arvieu 1993), and optical cavities (Nöckel & Stone 1997), and in each case performed
well despite violating the strict constraints on its use. One justiﬁcation is the possibility thatswitched systems spend a suﬃciently small amount of their time in regions where the procedure
is formally invalid that the ﬁnal results are only slightly perturbed. Still, the strongest support
comes from its long-standing success, and it is in that spirit that we use it to explore the eﬀects
of convective motions on solar eigenfrequencies.

Chapter 3

Implementation and Tests

Having established the underlying theory, we now turn to a more detailed description of
its application to the determination of convective eﬀects on p-mode eigenfrequencies. In §3.1 weexamine our reference state in some detail. Section 3.2 documents the equations necessary to
implement adiabatic switching as well as the computer program we have written to solve them.
Two tests of the method and our code are discussed in §3.3 and §3.4.

3.1

The Reference StateWe are interested in ﬁnding the eigenfrequencies of non-integrable systems which can be

written in the form of equation (2.45). Our reference state, corresponding to H0 in equation(2.45) and which we will show below to be integrable and non-degenerate, is a two-dimensional,
adiabatically stratiﬁed, plane-parallel polytrope. The horizontal and vertical unit vectors are x̂
and ẑ, respectively. We also deﬁne z = 0 as the level where p = ρ = 0 and take the vertical axis
to increase inwards, opposite the radial axis in a spherical polar coordinate system. Polytropic
systems are those where the pressure p and density ρ are assumed to be related by
p
=
p0



ρ
ρo

1+ µ1,(3.1)

where p0 and ρ0 are constants and µ is the polytropic index, which for an adiabatically stratiﬁed
system may be written in terms of the adiabatic exponent deﬁned in equation (2.5) as µ =1 + 1/Γ1 . For models in radiative equilibrium µ = 3, although observations suggest µ = 3.5

40
is a better ﬁt for the solar photosphere. We adopt the former value in this work. Enforcing
hydrostatic equilibrium leads to the relations
p
=
p0



z
z0

µ+1

ρ
=
ρ0

and



z
z0

µ
,

(3.2)

where z0 is a constant. Finally, the adiabatic sound speed and cutoﬀ frequency can be computedfrom equations (2.4) and (2.28) respectively, and are given by
c2 =

gz
µ

and

ωc2 =

µ(µ + 2) 2
c ,
4z 2

(3.3)

where g is the gravitational acceleration. Assuming the ideal gas law holds, equations (3.2) and(3.3) imply that the square of the sound speed is proportional to the temperature.
For this system (compare with equation (2.35)), the Hamiltonian isω(x, z, kx , kz ) =gz 2g(µ + 2)
(kx + kz2 ) +
,
µ
4z

(3.4)

where kx and kz are the wavenumbers in the horizontal and vertical directions. As this twodegree of freedom Hamiltonian is independent of time as well as the coordinate x, there are twoconstants of the motion: ω and kx . These constraints are enough to ensure that this Hamiltonianis integrable. A raypath for this system is shown in Figure 3.1. Although this geometry is planeparallel rather than spherically symmetric, the raypaths are obviously quite similar to those of
Figure 2.1. In this ﬁgure the upper turning point at z ≈ 0.3, an eﬀect produced by the acoustic
cutoﬀ frequency, is clearly visible.
We turn next to the formulation of the reference system in action-angle coordinates.
Combining the expressions for c2 and ωc2 from equation (3.3) with the dispersion relation ofequation (3.4), we may solve for kz :kz2 =

ω 2 − ωc2
µω 2
µ2
− kx2 =
− 2
2
c
gz
4z



2
1+
− kx2 .
µ

(3.5)

In the extended periodic domain described in §2.4.2 we can integrate along a curve at constantx and write equation (2.43) as
Iz =

1
2π


kz dz =

1π

 z2
kz dz,
z1

(3.6)

41

Raypaths in a Polytropez (Scaled Units)

0.0

0.5

1.0

1.5
0

100

200
300
x (Scaled Units)

400

Expanded Viewz (Scaled Units)

0.0

0.5

1.0

1.5
180

200

220

240

x (Scaled Units)
Figure 3.1: A sample raypath in an adiabatically stratiﬁed, two-dimensional, polytropic envelopewith µ = 3. This ray corresponds to a mode with n = 1,  = 200. The top panel shows theentire horizontally periodic cavity with length 2πR . The large value of  produces the many
oscillations in the horizontal direction. The bottom panel shows an expanded version of the top
panel. Note the two caustic surfaces at the vertical extremes of the cavity.

42
where z1 and z2 are the upper and lower turning points of the ray. With kz given by equation(3.5), and since kx is a constant of the motion, the integral in equation (3.6) can be analyticallycomputed. The result is
Iz =

1
2



2
µω 2
−µ 1+
gkx
µ


.

(3.7)

A similar integral for the horizontal action variable gives
Ix = kx R ,(3.8)

and thus the reference Hamiltonian in action-angle coordinates is



1 
 2g
2 2
µ

H =ω=
Ix Iz +
1+
.
µR
2
µ

(3.9)

The characteristic frequencies, derived from equation (2.42), are
νx =

∂ω
ω
=
∂Ix
2Ix

and

νz =

∂ω
=
∂Iz

ω


1 .
µ
2 2
2 Iz + 2 1 + µ(3.10)

Although this Hamiltonian violates the stringent degeneracy criterion of equation (2.62), it
satisﬁes the looser standard of equation (2.63). Hence the KAM theorem applies to this reference
Hamiltonian.

3.1.1

Reference State Results
Several properties of this polytropic reference state will be of use in later sections. Per-

haps most importantly Christensen-Dalsgaard (1980), beginning with the perturbed structure
equations (2.7)–(2.10), solved for the eigenfrequencies of the system. Their exact values areω2 =2gkx 
µ
n+
µ
2

for

n = 1, 2, ...

(3.11)

where kx = /R and  is an integer which we associate with the spherical harmonic degree.1
The upper and lower turning points of a ray occur when kz = 0. Using equation (3.4) andthe condition given by equation (3.11), the cavity for an eigenmode is bounded in the vertical
1 In the spherically symmetric case, kx =p( + 1)/R .
43
direction byµ  R
z1,2 = n +
1∓
2




µ(µ + 2)1−.
(µ + 2n)2(3.12)The midpoint of a modal cavity is merelyµ  R
,zmid = n +2(3.13)

which demonstrates, as noted earlier, that the depth of the cavity scales roughly as n/.
The evolution of a ray is given by Hamilton’s equations, the general form of which are
given in equation (2.38). Specializing to a polytrope, we have
dx
c2 kx
=
dt
ω
dz
c2 kz
=
dt
ω
dkx
=0
dt
dkz
ω
=− ,
dt
2z

(3.14)(3.15)(3.16)(3.17)

which explicitly contains the invariability of kx .
We can also determine the characteristic periods of motion in each coordinate from equations (3.10) and (3.11). Doing so, we have
Tz =

2π
4π 
µ
=
n+
νz
ω
2

(3.18)

2π
4π
=
.
νx
ω

(3.19)andTx =Finally, we arrive at the ratio of the horizontal and vertical periods:
Tx

=
.
Tz
n + µ2(3.20)For the high-degree modes under consideration the horizontal period is substantially longer thanthe vertical. In Figure 3.1, the ratio is 80 : 1.

44
3.1.2

Quantization Conditions
Since the reference state is integrable we can apply to it the quantization procedure of§2.4. As expected, in the appropriate limit it reproduces equation (3.11). Starting with equation(3.6), and noting that the integration path touches two caustics, equation (2.47) becomes
1
Iz =
2π



1
kz dz =
π

 z2
z1



1
kz dz = n −2

forn = 1, 2, . . .

(3.21)

with a ﬁnal result of


1 
2gkx
2 2
1 µ
ω =
n− +
1+
µ
2 2
µ
2

for

n = 1, 2, . . .

(3.22)

If we also evaluate equation (2.47) on a curve of constant z we encounter no caustics and ﬁnd
that
Ix =

1
2π


kx dx =

12π

 2πR
kx dx = for

 = 0, 1, . . .

(3.23)

0for which
kx =


R

for

 = 0, 1, . . .

(3.24)

Although equation (3.24) is the correct horizontal quantization condition, the eigenfrequencies
of equation (3.22) are not the same as those of equation (3.11). Recall, however, that thederivation of the quantization conditions relies on formulating the solutions as local plane waves,an approximation which is only valid when the wavelength of the variations in the reference stateis larger than the wavelength of the solution. The limit in which we expect the two results toagree, then, is that where the eigenfunction rapidly oscillates in the vertical direction; in otherwords, when n

1. In this limit the results agree. It should be noted that even with n = 1 andµ = 3 the relative error in the eigenfrequencies is ≈ 3%, with agreement improving as either µor n increases.
In this case we were able to evaluate the action variables due to the simple form of
the sound speed and acoustic cutoﬀ frequency given in equation (3.3). However, even when
the action variables cannot be derived analytically, the quantization procedure will still work,

45
provided the system is integrable. For solar p modes the most realistic condition guaranteeing
integrability is for the sound speed and cutoﬀ frequency to be independent of x, implying that kx
is a constant of the motion. The quantized action variables can then be evaluated numericallyto determine the eigenfrequency ω.

3.2

Adiabatic Switching
Our work in Chapter 2 demonstrated that the dispersion relation of equation (3.4) gov-

erns ray motion. However, this derivation was rather restrictive in that it ignored interesting
eﬀects such as advective motions, magnetic ﬁelds, and rotation, all of which contribute terms
to the dispersion relation and can aﬀect ray trajectories. Of the three, only advective motions
are necessary for describing convective structures, and so in this initial examination we ignore
magnetism and rotation. In a frame of reference moving with the ray, advective motions willcontribute a Doppler shift to the ray frequency implying that the correct dispersion relation is:(ω − k · v)2 = k 2 c2 + ωc2(3.25)

where v is the velocity ﬁeld. We adopt equation (3.25) as the dispersion relation for a ray propagating in a convective atmosphere in the remainder of this work. Note that for the polytropic
reference state of §3.1 equation (3.25) reduces to equation (3.4) if v = 0.
After some algebraic manipulation we can write this Hamiltonian in the form of equation(2.59):



1

1
µ(µ + 2) 2µ(µ + 2) 222
ω = c0 k ++ λ(t) (c − c0 ) k ++k·v ,
4z 2
4z 2
where c0 =(3.26)

gz/µ is the polytropic sound speed, c and v are the sound speed and velocity

ﬁeld of the convective model, and the cutoﬀ frequency of equation (3.3) replaces ωc . The true
deﬁnition of the cutoﬀ frequency, given in equation (2.28), depends on the density scale heightHρ which, in a self-consistent model, will almost certainly undergo perturbations somewhat
diﬀerent from those aﬀecting the sound speed. However, this is a minor eﬀect and we assume in46
this work that the only perturbations to the cutoﬀ frequency arise from variations in the soundspeed in the manner shown.
For this more complicated dispersion relation, Hamilton’s equations are
 2 


dx
c20 kxc0 kx c=+ λ(t)− 1 + vx
(3.27)
dt
Ω
Ω
c0
 2 


c2 kz
c kz cdz= 0 + λ(t) 0(3.28)− 1 + vz
dt
Ω
Ω
c0


dkx
Ω ∂c
∂vx
∂vz
= −λ(t)+ kx+ kz(3.29)
dt
c0 ∂x
∂x
∂x
 


dkz
Ω
1 ∂c
1
∂vx
∂vz
= − − λ(t) Ω
−
+ kx
+ kz
,
(3.30)
dt2z
c0 ∂z 2z∂z∂zwhere Ω =c20 (kx2 + kz2 ) + ωc2 is the instantaneous polytropic frequency. Note that in thepolytropic limit where c = c0 and v = 0, we return to equations (3.14)–(3.17).
Although we include the function λ(t) in our equations, we have not yet addressed our
choice for the switching time T of equation (2.61). In order to ensure adiabaticity, we mustchoose T to be much longer than the longest dynamical timescale associated with the system.
Heuristically, this ensures that the ray has a chance to explore the available phase space during
the switching process. The relevant timescales in this case are the time for the ray to undergo
an oscillation in the vertical and horizontal directions — taking our cue from equation (3.20),
we choose T

Tx . As we discussed in §2.5, the presence of Arnold diﬀusion in non-integrable

systems sets an upper constraint on the choice of T . Our numerical experiments show thatchoosing T ≈ 40Tx oﬀers a satisfactory compromise between the two constraints. This choiceagrees with that made in other implementations of adiabatic switching.
Beginning from diﬀerent initial conditions should, in principle, lead to identical eigenfrequencies for the ﬁnal state. For suﬃciently complex systems our numerical experiments ﬁnd ascatter in the results, presumably arising from the combined eﬀects of a ﬁnite switching timeand Arnold diﬀusion. This scatter is observed by other investigators and it has been foundempirically that better results are obtained by implementing the procedure for several randomly
distributed initial conditions and averaging the ﬁnal eigenvalues (although see the mathematicaljustiﬁcation presented in Skodje & Cary (1988)). This averaging is particularly necessary for the47systems of §4.3 and §4.4 where we typically ﬁnd the eigenfrequencies for 20 initial conditions, achoice which we ﬁnd suitably constrains the result.
With these considerations, we have written a computer program to implement adiabatic
switching by numerically integrating equations (3.27)–(3.30). As our reference state we choose
an adiabatically stratiﬁed, plane-parallel polytrope with an index µ corresponding to the average
stratiﬁcation of the photospheric model. We approximate the photospheric model itself as a gridof sound speeds c and velocity vectors v = (vx , vz ). Although our current code requires equally
spaced grids, they are not a necessity. Indeed, since the upper portion of the solar convection
zone contains many scale heights of material an adaptive grid spacing would be computationallybeneﬁcial. The spacing of our grid elements represents, as usual, a trade-oﬀ between accuracy
and computational speed. Currently we only work with two-dimensional, plane-parallel convective models with thermal and velocity structures, though there is no fundamental barrier toextending the method to explore other structures including, for example, those with magnetic
ﬁelds or three dimensions.
To integrate equations (3.27)–(3.30), the code uses a slightly modiﬁed version of the DE
package of Shampine & Gordon (1975): a variable-order, variable-timestep Adams-BashforthMoulton PECE (predict-evaluate-correct-evaluate) method. The modiﬁcations include an updating of the code to FORTRAN 90 and an adaptation for horizontally periodic boundary
conditions. Fourth-order central diﬀerences are used to evaluate the partial derivatives of thesound speed and velocity ﬁelds everywhere within the cavity, except near the top and bottom
boundaries where fourth-order forward and backward diﬀerences are used. The DE package
allows the user to specify tolerances for both the relative and absolute errors. We found thatthe integration was suﬃciently quick and accurate if both were set to 10−9 . Integration times
for a trajectory greatly vary depending on the parameter range and the convective model. For
a n = 1,  = 200 mode propagating in the perturbation of §3.3, integration of one trajectorytakes ≈ 5 minutes to complete on a 194 MHz SGI Power Onyx. A further description of thecode’s history and structure is presented in Appendix D.

48
Finally, we have scaled our variables in order to simplify the calculations. Unless otherwise
noted, all lengths are quoted in units of 109 centimeters and all times in units of 102.5 seconds.
So, the solar surface gravity is g = 2.7397, the Sun’s radius is R = 69.599 and a canonical5-minute oscillation has a frequency of ω/(2π) = 1.0541.

3.3

A Test of an Integrable System
To test the methodology and our code, we examine an integrable system whose eigen-frequencies can be computed both with adiabatic switching and through EBK quantization.
Consider the simple, but not particularly realistic, case of a photosphere with no ﬂuid motions
and a sound speed closely related to the polytropic value of equation (3.3):

c=

gz
µ

 12
[1 +  sin(Kz)] ,(3.31)

where  is the amplitude and K the wavenumber of the perturbation. Horizontal invarianceimplies that kx is a constant of the motion and hence that the system is integrable at every
point during the switching. For this perturbation, the acoustic cutoﬀ frequency of equation (3.3)
becomesωc =g(µ + 2)4z

 12
[1 +  sin(Kz)] .

(3.32)As mentioned above, the true expression for ωc , equation (2.28), directly depends on the density
scale height Hρ as well as the sound speed. Equation (3.32) is not fully self-consistent, but we
ignore the resulting discrepancies.
In order to derive the eigenfrequencies from EBK quantization, we must solve the integral
equation z2 

 z2
kz dz =
z1

z1

µω 2
µ2
− 2
2
gz[1 +  sin(Kz)]
4z



 12
2
π
1+
− kx2
dz, =
µ
2

(3.33)

for the eigenfrequency ω where z1 and z2 are the depths at which the integrand vanishes.
This equation is only valid for the lowest (n = 1) mode; for higher modes the right-hand sideincreases by the necessary factors of π. Perhaps the most straightforward method of evaluating

49
this equation, and the one we use to derive the results presented below, is a bisection-liketrial-and-error method.
We compute relative frequency shifts ∆ω/ω arising from the sound speed perturbationof equation (3.31) via adiabatic switching and EBK quantization and list the results in Table
3.1. The unperturbed frequencies, which we list for three of the modes in the ﬁrst rows, are
not identical for the two methods. For adiabatic switching we use the true eigenfrequencies ofequation (3.11) while for EBK quantization we use the results of equation (3.22). The fractionalfrequency shifts do agree, however. The results suggest three separate trends in the frequency
shifts:(a) The fractional frequency shift increases as n increases.
(b) The fractional frequency shift decreases as  increases.
(c) The fractional frequency shift increases as  increases.
In order to explain these results we turn to an analysis based on the work of Gough (1993),
wherein we expand the EBK quantization condition of equation (3.6) around the polytropic
reference state.
From equation (3.31),
c = c0 + ∆c =

gz
µ



 12
+

gz
µ

 12
sin(Kz)(3.34)

where c0 is the polytropic sound speed. Making similar expansions to ω and ωc , the EBK
quantization condition for the perturbed state is z2 

 z2
kz dz =
z1

z1

 12
(ω0 + ∆ω)2 − (ωc,0 + ∆ωc )2π
2− kx dz = ,(c0 + ∆c)22

(3.35)

where quantities with a 0 subscript correspond to the polytropic reference state. Next, assumingthat ∆ω, ∆c, and ∆ωc are small parameters, we expand and discard higher-order terms. Afterrearranging, we have


 z2
kz,0
z1

2
1+ 2 2c0 kz,0


 12
π
2 ∆c
2 ∆c
ω0 ∆ω − ω0+ ωc
− ωc ∆ωcdz = ,c0c02(3.36)

50

Table 3.1: Frequency shifts for a vertical sound speed perturbation. Some entries have been
repeated in order to display trends.

na

b

1f
2
1
1
2
3
4
1
1
1
1
1
1
1
1
1

200
200
400
200
200
200
200
100
200
400
800
200
200
200
200
200

c

0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.01
0.02
0.04
0.08
0.16

Kc

[∆ω/ω]AS d[∆ω/ω]EBK e

0.05
0.05
0.05
0.05
0.05
0.05
0.05
0.05
0.05
0.05
0.05
0.05
0.05

3.622
4.286
5.123
1.7 × 10−3
2.4 × 10−3
3.1 × 10−3
3.8 × 10−3
3.4 × 10−3
1.7 × 10−3
8.7 × 10−4
4.4 × 10−4
4.1 × 10−4
8.4 × 10−4
1.7 × 10−3
3.4 × 10−3
6.8 × 10−3

3.576
4.247
5.057
1.7 × 10−3
2.4 × 10−3
3.0 × 10−3
3.7 × 10−3
3.4 × 10−3
1.7 × 10−3
8.7 × 10−4
4.4 × 10−4
4.2 × 10−4
8.5 × 10−4
1.7 × 10−3
3.4 × 10−3
6.8 × 10−3

a Vertical mode number.
b Horizontal mode number.
c Perturbation parameter from equation (3.31).
d Relative frequency shifts from adiabatic switching.
e Relative frequency shifts from EBK quantization.
f The ﬁrst three rows give unperturbed frequencies for comparison purposes.

51
where kz,0 is the positive root of equation (3.5) and is thus a function of only unperturbed
quantities. We expand again, ignore any perturbations to the cavity boundary (which we will see
below are small), and cancel the zeroth-order term leaving, after some algebraic simpliﬁcations,
∆ω
=
ω0

 z2
z1

∆c dz
c0 kz,0 c20

  z2
z1

dz
kz,0 c20

−1.

(3.37)

Note that we can derive the same result by expanding the horizontal quantization condition
as well. In this casekz dz is simpler to evaluate thankx dx since the former includes theconstant kx2 within the integrand while the latter would include terms proportional to ∆kz .
Thus, the fractional frequency shift is the weighted average of the fractional sound speed
perturbation over the cavity. A mathematically satisfactory explanation for the trends noted
above could come from integrating equation (3.37) analytically for the perturbation given inequation (3.34). However that approach is only possible for simple perturbations and so, in
preparation for the results of Chapter 4, we cast our arguments in a more physically motivatedform.
First, ∆c/c0 =  sin(Kz) ≈ Kz when, as is the case for the parameters we consider,Kz is small. Since the perturbation is positive and depends linearly on , it is not surprisingthat the fractional frequency shifts seen in Table 3.1 are always positive and scale linearly withthe perturbation strength. The other two eﬀects are somewhat more subtle as they arise fromchanges in the depth of the acoustic cavity. From equation (3.13), we see that the location of acavity varies inversely with  and linearly with n + µ/2. As n and  change, a ray encounters
perturbations of diﬀerent strengths which impart diﬀerent frequency shifts. In this case theperturbation strength depends linearly on depth and we expect the fractional frequency shift
to decrease inversely with  and increase as n + µ/2. These trends are reproduced by adiabatic
switching.
We conclude that for vertically dependent sound speed perturbations, the fractional frequency shift depends linearly on the strength of the perturbation. The variation with the modeparameters depends on the functional form of the perturbation. For ∆c/c0 ∝ z p , we expect the52Raypath Comparisonz (Scaled Units)

0.0

0.5

1.0

1.5
290

300

310
320
330
x (Scaled Units)

340

350

Frequency During Shift
3.650

ω (Scaled Units)

3.645
3.640
3.635
3.630
3.625
3.620
0.0

0.2

0.4
0.6
Switching Time (t/T)

0.8

1.0

Figure 3.2: A raypath and the evolution of he eigenfrequency for the perturbation of §3.3.
Results for a n = 1,  = 200 mode and the sound-speed perturbation of equation (3.31) with
 = 0.16, K = 0.05 are plotted. See the text for further discussion.

53
fractional frequency shifts to vary as (n + µ/2)p and −p .
Finally, we present the perturbed raypaths and the modiﬁcation of the frequency during
adiabatic switching in Figure 3.2. The top panel shows two raypaths in a portion of the acoustic
cavity. The dotted line is a raypath in the polytropic reference state while the solid line is a
raypath in the perturbed model whose structure is given by equation (3.31). The cavity in theperturbed model is shifted slightly toward the surface due to the increased sound speed. The
second panel shows the evolution of the ray frequency as a function of switching time. The
smooth transition occurs when switching between integrable systems and, as we will see, does
not occur for more complicated perturbations.

3.4

Doppler Shifts
As a ﬁnal test before turning to non-integrable systems, we consider the eﬀect of aconstant horizontal ﬂow on the eigenfrequencies. From equation (3.29) we see that a suﬃcientcondition for kx to be constant (thus maintaining integrability) is ∂c/∂x = ∂vx /∂x = ∂vz /∂x =0. In particular this occurs for the stringent criteria that vx is a non-zero constant, vz = 0,
and c = c0 . Other selections also maintain integrability, but this choice of parameters leads toresults which are easily interpreted.
Table 3.2 shows the frequency shifts from EBK quantization as well as those from adiabatic switching. We see that:(a) The fractional frequency shift decreases as n increases.
(b) The fractional frequency shift increases as  increases.
(c) The fractional frequency shift increases as vx increases.
Notice that the shifts for  and n are of the opposite sense than for the sound speed perturbationof §3.3.
To explain these trends we can expand the EBK quantization condition as we did in theprevious section. In this case, however, the simplicity of the perturbation (a constant vx acting

54

Table 3.2: Frequency shifts for a horizontal velocity perturbation. Some entries have been
repeated in order to display trends.

na

b

1f
2
1
1
2
3
4
1
1
1
1
1
1
1
1
1

200
200
400
200
200
200
200
100
200
400
800
200
200
200
200
200

vx c

[∆ω/ω]AS d[∆ω/ω]EBK e

0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.005
0.01
0.02
0.04
0.08

3.622
4.286
5.123
7.9 × 10−3
6.7 × 10−3
5.9 × 10−3
5.3 × 10−3
5.6 × 10−3
7.9 × 10−3
1.1 × 10−2
1.6 × 10−2
4.0 × 10−3
7.9 × 10−3
1.6 × 10−2
3.2 × 10−2
6.3 × 10−2

3.576
4.247
5.057
7.9 × 10−3
6.7 × 10−3
5.9 × 10−3
5.3 × 10−3
5.6 × 10−3
7.9 × 10−3
1.1 × 10−2
1.6 × 10−2
4.0 × 10−3
7.9 × 10−3
1.6 × 10−2
3.2 × 10−2
6.3 × 10−2

a Vertical mode number.
b Horizontal mode number.
c Horizontal velocity.
d Relative frequency shifts from adiabatic switching.
e Relative frequency shifts from EBK quantization.
f The ﬁrst three rows give unperturbed frequencies for comparison purposes.

55
on a constant kx ) allows us to write down the correct result. The fractional frequency shift isgive by a Doppler shift:

− 1
kx vx
µ 2
∆ω
2gkx 
== kx vx
n+
,
ω
ω
µ
2

(3.38)

where we have used equation (3.11) to substitute for ω. Since kx ∝ , the dependencies noted
above are correctly explained by equation (3.38).
In Figure 3.3 the top panel shows two raypaths in a small portion of the acoustic cavity.
The dotted line is a raypath in a polytropic atmosphere while the solid line is a raypath in themodel including the horizontal ﬂow. The eﬀect of the ﬂow in advancing the raypath motion
is clearly visible. The second panel shows the evolution of the ray frequency as a function ofswitching time. Again, the shift between integrable systems is quite smooth.

56

Raypath Comparisonz (Scaled Units)

0.0

0.5

1.0

1.5
290

300

310
320
330
x (Scaled Units)

340

350

Frequency During Shift
3.90

ω (Scaled Units)

3.85
3.80
3.75
3.70
3.65
3.60
0.0

0.2

0.4
0.6
Switching Time (t/T)

0.8

1.0

Figure 3.3: A raypath and the evolution of the eigenfrequency for the perturbation of §3.4.
Results for a n = 1,  = 200 mode and a Doppler shift perturbation with vx = 0.08. See thetext for further discussion.

Chapter 4

Convective Structures

The sound speed and velocity ﬁelds of the solar convection zone do not form an integrable
system. Within this chapter we introduce several models of convective structures containing nonintegrable perturbations and examine their eﬀects on p-mode eigenfrequencies. Although themodel of §4.1 contains only velocity perturbations, its simple structure provides insight into the
behavior of more complex convective structures. In §4.2 we examine a two-dimensional RayleighBénard convective model which contains both sound speed and velocity perturbations. Morerealistic convective simulations appear in §4.3 where we treat a model of a cool, compressible
plume and §4.4 where we consider a two-dimensional slice of a three-dimensional simulation ofcompressible convection.

4.1

Cell-like Velocity Perturbation
As a ﬁrst attempt we choose to model convective cells as velocity perturbations super-

imposed on a polytropic temperature stratiﬁcation. Our velocity ﬁeld is derived from a stream
functionψ =  sin

2πx
L




sin

2πz
L


,

(4.1)

where L = 2πR is the horizontal length of the cavity and  is the perturbation strength.
Although this parameterization unphysically assumes an incompressible ﬂuid, we are more concerned with the physical insight we can derive from this model than its self-consistency. With

58
the normal stream function relations (v = ∇ × ψŷ where ŷ is a unit vector perpendicular toour other coordinates), equation (4.1) gives the velocity ﬁeld

sin
vx =
R



x
R




cos

z
R


and


cos
vz = −
R



x
R




sin

z
R


.

(4.2)

We plot streamlines for this ﬂow in Figure 4.1.
Given this perturbation, we use adiabatic switching to investigate the dependence of
eigenfrequencies on the parameters , n, and . Our results are shown in Table 4.1 and suggest
the following trends:(a) The fractional frequency shift decreases as n increases.
(b) The fractional frequency shift increases as  increases.
(c) The fractional frequency shift increases as  increases.
Although the sense of these trends are the same as for the Doppler shift of §3.4, the functional
form of the shifts is diﬀerent. Notably, the frequency shift does not linearly depend on ; instead,
it is quadratic in .
To explain these relations, we turn, with an important caveat, to the mathematicalframework developed in §3.3. We are no longer working with an integrable system and hence
the quantization conditions of equation (2.47) are not well deﬁned. For example, consider thevertical quantization condition for the lowest mode: kz dz =

 12
(ω − k · v)2 − ωc22− kx dz = π.
c2

(4.3)

In earlier applications the integrand depended solely on z and the integral was performed at anarbitrary value of x. Now kx depends on both x and z and is no longer a constant of the motion;
hence, the integral depends on the x at which it is evaluated. The condition is ill-deﬁned, but
the KAM theorem and the empirical success of adiabatic switching suggest that the quantization
conditions are accurate for small perturbations. In this spirit, we average equation (4.3) overthe horizontal dimension while acknowledging that we should use the results for guidance, andnot absolute statements, concerning the eigenfrequency shifts.

59

Streamfunction for velocity perturbation
0z (Scaled Units)

1

2

3

4

5
0

100

200
300
x (Scaled Units)

400

Figure 4.1: Streamlines for the velocity perturbation of §4.1. The two dashed lines near the topof the plot delineate the resonant cavity of the n = 1,  = 200 polytropic mode shown in Figure
3.1. The WKB criterion is clearly met in both directions. In order to resolve the cavity, theclosure of the streamlines at lower depths is not shown.

60

Table 4.1: Frequency shifts for convective cells composed of velocity perturbations. Some entries
have been repeated in order to display trends.

na

b

c

[∆ω/ω]d

1
2
3
4
5
1
1
1
1
1
1
1
1
1
1
1

200
200
200
200
200
100
200
300
400
500
800
200
200
200
200
200

0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.01
0.02
0.04
0.05
0.08

−1.3 × 10−3
−9.0 × 10−4
−7.0 × 10−4
−5.7 × 10−4
−4.8 × 10−4
−6.3 × 10−4
−1.3 × 10−3
−1.9 × 10−3
−2.5 × 10−3
−3.2 × 10−3
−5.1 × 10−3
−7.9 × 10−5
−3.2 × 10−4
−1.3 × 10−3
−2.0 × 10−3
−5.1 × 10−3

a Vertical mode number.
b Horizontal mode number.
c Perturbation strength.
d Relative frequency shifts from adiabatic switching.

61
Our system contains no perturbation to the sound speed, and thus none to the acoustic
cutoﬀ frequency, but now we must include both the Doppler term of equation (3.25) as well asa perturbation to kx . So, kz dz =

2
(ω0 + ∆ω − k0 · v)2 − ωc,0
− (kx,0 + ∆kx )2
2
c0

 12
dz = π,(4.4)

where quantities with a 0 subscript correspond to the polytropic reference state and are thus
independent of x. Averaging in the horizontal direction, treating k0 · v as a small parameter,and performing the same expansion as in §3.3, yields
1
∆ω
=
ω0
L




dx

k0 · v dz+
ω0 kz,0 c20

kx,0∆kx dz
kz,0



1
L





dz
kz,0 c20

dx

−1
,

(4.5)

where L is the horizontal extent of the domain. From equations (3.14)–(3.17), Hamilton’s
equations for the unperturbed state, we can rewrite the ﬁnal term in the ﬁrst set of brackets as
kx,0∆kx dz =kz,0∆kx dx.

(4.6)

From the quantization conditions this term is identically zero, leaving as a ﬁnal result
∆ω
=
ω0




dx

k0 · v dz
ω0 kz,0 c20

 


dx

dz
kz,0 c20

−1
(4.7)

which should be compared to equations (3.37) and (3.38). Again, we ﬁnd that the frequency
shift depends on the weighted average of the perturbation over the cavity. As an aside, we notethat a similar treatment beginning with the horizontal quantization condition yields
∆ω
=
ω0




dz

k0 · v dx
ω0 kx,0 c20

 


dz

dx
kx,0 c20

−1
.

(4.8)

As we now show, these quantization conditions predict that the velocity perturbation ofequation (4.2) produces no ﬁrst-order frequency shift. The sole x-dependence in the integrandsof equation (4.7) arises within v. This dependence is sinusoidal for the perturbations we areconsidering and hence the horizontal integral of equation (4.7) is identically zero. Somewhat
more generally, we see that this result holds for any convective model where, at every depth,
the velocity components satisfy
 L

 L
vx dx =

0vz dx = 0.
0

(4.9)

62
In such a ﬂow, a ray whose frequency is boosted at a given point suﬀers a frequency retardation
at a complementary point elsewhere in the domain. To ﬁrst order, there is no net eﬀect.
Only a small subset of ﬂows exhibit the symmetry necessary for equation (4.9) to be true.
However, equation (4.7) allows us to make a much more general statement about the eﬀects ofﬂows on eigenfrequencies. The shift produced by the vertical component of a ﬂow is
∆ω
=
ω0




dx

vz 1
dz
ω0 c20

 


dx

dz
kz,0 c20

−1,(4.10)

where the dependence on kz,0 has been cancelled in the ﬁrst term. This cancellation is vital as
the path integral in the numerator is now


vz 1
dz =
ω0 c20

 z2
z1

vz 1
dz +
ω0 c20

 z1
z2

vz 1
dz = 0,ω0 c20(4.11)

where z1 and z2 are the vertical boundaries of the acoustic cavity. Thus, the vertical componentof any advective ﬂow produces no net ﬁrst-order eﬀect. This is not the case for horizontal ﬂows
which, as we saw in §3.4, can produce ﬁrst-order shifts.
The asymmetry arises from the observation that rays are trapped in the vertical direction
but not in the horizontal.1 Since a ray never reverses its horizontal motion, the horizontal
wavenumber kx of any ray maintains the same sign along its trajectory. Conversely, the verticalwavenumber kz changes sign at each caustic. Because of this asymmetry, the expressions
and





dx

dz have very diﬀerent mathematical meanings:


 L
dx =


dx = Lwhiledz = 0.

(4.12)

0

For a more physical explanation, consider that during an ideal switching (of inﬁnite extent intime) a ray comes arbitrarily close to any point in the cavity. The horizontal velocity ﬁeld at agiven point interacts with rays with only one sign of kx and imparts a net frequency shift. At
the same point the vertical ﬂow will encounter two manifestations of every ray, diﬀering only inthe sign, but not the magnitude, of kz . To ﬁrst order the eﬀect of the vertical ﬂow on these two
1 It is possible to construct ﬂows where a ray is trapped in the horizontal direction as well, however they do
not correspond with reasonable solar ﬂows and we do not consider them here.

63
rays will be equal and opposite, leaving no net frequency shift.2
For the ﬂow under consideration, the horizontal contribution to the frequency shift iszero, to ﬁrst order, from the argument summarized by equation (4.9). The eﬀects of verticalﬂows are negated, to ﬁrst-order, by either equation (4.9) or (4.11). Regardless of the causal
agent, we expect a non-zero second-order frequency shift, as can be shown by the following
argument, originally due to Brown (1984).
A region containing anti-symmetric ﬂows produces no ﬁrst-order frequency shift. However, a ray will take more time to traverse a region when traveling against the ﬂow (k · v < 0)than to cross a region where the ray and the ﬂow are co-directional (k · v > 0). The diﬀerencein crossing times between the two regions is of order k · v/ω, while the frequency shift in the
counter-propagating region is −k · v/ω. The result is a fractional frequency shift
∆ω
∝−
ω



k0 · vω

2.

(4.13)

With this result we can anticipate the frequency shifts from a domain including the ﬂow
of equation (4.2). First, we expect to observe a fractional frequency shift proportional to 2 .
From equations (3.24) and (3.11) we note that k is proportional to  while ω 2 is proportionalto  and to (n + µ/2). Hence, we also expect a fractional frequency shift which is proportionalto , and inversely dependent on n. As we see in Table 4.1, adiabatic switching conﬁrms these
expectations.
Although we have worked solely with advective ﬂows in this section, our results can beextended to sound speed perturbations. The shift from a perturbation with both horizontal andvertical dependencies can be written as
∆ω
=
ω0




dx

∆c dz
c0 kz,0 c20

 


dx

dz
kz,0 c20

−1,(4.14)

which reduces to equation (3.37) in the appropriate limit. Thus, any antisymmetric perturbation
2 In earlier calculations involving integrals over z, such as equation (3.6), the positive and negative branchesof kz prevent the cancellation of the integral.

64
to the sound speed satisfying
∆c dx = 0

(4.15)

at every depth produces no ﬁrst-order frequency shift. We saw in §3.3 that ﬁrst-order shifts fromsound speed perturbations are not impossible. However, if symmetry dictates that the frequency
shift vanishes to ﬁrst order, a second-order downward shift will arise for a reason analogous tothat described above: rays propagate faster in, and hence spend less time traversing, regionswith higher sound speeds.
In the top panel of Figure 4.2 we show a sample raypath after the system has undergone
adiabatic switching. The large perturbation to the acoustic cavity occurs despite the smallchange in the eigenfrequency (compare this to the behavior in Figure 3.2 where the opposite istrue). Our EBK analysis of the frequency shifts does assume that perturbations to the acoustic
cavity have negligible eﬀects. In this case, the horizontal average of the perturbation to thecavity is zero and we ignore it in the calculations leading to equation (4.7). Presumably there
are higher-order eﬀects which we do not consider.
The bottom panel of Figure 4.2 shows the evolution of the ray frequency as a function ofswitching time. The large oscillations are a characteristic of adiabatic switching in non-integrable
systems.

4.2

Rayleigh-Bénard Boussinesq Convection
While the cell studied in the previous section was helpful in shaping our intuition, it

was unrealistic in that it did not include perturbations to the sound speed. In this section we
explore a model which includes horizontal and vertical variations in both the sound speed andthe velocities. Our reference state remains the adiabatically stratiﬁed, plane-parallel polytrope
described above. We take our parameterization for the convective state from Shirer (1987). It
is a two-dimensional, nonlinear model of steady Rayleigh-Bénard convection derived from theBoussinesq ﬂuid equations. We consider the solution with the lowest mode numbers, two cells

65

Velocity Cell Raypathz (Scaled Units)

0.0

0.5

1.0

1.5

2.0
0

100

200
300
x (Scaled Units)

400

Frequency During Shift
3.624

ω (Scaled Units)

3.622
3.620
3.618
3.616
3.614
0.0

0.2

0.4
0.6
Switching Time (t/T)

0.8

1.0

Figure 4.2: A raypath and the evolution of the eigenfrequency for the perturbation of §4.1.
Results for a n = 1,  = 200 mode and the velocity perturbation of equation (4.2) with  = 0.05.
See the text for further discussion.

66
in the horizontal direction and one in the vertical.
After converting Shirer’s parameterization to our notation, the sound speed isc2 =gz+ α1 cos
µ



2πx
L




 

 
z
z
− α2 sin 2π 1 −,sin π 1 −zT
zT

(4.16)

where L = 2πR is the length of the domain, zT is the vertical extent of the cell, and α1 andα2 are functions of the product of the thermal diﬀusivity and the kinematic viscosity, κν. The
velocity perturbations are given byvx = γ1 sin

2πx
L



 

z
cos π 1 −
zT


and

vz = γ2 cos

2πx
L



 

z
sin π 1 −,zT
(4.17)

where γ1 and γ2 are functions of κ and the product κν. The nonlinear nature of the convection
arises from the ﬁnal term in equation (4.16) which is independent of x. For a linear cell, theperturbation would oscillate in x. Figure 4.3 shows streamlines for the velocity proﬁle given byequation (4.17) as well as the isotherms for the temperature perturbation given by the ﬁnal two
terms of equation (4.16). The standard picture of rising warmer ﬂuid and sinking cooler ﬂuid isseen to hold.
The coeﬃcients of equations (4.16) and (4.17) are√1
gzT
8 a2 + 1
α1 =(Ra − Rac ) 2µ π Ra a(4.18)

where
a=

2zT
,
L

Ra =

gzT3
,
κνπ 4

and

Rac =

(a2 + 1)3a2

(4.19)

are twice the aspect ratio of the convective domain, the Rayleigh number, and the critical
Rayleigh number for the onset of convection, respectively. The parameter α2 is deﬁned as
α2 =

gzT 1
(Ra − Rac )µ π Ra(4.20)We also have
√
1
π 8
γ1 =
κ (Ra − Rac ) 22zT (a + 1)

(4.21)

67

Streamfunction for velocity perturbation
0z (Scaled Units)

1
2
3
4
5
0

100

200
300
x (Scaled Units)

400

Isotherms of temperature perturbation
0z (Scaled Units)

1
2
Cold

Hot

Cold

3
4
5
0

100

200
300
x (Scaled Units)

400

Figure 4.3: Streamlines and isotherms for the convective perturbation of §4.2. The top panel
shows streamlines and the bottom panel displays isotherms. Both connect at lower depths. The
boundaries of the polytropic cavity of Figure 3.1 are shown as dashed lines.

68
and
√
1
aπ 8
κ (Ra − Rac ) 2 .
γ2 =2zT (a + 1)

(4.22)

From equations (4.18)–(4.22) we can see that a critical parameter governing the size of theperturbations is(Ra − Rac ) =
gzT3 1(a2 + 1)3−.
π 4 κν
a2(4.23)

To continue, we must specify three free parameters which we take to be zT , κ, and(Ra − Rac ). Mindful of the constraint that the lengthscales of the ﬂuid motions should bemuch larger than the lengthscales of the oscillation, we work with the largest possible convective
cells, ones the approximate depth of the solar convection zone, and take zT = 0.3R . These
convective motions, then, have similar dimensions to solar giant cells.
Estimates of κ and (Ra − Rac ) in the solar convection zone vary widely over manyorders of magnitude. Noting this uncertainty, we choose these parameters so that the size ofthe perturbations when compared to the ambient sound speed is small. In justiﬁcation, we notethat the velocities of giant cells in the solar convection zone are, if they exist at all, much smaller
than even those of supergranules and hence signiﬁcantly subsonic. Since we are exploring theparameter regime near the boundary of convective instability, small changes in the parameters
will have large eﬀects on the size of the perturbations.
From the discussions in §3.3 and §4.1 we can draw some general conclusions about thedependency of the frequency shifts on the parameters. We ﬁrst consider the velocity perturbations. The sinusoidal structure exhibited in equation (4.17) leads us to expect that the advective
perturbations will generate a second-order frequency shift. Since equations (4.21) and (4.22)
demonstrate that the magnitude of the vertical velocity is a ≈ 0.1 times the horizontal velocity,
we may safely neglect any shift arising from the vertical component. From the arguments of§4.1, the fractional frequency shift must be of the form
∆ω
∼−
ω



kx vx
ω

2
.

(4.24)

69
Equation (3.11) allows us to write the dependency in terms of the mode parameters, yielding∆ωk2∼ − x2 κ2 (Ra − Rac ) ∼ −κ2 (Ra − Rac ) .
ωωn + µ/2(4.25)

Again, as in §4.1, the fractional frequency shift depends quadratically on the strength of the1

perturbation (κ (Ra − Rac ) 2 ), linearly on , and inversely on n.
The frequency shift arising from the sound speed perturbation is somewhat more complicated as it combines a strong perturbation whose periodic horizontal structure ensures its
eﬀects vanish to ﬁrst order — the second term of equation (4.16) — and a weak perturbation
which will enter at ﬁrst order — the third term of that equation. After expanding equation(4.16) we ﬁnd that the shifts produced by the two terms are of equal magnitude and produce a
net downshift of∆ωRa − Rac∼−,ω
Ra(4.26)

independent of  and n. A further expansion shows that there is an additional downshift oflesser magnitude given by
Ra − Rac(n + µ/2)2 Ra − Rac∆ω∼ −z 2
∼−
,
ω
Ra
2
Ra

(4.27)

where we have used equation (3.13) to write the dependence of the depth of the acoustic cavity
on the mode parameters. Equations (4.25) and (4.26) exhibit diﬀerent dependencies on thestrength of the convective perturbation. Hence, by varying κ and keeping Ra − Rac constant,
we can choose whether the sound speed or the velocity perturbations make the predominant
contribution to the frequency shift.
We ﬁrst choose parameters which make the horizontal velocity perturbation the dominantterm and use adiabatic switching to ﬁnd the eigenfrequencies. The results are shown in Table
4.2. As expected, we see a linear dependence on  and an inverse dependence on n. Also, as
expected from equation (4.25), there is a linear dependence on (Ra − Rac ) and a quadratic
dependence on κ.

70

Table 4.2: Frequency shifts for convection cells with dominant horizontal velocity perturbations.
Some entries have been repeated in order to display trends.

na

b

κc

(Ra − Rac )d[∆ω/ω]e

1
2
3
4
5
1
1
1
1
1
1
1
1

200
200
200
200
200
100
200
300
400
500
200
200
200

1
1
1
1
1
1
1
1
1
1
1
1
1

0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.03
0.1
0.3

−1.5 × 10−2
−1.1 × 10−2
−8.7 × 10−3
−7.3 × 10−3
−6.2 × 10−3
−7.9 × 10−3
−1.5 × 10−2
−2.2 × 10−2
−3.1 × 10−2
−3.8 × 10−2
−4.9 × 10−3
−1.5 × 10−2
−4.7 × 10−2

1
1
1

200
200
200

0.5
1
2

0.1
0.1
0.1

−4.4 × 10−3
−1.5 × 10−2
−6.2 × 10−2

a Vertical mode number.
b Horizontal mode number.
c Thermal diﬀusivity. Scales strength of convective perturbation.
d See equation (4.23). Scales strength of convective perturbation.
e Relative frequency shifts from adiabatic switching.

71
As a second example we consider parameters such that the sound speed and velocity
perturbations impart approximately equal frequency shifts. The results are shown in Table 4.3.
From our analysis above, we expect the fractional frequency shifts to behave as a combinationof the dependencies in equations (4.25), (4.26), and (4.27). This is the case, as is most strikingly
seen in the dependence on  in which the fractional frequency shift reaches a minimum value. Twoeﬀects contribute to this behavior. The velocity perturbation, from equation (4.25), produces
shifts proportional to  and the sound speed perturbation, see equation (4.27), generates shifts
proportional to −2 . The location of the minimum, at  ≈ 250, arises from our particular
choices for κ and (Ra − Rac ) and can be roughly predicted when more careful track is kept of
the constants in equations (4.25) and (4.27). A similar dependence occurs for n, although inthis case the minimum occurs at n ≈ 1.5 and is hence not apparent in Table 4.3.
Finally, in the top panel of Figure 4.4 we show a sample raypath after the system has
undergone adiabatic switching. The acoustic cavity has again undergone a signiﬁcant perturbation which we have ignored in our analysis. The bottom panel shows the evolution of the rayfrequency as a function of switching time.

4.3

Compressible Plume
Solar convection only superﬁcially resembles the Rayleigh-Bénard model of the previous

section. Rather than symmetric velocity structures, rapid changes in the pressure and density
scale heights near the solar photosphere produce convective structures with broad, gentle upﬂows
and narrow, fast downﬂows. To gain some insight into the eﬀects of such structures, we consider
the interactions of p modes with a model of a cool plume descending into an adiabaticallystratiﬁed compressible envelope (Rast 1998). In the Sun such plumes are believed to be genericfeatures found at the edges of convective cells of all scales.
Pictured in Figure 4.5 are the sound-speed perturbations for the simulation we choose
as our convective model. The vertical velocities have nearly an identical structure and thehorizontal velocities are signiﬁcant only at the bottom of the plume. Although we work with

72

Table 4.3: Frequency shifts for convection cells with equal strength horizontal velocity and soundspeed perturbations. Some entries have been repeated in order to display trends.

na

b

κc

(Ra − Rac )d[∆ω/ω]e

1
2
3
4
5
1
1
1
1
1
1
1
1
1
1

200
200
200
200
200
100
200
300
400
500
600
700
800
900
1000

0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03

0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1

−8.94 × 10−4
−9.06 × 10−4
−9.25 × 10−4
−9.48 × 10−4
−9.74 × 10−3
−9.24 × 10−4
−8.94 × 10−4
−8.94 × 10−4
−9.00 × 10−4
−9.06 × 10−4
−9.10 × 10−4
−9.15 × 10−4
−9.21 × 10−4
−9.26 × 10−4
−9.30 × 10−4

a Vertical mode number.
b Horizontal mode number.
c Thermal diﬀusivity. Scales strength of convective perturbation.
d See equation (4.23). Scales strength of convective perturbation.
e Relative frequency shifts from adiabatic switching.

73

Rayleigh-Benard Raypathz (Scaled Units)

0.0

0.5

1.0

1.5

2.0
0

100

200
300
x (Scaled Units)

400

Frequency During Shift3.6240ω (Scaled Units)3.62303.62203.62103.62003.61903.6180
0.0

0.2

0.4
0.6
Switching Time (t/T)

0.8

1.0

Figure 4.4: A raypath and the evolution of the eigenfrequency for the perturbation of §4.2.
Results for a n = 1,  = 200 mode with κ = 0.05, (Ra − Rac ) = 0.1. See the text for further
discussion.

74
a snapshot of a time-dependent structure, by the time this conﬁguration has been reached theneck of the plume has nearly reached a steady-state. The majority of the interaction between
the ray and the plume is with this feature rather than the head which is still rapidly evolving
on oscillatory timescales.
Our model is one realization of Rast’s Case B1 at a time, in the non-dimensional units of
his paper, of 112.4. Rast’s Table 1 fully documents the model parameters, but we note that forthis simulation the Reynolds number is 100 and the Prandtl number is 0.1. The calculation was
performed on a stretched grid of 1024×1024 data points covering a 60×60 domain where theunits of length are measured in terms of the full width at half maximum of the applied Gaussian
temperature perturbation at the top of the atmosphere. We interpolate onto a evenly spaced
grid for our use. The polytropic index is µ = 3/2 and the acceleration of gravity is g = 5/8, tobe compared to the values of 3 and 2.7397, respectively, which were used in earlier sections.
The top boundary of Rast’s model occurs at a temperature of T = 1 where a cold Gaussian
temperature perturbation with a 30% amplitude is maintained. This level corresponds to z = 4
for a polytrope with the above stratiﬁcation and zero boundary conditions. To complete theatmosphere, we attach to the top of the simulation a µ = 3/2 polytrope extending from z = 0
to z = 4. This addition generates discontinuities in the sound speed and horizontal velocity
perturbations at z = 4, but they are primarily conﬁned to a small region near the plume center.
In order to investigate the eﬀects of the plume strength on the eigenfrequencies, we isolate
the perturbations (both velocity and sound speed) by subtracting oﬀ the ambient polytropic
stratiﬁcation. After scaling all of the perturbations by a constant factor we then add them backto the ambient stratiﬁcation.
Our choice of mode parameters is somewhat limited by the characteristics of the model.
The plume extends vertically from z = 4 to z ≈ 30 and we must choose n and  so that the
acoustic cavity encompasses these depths. Within this constraint we must also choose an  large
enough to satisfy the WKB approximation. The combination of these constraints is restrictive
and dictates our choice of n = 10 and  = 10.

75

Figure 4.5: Sound speed perturbations for the plume of §4.3. Cooler regions have slower soundspeeds and are darker in color. The background stratiﬁcation has been subtracted and the entire
60×60 computational domain is shown.

76
Unfortunately the narrow width of the plume implies that  = 10 only marginally satisﬁes
the WKB requirement. Although we are strictly violating one of the approximations we madein changing to a ray description, we suggest that employing the method of adiabatic switchingallows us, in a loose sense, to do so. Recall that we made the WKB approximation in order topass from a modal to a ray description of the eigenmodes. Since a mode is a three-dimensional
object while a ray is merely one-dimensional, the WKB approximation is necessary to ensurethat a given ray fully samples the background state. If the background varies on lengthscales
which are small compared to the wavelength of the ray, those features which are not directlyin the path of the ray will not be included in the determination of the mode eigenfrequency.
However, when we employ adiabatic switching a single ray passes through a cavity multiple
times. With an inﬁnite switching time a ray will sample every point in the cavity; for the ﬁnite
switching times we consider the ray will still explore structures on much smaller wavelengths
than is possible in any one traversal of the cavity. In essence, adiabatic switching allows the ray
to construct a detailed picture of the entire cavity and adjust its eigenfrequency accordingly.
Armed with this heuristic argument, we apply adiabatic switching despite the formally invalid
nature of the WKB approximation.
Before turning to the eigenfrequency shifts, in the top panel of Figure 4.6 we show asample raypath in the atmosphere including the plume. The plume itself is not shown, but iscentered at x = 30 and extends from z = 4 through the bottom of the displayed region. As we
can see, the raypaths are quite diﬀerent than in previous examples. Unlike Figures 4.2 and 4.4,
the raypaths do not possess well deﬁned caustics. Moreover, in several locations there appear
to be more than two rays passing through a single point.
What causes this behavior? Recall the KAM theorem, which states that for suﬃciently
large perturbations a non-integrable system will contain no invariant tori. When a trajectory
passes through a destroyed torus, it is left free to wander in phase space. We believe that the
perturbation of the plume is strong enough for that to be the case here. Since the ray is no longer
conﬁned to an acoustic cavity, given enough time it will eventually diﬀuse onto a path which

77
departs the computational domain. In the bottom panel of Figure 4.6 we show the evolution of
the frequency during the switching process. The jagged transition is also much diﬀerent than
we have seen for earlier examples and is also a consequence of the large perturbation provided
by the plume.
In Figure 4.7 we plot the dependence of the fractional frequency shift on the plumestrength which was adjusted by the scaling process discussed above. The error bars are adirect result of the diﬀusion of the ray in phase space. Since the ﬁnal state does not resideon an invariant torus, the ray undergoes Arnold diﬀusion during the switching process. When
not being switched, the eigenfrequency is a constant of the motion, however during adiabatic
switching the ray frequency is not constant, and can also drift in phase space. As we discussedin §3.2, we deal with this possibility in our computer code by beginning with several randomly
chosen initial conditions and averaging the ﬁnal eigenfrequencies. The error bars are merely
the standard deviations of the mean. Increasing the strength of the plume increases the rate ofdiﬀusion, thus explaining the trend in the size of the error bars. Presumably such errors also
occur for the non-integrable systems of §4.1 and §4.2, however, we neglected them in those cases
because they are much smaller than our quoted results. We also note that the probability thata ray will diﬀuse out of the computational domain during the switching process increases as
the perturbation strength increases. If this occurs, we ignore that ray and proceed to the next
initial condition. Hence for stronger plumes we must sample more initial conditions to achieve
the same number of completed switchings.
In Figure 4.8 we plot the perturbed quantities from a horizontal slice of the plume.
This cut is taken within the acoustic cavity of the unperturbed ray and slices at other depths
within the cavity are similar, diﬀering primarily in amplitude. As can be seen, the horizontal
velocity proﬁle possesses the symmetry necessary to make the associated frequency shifts cancel
to ﬁrst order. From our earlier arguments, we expect the vertical velocity, despite its signiﬁcanthorizontal integral, to also produce no ﬁrst-order shift. Both velocity perturbations will produce
second-order downshifts. Since the horizontal integral of the sound-speed perturbation is non-

78

Plume Raypaths
0

z

5

10

15
0

10

20

30
x

40

50

60

Frequency During Shift
3.0629ω (Scaled Units)

3.06283.06273.06263.06253.06243.0623
0.0

0.2

0.4
0.6
Switching Time (t/T)

0.8

1.0

Figure 4.6: A raypath and the evolution of the eigenfrequency for the perturbation of §4.3.
Results are for a n = 10,  = 10 mode. The plume is 0.3 times its original strength and iscentered on x = 30 in the top panel. See the text for further discussion.

79

Frequency Shift versus Plume Strength
0.0000-0.0002∆ω/ω

-0.0004-0.0006-0.0008-0.0010
0.00

0.10

0.20

0.30
0.40
Plume Strength

0.50

0.60

0.70

Figure 4.7: Fractional frequency shift versus strength of the plume of §4.3. Results are for ann = 10,  = 10 mode. The source and signiﬁcance of the error bars is discussed in the text.

80

Sound Speed Perturbation
0.000-0.010

Amplitude-0.020-0.030-0.040-0.050-0.060
0

20

40

60

40

60

40

60

x

Horizontal Velocity Perturbation
0.04

Amplitude

0.02

0.00

-0.02

-0.04
0

20
x

Vertical Velocity Perturbation
0.8

Amplitude

0.6

0.4

0.2

0.0

-0.2
0

20
x

Figure 4.8: Horizontal slices of the plume of §4.3. The slice is at z = 10 and the three plotsshow the sound speed, horizontal velocity and vertical velocity perturbations, respectively.

81
zero, we expect it to contribute a ﬁrst-order frequency downshift. For the parameters we areconsidering, the ﬁrst-order fractional sound-speed perturbation is larger than the second-order
vertical velocity perturbation and so we expect the former to be the dominant term. Since theperturbation is negative, we expect any ray interacting with this plume to undergo a ﬁrst-ordershift which scales with the strength of the plume and is always downward. This expectation is
mirrored in Figure 4.7.

4.4

Turbulent Compressible Convection
For our ﬁnal example, we consider a two-dimensional slice of a snapshot of a time-dependent, three-dimensional model of turbulent compressible convection due to Brummell et al.
(1996). It incorporates all of the features we have explored in earlier sections and is the closest
approximation to solar convection that we consider. The perturbations for this slice are shownin Figure 4.9. Plumes are seen at either end of the horizontally periodic domain while the centercontains a gentler region of upﬂow.
The complete parameters for this model are listed in Table 1 of Brummell et al. (1996) ascase R0,3 although we note that the Prandtl number is 0.1 and the root-mean-square Reynolds
number is 932. The computational domain is a plane-parallel box of 192×192×130 equally
spaced grid points covering an aspect ratio of 4:4:1 in the horizontal and vertical coordinates,
respectively. To set the length scale, the depth of the box is taken to be unity. The initialstratiﬁcation corresponds to a polytropic index of µ = 1 and an adiabatic exponent of Γ1 = 5/3,
which here is not equal to 1 + 1/µ. The acceleration of gravity is g = 20. As was the case forthe plume, the top boundary of the model lies at a ﬁxed temperature of T = 1 and supports novertical velocity perturbations.
Although the initial stratiﬁcation of the simulation corresponds to a µ = 1 polytrope,
convective motions alter the temperature structure so that in this realization it is better characterized by µ = 1.42. We use the latter value for the polytrope we attach to the top of the domain
3 A lower resolution version of the same simulation appears as Case 3 in Cattaneo et al. (1991).

82

Figure 4.9: Sound speed and velocity perturbations of the convective model of §4.4. From top
to bottom, the sound speed, horizontal, and vertical velocity perturbations in a slice at y = 0.26
in the dimensions of the model. Light colors represent hotter material and stronger ﬂows. The
true aspect ratio is 4:1.

83
in order to extend it to z = 0. Since the temperature stratiﬁcation is much larger in this model
than for the plume of §4.3, the addition is much smaller in extent, covering 0 < z < 0.12. Unlike
the plume, the top boundary for this model has a constant temperature and hence there are no
discontinuities in the sound speed or vertical velocity perturbations. The horizontal velocity is
discontinuous at z = 0.12, but the eﬀects are small. Finally, just as for the plume, we isolate,scale, and return the perturbations to the background stratiﬁcation to test the dependence ofthe frequency shift on the amplitude of the motions.
The strong temperature stratiﬁcation allows us be somewhat more lenient in choosing
the mode parameters than we were in §4.3. Due to the the highly turbulent nature of theconvection we have almost no hope of strictly satisfying the WKB approximation for all scalesof the convective structures. However, the argument made in §4.3 suggests that adiabatic
switching will allow us to sample smaller structures than might be expected. As a compromise
between sampling a signiﬁcant portion of the cavity and satisfying the WKB condition, we workwith an n = 1,  = 10 mode. A sample raypath and the evolution of the switched frequency are
shown in Figure 4.10. Qualitatively, the results are similar to those presented in §4.3 leading us
to conclude that the perturbation is strong enough to destroy all of the invariant tori.
We plot horizontal slices of the plume showing the perturbed quantities at a depth ofz = 0.3 in Figure 4.11. Unlike the plume of §4.3, slices at other depths diﬀer both in amplitude
and in structure from those shown. However, we can still make some general statements about
the expected frequency shifts. The horizontal velocity perturbation always has a nearly zero
horizontal integral implying that it contributes a second-order downshift. For reasons we havepreviously discussed (see §4.1), the vertical velocity perturbation, regardless of form, will also
contribute a second-order downshift. The form of the sound speed perturbation varies from
depth to depth and we expect that over the entire convective domain it will average to zero.
However, for the particular acoustic cavity we consider, the integral of the perturbation ispositive. Hence, it will produce a ﬁrst-order frequency upshift, although other modes occupying
other acoustic cavities will undergo frequency downshifts.

84

Raypaths in Turbulent Convection
0.00
0.10

z

0.20
0.30
0.40
0.50
0

1

2
x

3

4

Frequency During Shift
19.2380

ω (Scaled Units)19.2370
19.2360
19.2350
19.2340
19.2330
0.0

0.2

0.4
0.6
Switching Time (t/T)

0.8

1.0

Figure 4.10: A raypath and the evolution of the eigenfrequency for the perturbation of §4.4.
Results are for an n = 1,  = 10 mode with the convection at 0.1 times its original strength. See
the text for further discussion.

85Sound Speed Perturbation
0.08

Amplitude

0.06

0.04

0.02

0.00

-0.02
0

1

2
x

3

4

3

4

3

4

Horizontal Velocity Perturbation
1.0

Amplitude

0.5

0.0

-0.5

-1.0
0

1

2
x

Vertical Velocity Perturbation
0.6
0.4

Amplitude

0.2
0.0
-0.2
-0.4
-0.6
-0.8
0

1

2
x

Figure 4.11: Horizontal slices of the plume of §4.4. The slice is at z = 0.3 in the units ofthe model and the three plots show the sound speed, horizontal velocity, and vertical velocity
perturbations, respectively.

86
The net frequency shift is a combination of these diﬀerent shifts. A rough order-ofmagnitude estimate suggests that the velocity perturbations contribute the dominant terms.
This is not surprising since the Prandtl number (the ratio of the viscosity to the thermal diﬀusivity) of this simulation is less then unity — in this case, 0.1. Low Prandtl numbers imply thattemperature perturbations (and hence sound speed perturbations) are smoother and weaker
than velocity perturbations. Since the velocity perturbations contribute the dominant terms,
we expect the fractional frequency shift to be downward and quadratic in the strength of theperturbation.
In Figure 4.12 we plot the dependence of the fractional frequency shift on the strengthof the convection for an n = 1,  = 10 mode. The error bars again arise from averaging theresults from diﬀerent initial conditions and signal the presence of destroyed phase-space tori.
The frequency shift, while always negative, scales approximately as 1.5 , where  is the strengthof the perturbation. This dependency is likely a combination of the linear frequency shift arisingfrom the sound speed perturbation and the quadratic perturbation arising from the advective
motions.

87

Frequency Shift versus Convective Strength
0.000-0.002-0.004∆ω/ω

-0.006-0.008-0.010-0.012-0.014-0.016
0.00

0.10

0.20
0.30
Convective Strength

0.40

0.50

Figure 4.12: Fractional frequency shift versus strength of the convective model of §4.4. The
results are for a mode with n = 1,  = 10.

Chapter 5

Conclusions

Adiabatic switching provides a powerful tool for predicting the eﬀects of convective structures on p-mode eigenfrequencies. In §5.1 we summarize a few of our key results and in §5.2 weoutline several possibilities for future applications and extensions.

5.1

Results
In the ray approximation, p modes are Hamiltonian systems. The canonical positions

and momenta are the position and wavenumber, respectively, while the dispersion relation is
the functional form of the Hamiltonian. We have shown that time-independent, horizontallyinvariant convective models form integrable systems. Their eigenfrequencies can be predictedwith adiabatic switching as well as an independent method, semi-classical EBK quantization.
The results from the two methods agree. However, we have also shown that adiabatic switchingis a viable method of ﬁnding the eigenfrequencies of non-integrable systems (in this context,
those with horizontal variations), a domain in which EBK quantization can only oﬀer guidance
and not precise results.
From our investigations we draw some general conclusions concerning the eﬀects of convection on p-mode eigenfrequencies. Foremost, the fractional frequency shift is the weighted
integral of the perturbation over the resonant cavity of the eigenmode. For cases such as thesound-speed perturbation of §3.3 and the Doppler shift of §3.4 where the perturbation is monotonic over the acoustic cavity, the shift is linear in the perturbation strength.

89
Of more physical interest are examples where the perturbations oscillate in sign. For
horizontally antisymmetric perturbations, such as those of §4.1 and §4.2, we have argued on both
physical and mathematical grounds that the frequency shift vanishes to ﬁrst order. Adiabatic
switching conﬁrms this result. Furthermore, it demonstrates that a second-order perturbation
does exist, resulting from the ray spending extra time in regions of retarding material. We havealso shown that vertical velocity perturbations always produce second-order fractional frequency
shifts, independent of the perturbation structure.
We also treated even more complex perturbations such as a plume in §4.3 and a twodimensional slice of a convective model in §4.4. In these complex models, over large enoughscales, we expect every perturbation to average to zero — cold, narrow plumes are balancedby broad, warm upﬂows. If p-mode cavities extended over large scales we would always expect
second-order downshifts as any ﬁrst-order shifts will cancel. However, acoustic cavities have alimited vertical extent and hence may sample a region in which the perturbations do not fully
cancel. In that case, the dominant perturbation cannot necessarily be predicted. But, it is
known that solar convection is characterized by low Prandtl numbers, implying that velocity
perturbations are generally larger in amplitude than sound speed perturbations. Within atypical acoustic cavity, velocity perturbations will be the dominant eﬀect. In the majority ofcases, then, we expect second-order frequency downshifts arising from velocity perturbations.
For the simpler perturbations of §4.1 and §4.2 we investigated the dependence on radialmode number n and angular degree . First-order velocity perturbations vary roughly as√√and 1/ n and second-order perturbations as the square of these quantities. For sound speed
perturbations whose amplitude increases with depth, the fractional frequency shifts generally
increase with n and decrease with . The exact dependency of the shift, however, depends onthe functional form of the perturbation.
In addition, we have demonstrated that even mild perturbations can signiﬁcantly change
the structure of the resonant cavity of an eigenmode despite producing only small frequency
shifts. As the perturbation strength increases for realistic convective models, the ray trajectories

90
diﬀuse and are no longer conﬁned within an obvious cavity. It is an open question as to what
extent this occurs in the Sun, particularly as results from time-distance helioseismology imply
that ray-like entities which follow predictable trajectories do exist.
It is intriguing to note, as discussed in Gough et al. (1996) and seen in Figure 1.1, thatcurrent solar models predict larger sound speeds in the outer portion of the convection zone thanare suggested by p-mode observations. Since p-mode eigenfrequencies are proportional to thesound speed, the implication is that theory predicts higher eigenfrequencies than are actually
produced by the Sun. We have shown in this work, however, that one eﬀect of convective structures is to lower p-mode eigenfrequencies. This eﬀect is not accounted for in current theoreticalsolar models and is a potential source of the discrepancy.
Not only is the eﬀect of the proper sign, it has the proper magnitude to explain thediﬀerence. As we saw in Chapter 4 the scaling of the frequency shift from convection depends
on the integral of the perturbation over the cavity. To make a conservative estimate, we assume
that any frequency shift is second-order in the velocity:
−

∆ω
v2
∼ 2.
ω
c

(5.1)

Evaluating c at the midpoint of an acoustic cavity gives− v 2
v2
∆ω
∼
∼
,
ω
gR
400

(5.2)

where v is the convective velocity in km/s,  is the mode degree, and we have suppressed
terms of order 1. For typical supergranular photospheric outﬂow velocities of 0.5 km/s, thefractional frequency shift is ≈ 3 × 10−4. This result is a conservative lower limit: higher  values,
higher velocities, and lower-order frequency shifts will all produce larger shifts. Observational
determinations of p-mode eigenfrequencies quote error bars of ∆ω/ω ≈ 3 × 10−4 . Hence weconclude that convective structures will, at the minimum, produce detectable shifts.

91

5.2

Future Work
Although adiabatic switching is an attractive method for ﬁnding p-mode eigenfrequencies,

care must be taken to apply it to problems for which it is valid. In the Sun, the true excited
objects are global modes. But, inherent in our methodology is the assumption that in the smallwavelength limit these modes may be approximated as one-dimensional rays. We did not discuss
the transition between mode and ray. However, one may form a ray by superposing eigenmodes
into a wavepacket whose dynamics are then governed by the ray equations (2.38). Bogdan (1997)
has shown that, due to the ﬁnite number of excited solar p modes, a truly one-dimensional ray
cannot be formed. Instead, the smallest wavepacket which can be constructed is ≈ 30 Mm insize. This dimension, characteristic of supergranules, suggests a lower size limit to the structures
which can be explored using adiabatic switching. Unfortunately this limit is somewhat more
stringent than the WKB criteria, which sets the lower size limit at ≈ 400/ Mm. It may bepossible, as was discussed in §4.3, to show that the method of adiabatic switching loosens these
constraints. Although it is not clear how to frame this argument in a rigorous manner, further
investigation into the limitations of the WKB approximation is certainly warranted. In any case,the smallest convective structures cannot strictly be treated within our present formulation of
adiabatic switching. As discussed in §1.3, other methods, such as the inclusion of a turbulent
pressure in the governing equations, must be employed.
A conceptually simple, although perhaps computationally expensive, extension to ourwork would be the incorporation of a third dimension into the plane-parallel geometry. The
system would gain another degree of freedom, however the reference state would remain integrable since another constant of the motion (ky ) would also be added. The transition fromtwo to three dimensions is accompanied by two important geometric eﬀects. As an exampleof the ﬁrst, consider the plume of §4.3. In two dimensions it is impossible for rays to avoid
the plume; with every traversal of the cavity they pass through the heart of the disturbance.
In three dimensions a one-dimensional ray can easily bypass a plume-like structure. Second,

92
convection in three dimensions is vastly diﬀerent than in two dimensions. Even highly turbulent
two-dimensional convection is organized in cell-like structures such as those of §4.2. Addition ofa third dimension results in a much more complicated system. In §4.4 we attempted to accountfor the second of these diﬀerences by taking two-dimensional slices of a three-dimensional model,
but a full treatment involves extending our method to include the propagation of rays in three
dimensions.
Working in a spherical, rather than plane-parallel, geometry would complement the addition of a third dimension. The solar convection zone extends to a depth of ≈ 0.7R and so
we expect the eﬀects of curvature to be small. However, we would be able to explore recent
spherical shell simulations of the convection zone (Elliott et al. 1998, for example).
Both rotation and magnetic ﬁelds can be expressed as additional terms in the governing
dispersion relation (Gough 1993). The primary eﬀect of rotation is known: it breaks the sphericalsymmetry and introduces a dependence of the frequency on m, the azimuthal order in thespherical harmonic decomposition. Magnetic eﬀects are not as well understood. One of the mostintriguing, as well as most diﬃcult, questions is the eﬀect of active regions on the eigenfrequency
spectrum. Magnetic ﬁeld strengths reach equipartition in sunspots, becoming strong enough tosuppress convective motions as well as have signiﬁcant impacts on ray propagation. We can also
investigate the signatures of a buried magnetic ﬁeld on the p-mode eigenfrequencies. Detectionsof such signals would permit predictions of the appearance of magnetic ﬂux at the solar surface.
Finally, power spectra of the solar oscillations show that p modes possess intrinsic line
proﬁles, the source and shape of which is not completely understood. Previous works (Kumar
et al. 1994; Roxburgh & Vorontsov 1995; Rast & Bogdan 1998) have treated the eﬀects of
intrinsic damping, noise, and source structure on the line proﬁles. Intriguingly, the error barsof Figures 4.7 and 4.12 can also be interpreted as linewidths. They arise from the drift of theeigenfrequency in phase space as the ray passes through time-dependent convective structures(in this case, during the switching process). In the Sun, time-dependent convective motions,even those lasting for timescales which are long compared to a p-mode period, will broaden the93mode eigenfrequencies via the process of Arnold diﬀusion discussed in §2.5. A closer explorationof the error bars we discuss in §4.3 and §4.4 will likely yield insight into the nature of p-mode
linewidths.

Bibliography

Arnold, V. I. 1978, Mathematical Methods of Classical Mechanics (New York: Springer-Verlag)
Baker, N. H. & Gough, D. O. 1979, ApJ, 234, 232
Balmforth, N. J. 1992, MNRAS, 255, 632
Beck, J. G., Duvall Jr, T. L., & Scherrer, P. H. 1998, Nature, 394, 653
Bogdan, T. J. 1997, ApJ, 477, 475
Brillouin, L. 1926, J. Phys. Radium, 7, 353
Brown, T. M. 1984, Science, 226, 687
Brummell, N., Cattaneo, F., & Toomre, J. 1995, Science, 269, 1370
Brummell, N. H., Hurlburt, N. E., & Toomre, J. 1996, ApJ, 473, 494
Brut, F. & Arvieu, A. 1993, J. Phys. A, 4749
Cattaneo, F., Brummell, N. H., Toomre, J., Malagoli, A., & Hurlburt, N. H. 1991, ApJ, 370,
282
Chandrasekhar, S. 1987, Ellipsoidal Figures of Equilibrium (New York: Dover Publications,Inc.), 28–30
Christensen-Dalsgaard, J. 1980, MNRAS, 190, 765
Christensen-Dalsgaard, J. 1997, Lecture Notes on Stellar Oscillations, 4th edn., available at
http://www.obs.aau.dk/∼jcd/oscilnotes/
Cowling, T. G. 1941, MNRAS, 101, 367
Cox, J. P. 1980, Theory of Stellar Pulsation (Princeton: Princeton University Press)
Duvall, Jr., T. L., Kosovichev, A. G., & Murawski, K. 1998, ApJ, 505, L55
Ehrenfest, P. 1917, Philos. Mag., 33, 500
Einstein, A. 1917, Verhandl. Deut. Physik. Ges., 19, 82, as translated in Engel (1997)
Elliott, J. R., Miesch, M. S., Toomre, J., Clune, T., & Glatzmaier, G. A. 1998, in Structure
and Dynamics of the Interior of the Sun and Sun-like Stars, ed. S. Korzennik & A. Wilson(Noordwijk: ESA), 765–770
Engel, A. 1997, The Collected Papers of Albert Einstein, Vol. 6, The Berlin Years: Writings,
1914–1917 (Princeton: Princeton University Press), 434–443

95
Fetter, A. F. & Walecka, J. D. 1980, Theoretical Mechanics of Particles and Continua (NewYork: McGraw-Hill), 434–458
Ghosh, P., Antia, H. M., & Chitre, S. M. 1995, ApJ, 451, 851
Goldstein, H. 1965, Classical Mechanics (Reading: Addison-Wesley)
Gough, D. O. 1993, in Dynamique des Fluides Astrophysiques: Astrophysical Fluid Dynamics,
Les Houches No. XLVII (New York: North-Holland), 399–560
Gough, D. O., Kosovichev, A. G., Toomre, J., et al. 1996, Science, 272, 1296
Gruzinov, A. V. 1998, ApJ, 498, 458
Hénon, M. & Heiles, C. 1964, AJ, 69, 73
Johnson, B. R. 1985, J. Chem. Phys., 83, 1204
Keller, J. B. 1958, Ann. Phys. (N.Y.), 4, 180
Keller, J. B. & Rubinow, S. I. 1960, Ann. Phys. (N.Y.), 9, 24
Kumar, P., Goldreich, P., & Kerswell, R. 1994, ApJ, 427, 483
Lamb, H. 1932, Hydrodynamics, 6th edn. (Cambridge University Press)
Lavely, E. M. & Ritzwoller, M. H. 1992, Phil. Trans. R. Soc. Lond. A, 339, 431—. 1993, ApJ, 403, 810
Leibacher, J. & Stein, R. F. 1971, Astrophys. Lett., 7, 191
Leighton, R. B., Noyes, R. W., & Simon, G. W. 1962, ApJ, 135, 474
Lichtenberg, A. J. & Lieberman, M. A. 1983, Regular and Stochastic Motion (New York:
Springer-Verlag)
Murawski, K. & Goossens, M. 1993, A&A, 279, 225
Murawski, K. & Roberts, B. 1993a, A&A, 272, 595—. 1993b, A&A, 272, 601
Nöckel, J. U. & Stone, A. D. 1997, Nature, 385, 45
Noid, D. W. & Marcus, R. A. 1975, J. Chem. Phys., 62, 2119
Patterson, C. W. 1985, J. Chem. Phys., 83, 4618
Rast, M. P. 1998, J. Fluid Mech., 369, 125
Rast, M. P. & Bogdan, T. J. 1998, ApJ, 496, 527
Rosenthal, C. S. 1997, in SCORe’96 : Solar Convection and Oscillations and their Relationship,
145–160
Rosenthal, C. S., Christensen-Dalsgaard, J., Kosovichev, A. G., Nordlund, Å., Reiter, J., Rhodes
Jr, E. J., Schou, J., Stein, R. F., & Trampedach, R. 1998, in Structure and Dynamics of theInterior of the Sun and Sun-like Stars, 521–525
Roxburgh, I. W. & Vorontsov, S. V. 1995, MNRAS, 272, 850

96
Shampine, L. F. & Gordon, M. K. 1975, Computer Solution of Ordinary Diﬀerential Equations:
The Initial Value Problem (San Francisco: W. H. Freeman)
Shirer, H. N. 1987, in Lecture Notes in Physics, Vol. 271, Nonlinear Hydrodynamic Modeling:
A Mathematical Introduction, ed. H. N. Shirer (New York: Springer-Verlag), 22–46
Skodje, R. T. & Cary, J. R. 1988, Comput. Phys. Rep., 8, 221
Smith, R. S., Shirts, R. B., & Patterson, C. W. 1987, J. Chem. Phys., 86, 4452
Spiegel, E. 1971, ARA&A, 9, 323
Stix, M. & Zhugzhda, Y. D. 1998, A&A, 335, 685
Swisdak, M. & Zweibel, E. 1999, ApJ, 512, 442
Tabor, M. 1989, Chaos and Integrability in Nonlinear Dynamics: An Introduction (New York:
Wiley and Sons)
Ulrich, R. K. 1970, ApJ, 162, 993
Zhugzhda, Y. D. 1998, A&A, 332, 314
Zhugzhda, Y. D. & Stix, M. 1994, A&A, 291, 310

Appendix A

Table of Important Symbols

Table A.1: Important symbols. Referenced pages give the location of either the ﬁrst appearance
or a deﬁnition.

SymbolPage

Description

∇2h

15

Horizontal Laplacianˆ

17

Unit vector102

Eulerian perturbationa66

Aspect ratioa101Marker identifying a ﬂuid element

A

21

Varying plane-wave amplitude

B1 , B2 , B3 , B4

110,110,110,111

Functions of χc

13

Adiabatic sound speedc045

Polytropic adiabatic sound speed

C

25

Integration curvef

13

Force per unit mass

g, g

17,13

Scalar and vector gravitational acceleration

G

13

Newtonian gravitational constanth23

Planck’s constantH, H0 , H1

24, 35, 35

General, integrable, and non-integrable HamiltonianHρ

19

Density scale heightI, Ix , Iz

25, 42, 40

General, horizontal, and vertical action variable ofaction-angle coordinates

kr , kφ , kx , kz

28, 28, 40, 40

Radial, azimuthal, horizontal, and vertical wavenumbersk21

Wavevector
continued on next page98

Table A.1 continued
SymbolPage

Description

K

48

Perturbation wavenumber

16

Spherical harmonic degree

L

57

Horizontal dimension of computational domain

L1 , L2 , L3112

Diﬀerential operatorsm

16

Azimuthal order of spherical harmonicn

16

Radial mode numberN

18

Brunt-Väisälä frequencyp13

Pressurep

23

Canonical momentumq

23

Canonical position

Q

101

Generic variabler4

RadiusR

3

Solar radiusRa, Rac66

Rayleigh number, critical Rayleigh numbers28

Number of encountered caustics

S

31

Action

t

12

Time

T

35

Switching time

Tx , Tz43, 43

Periods of horizontal and vertical motion

vx , vz

46, 46

Horizontal and vertical velocityv

12

Velocityx

39

Horizontal coordinatex21

Position vector

Ym

16

Spherical harmonicz

39

Vertical coordinate

z1 , z2

43

Minimum and maximum depths of acoustic cavity

zmid

43

Average depth of acoustic cavity

zT

66

Depth of convective cell

α1 , α266, 66

Sound speed perturbation amplitudes

β18

Adiabatic stratiﬁcation parameter

γ1 , γ2

66, 68

Horizontal and vertical velocity perturbation amplitudes

Γ1

13

Adiabatic exponent
continued on next page

99

Table A.1 continued
SymbolPage

Descriptionδ102

Lagrangian perturbation∆

31

Generic diﬀerence

26

Small parameterθ

15

Polar angle

Θ25

Angle variable of action-angle coordinates

κ

66

Thermal diﬀusivity

λ35

Switching functionΛ

21

Ratio of lengthscalesµ

39

Polytropic index

ν66

Kinematic viscosityνk ,νx , νz25, 42, 42

Generic, horizontal, and vertical characteristic frequencies of motionξ

17

Radial component of the Lagrangian displacementρ12

Densityφ

15

Azimuthal angle

ϕ21

Varying plane-wave phaseΦ13Gravitational potentialχ

16

Divergence of the Lagrangian displacement

ψ57

StreamfunctionΨ

19

Scalar variable in wave equation

Ψ0

21

Constant plane-wave amplitudeω

21

Frequency

ω̃

22

Local frequency

ωc19

Acoustic cutoﬀ frequencyΩ

46

Polytropic frequency

Appendix B

Fluid Descriptions and Perturbations

At suﬃciently small scales all matter is discrete, separated into atoms and molecules.
However, for a large class of substances, termed ﬂuids, this discretization is not apparent. Instead, the mean free path of atomic and molecular motion is much smaller than any other relevant
lengthscale (for example, the density scale height or the diﬀusion length). This separation oflengthscales allows one to average over many particles and thus approximate a fundamentallydiscrete object as a continuous medium. In this appendix, we outline the Eulerian and Lagrangian descriptions of a ﬂuid. While both are valid representations, most of ﬂuid dynamics,
including the standard forms of the ﬂuid equations, is conceived and written using the Eulerian description. However, as is seen below, the Lagrangian description occasionally provides
simpliﬁcations. We also describe the Eulerian and Lagrangian perturbations (which are closely
related to their respective ﬂuid descriptions) and demonstrate their commutation properties
with respect to three mathematical operators.
As their names suggest, these descriptions are not new. Instead of breaking new ground,
this appendix is merely our own organization of the material. Other treatments may be foundin Chandrasekhar (1987) and in Chapter 5 of Cox (1980).

B.1

Descriptions
The choice of independent variables diﬀerentiates the Eulerian and Lagrangian ﬂuid de-scriptions. The position vector x and time t are the independent variables in an Eulerian ﬂuid

101
and any quantity Q is written as Q(x, t). This representation is completely general since, in
principle, the value of a variable at any point is uncorrelated with the value at a neighboring
point (ignoring, for the moment, constraints imposed by the ﬂuid equations).
In contrast, the Lagrangian description divides a ﬂuid into tiny parcels. The independent
variables are the time t and a quantity which uniquely identiﬁes a particular ﬂuid element. Often
this quantity is taken to be the position of the ﬂuid element at some reference time, but otherchoices are possible. To acknowledge this ambiguity, I follow the notation of Cox (1980) anddenote the identifying marker as a. Variables are associated with parcels of ﬂuid, not points inspace, and are denoted by Q(a, t).
By supposing the existence of ﬂuid parcels, the Lagrangian description appears, by theapplication of Ockham’s razor, to be inferior to the Eulerian. The need for ﬂuid parcels requires,at some level, that the ﬂuid be smooth enough for neighboring points to share the same physicaldescription; the Eulerian picture does not require this precondition. To resolve this apparent
diﬃculty, recall that, by deﬁnition, every ﬂuid possesses a lengthscale where averaging over many
particles is permissible. The Lagrangian description is valid when this lengthscale is larger than
the size of a ﬂuid parcel. As both descriptions are equally valid, other considerations will dictate
which to use for a given problem.

B.2

Mathematical Operations∂
Due to the choice of independent variables, the usual gradient operator ∇ = ∂x
is auseful operation in the Eulerian picture. This is not the case in the Lagrangian description
where position is not an independent variable; instead x = x(a, t). The appropriate spatial∂.
derivative is then ∇a = ∂a∂As time is exclusively an independent variable, ∂t
is useful in both the Lagrangian andEulerian pictures. However, it is not the same useful operator in both descriptions. In the

102
Eulerian frame,
∂
∂
,
=
∂t
∂t x

(B.1)

∂
∂
=
,
∂t
∂t a

(B.2)

while in the Lagrangian frame,with the notation |j indicating that j is to be held constant during diﬀerentiation. For theEulerian description the time derivative is taken at a point while in the Lagrangian description
it is taken with respect to a particular ﬂuid parcel.
The derivative of equation (B.2) is given a special name: the advective, or Stokes, derivaD 1
tive. In the Eulerian description it is often denoted Dt. Since the derivative is taken withrespect to a ﬂuid parcel, the position varies with time. Thus for any Eulerian variable Q associated with the element:
∂
∂Q
D
∂Q
∂x
Q(x, t) ≡Q(x(a, t), t) =+· ∇Q =+ v · ∇Q
Dt
∂t a
∂t x
∂t
∂t

(B.3)

where v is the ﬂuid velocity. In short, for derivatives taken along the motion of a ﬂuid parcel
D
∂
∂
=
=
+ v · ∇.
Dt
∂t a
∂t x

B.3

(B.4)

Perturbations
Of course, ﬂuids do not always occupy equilibrium conﬁgurations. Here we discuss theapplication of the Eulerian and Lagrangian descriptions to ﬂuid states which are slightly perturbed away from an reference state. Linear perturbation theory postulates that a variable Q
in the perturbed state may be expressed as
Q = Q0 + ∆Q

(B.5)

where Q0 is the variable in the reference state and ∆Q is a generic perturbation which is small
compared to Q0 .
1 Occasionally the more confusing d is used.
dt

103
Just as both the Eulerian and Lagrangian descriptions can be used to represent ﬂuids,
the perturbation of equation (B.5) can be written in diﬀerent ways. We will discuss two forms,
the Eulerian and Lagrangian perturbations, which are closely related to their respective ﬂuid
descriptions. As we will see, both perturbations can be considered in either the Eulerian or
Lagrangian description. Unfortunately, the literature has not settled on standard symbols forthe perturbations. We follow the notation of Cox (1980) and signify the Eulerian and Lagrangian
perturbations to a quantity Q by Q and δQ, respectively.
Figure B.1 illustrates both the Eulerian and Lagrangian perturbations. Position, normallya three-dimensional vector, has been collapsed to one dimension and plotted on the vertical axis;
time is shown on the horizontal axis. The slanted lines are streamlines for three parcels of ﬂuid.
The dotted lines are streamlines in the reference state (denoted by y0 ) while the solid line is a
streamline in the perturbed state (denoted by y). Each streamline is further delineated by the
position at time t1 of a ﬂuid parcel (serving as the a in the Lagrangian description). Thus, inthe reference state, a ﬂuid parcel which occupied position x2 at time t1 follows the streamline
labeled y0 (x2 ). Finally, representative parcels of ﬂuid are shown on the streamlines. So, y0 (x2 )and y(x2 ) are streamlines for the same ﬂuid parcel, the ﬁrst tracing the path of the parcel in thereference state, and the second in the perturbed ﬂuid. Although the reference and perturbed
states are shown occupying the same positions at time t1 , in general this need not be the case.
Note that the streamlines, although linear in this sketch, are, in general, curved.
An Eulerian perturbation to a variable Q is the diﬀerence at the same point and time
between the perturbed and reference ﬂuid states. We ﬁrst treat this perturbation in the Eulerian description. Suppose the two ﬂuid parcels following the streamlines y0 (x4 ) and y(x2 ) aredescribed by the variable Q. They are diﬀerent ﬂuid parcels, as indicated by their diﬀerent
positions at time t1 , yet at time t2 , they occupy the same position, the latter in the perturbed
and the former in the reference model. At this location and moment, the Eulerian variation isQ (x = x3 , t = t2 ) = Q(x = x3 , t = t2 ) − Q0 (x = x3 , t = t2 ),(B.6)

Position

104

x4
B
y(x2)
x3
A

y0(x4)x2
Ax1
y0(x2)

t1

t2

Time

Figure B.1: A pictorial demonstration of the Eulerian and Lagrangian perturbations. See thetext for details.

105
where we have included the independent variables to make it clear we are in the Eulerian
description. In general form:Q (x, t) = Q(x, t) − Q0 (x, t).

(B.7)

From equation B.7, we see that x = 0.
A Lagrangian perturbation is the diﬀerence between identical ﬂuid parcels, one in thereference state and one in the perturbed system. This perturbation is most naturally describedin the Lagrangian description (see below), but it proves useful to examine it in the Eulerian
description as well. In the speciﬁc case of Figure B.1, we have:δQ(x = x3 , t = t2 ) = Q(x = x3 , t = t2 ) − Q0 (x = x1 , t = t2 ).

(B.8)

Note the distinction between x3 and x1 in the ﬁnal two terms; it forms the crucial diﬀerence
between the Eulerian and Lagrangian perturbation. In the reference state the variable is evaluated at x = x1 because this is the location at t = t2 of the reference state parcel which, at
t = t1 , coincided with the parcel in the perturbed state.
Equation (B.8) gives δx = x3 − x1 . In general δx depends on both the time and positionof the observation, so extension to the general case yields
δQ(x, t) = Q(x, t) − Q0 (x − δx(x, t), t),(B.9)

where we again include the independent variables to make it clear these results apply in the
Eulerian description.
As both equations (B.7) and (B.9) are written in the Eulerian description, we may subtract to ﬁndδQ − Q = Q0 (x, t) − Q0 (x − δx(x, t), t) ≈ δx ·∂Q0 = δx · ∇Q0 ,∂x(B.10)

where the approximation is made to ﬁrst-order in smallness. At this level of accuracy, Q0 may
be replaced by Q.
Although the Lagrangian description is not as intuitive, it is possible to express both
perturbations in its formalism. Recall that the independent variables are no longer time and106
position, but time and a marker identifying the ﬂuid parcels (in the case of Figure B.1, initial
position). The Lagrangian perturbation in the Lagrangian description is, from Figure B.1,
δQ(a = x2 , t = t2 ) = Q(a = x2 , t = t2 ) − Q0 (a = x2 , t = t2 ).

(B.11)

Or, in general form
δQ(a, t) = Q(a, t) − Q0 (a, t).

(B.12)

Note the symmetry between this expression and the Eulerian perturbation in the Eulerian
description of equation B.7. A simple application of equation (B.12) demonstrates δa = 0.
Finally, consider an Eulerian perturbation in the Lagrangian description. From Figure
B.1
Q (a = x2 , t = t2 ) = Q(a = x2 , t = t2 ) − Q0 (a = x4 , t = t2 ).

(B.13)

Since the Eulerian perturbation (usually) compares diﬀerent ﬂuid parcels, the two terms on theright-hand side of equation (B.13) are evaluated at diﬀerent a. In order to generalize this result,notice that a = x2 − x4 in this example. In general, a is a function of both time and the initial
position of the ﬂuid parcel. So,Q (a, t) = Q(a, t) − Q0 (a − a (a, t), t).

(B.14)As might be expected, this result is quite similar to equation (B.9). Evaluating the diﬀerence
between equations (B.12) and (B.14) shows thatδQ − Q = Q0 (a − a (a, t), t) − Q0 (a, t) ≈ −a ·∂Q0 = −a · ∇a Q0∂a(B.15)

Again, Q may be substituted for Q0 in the last term without aﬀecting the result.

B.4

Commutation Relations
Using the formalism of the preceding sections, we examine the commutation propertiesof the perturbations with three mathematical operators. In particular, we demonstrate that[∇,  ] =



 

∂
D
, =
,δ =0
∂t x
Dt

(B.16)

107and
[∇, δ] = 0,




∂
D 
, = 0,
,δ =
 0,
∂t x
Dt

(B.17)

where the commutation symbol [ , ] is deﬁned for two operators F and G as [F ,G] ≡ F G − GF .
Although we can easily extend the results to treat the operator ∇a , it is rarely used in theliterature and so we neglect it here.

B.4.1

Eulerian Perturbation, Commuting Operators

First, consider the Eulerian perturbation. It is simple to prove the two equalities ofequation (B.16) using the Eulerian description of equation (B.7). Applying the two operators
to the perturbation yields∇(Q ) = ∇Q − ∇Q0

and

∂
∂
∂
(Q ) =
Q−
Q0 ,
∂t x
∂t x
∂t x

(B.18)

while applying the perturbation to the operators gives(∇Q) = ∇Q − ∇Q0

and

∂Q
∂t x


=

∂Q
∂Q0
−
.
∂t x
∂t x

(B.19)The right-hand sides of equations (B.18) and (B.19) are clearly identical, thus proving the ﬁrst
two commutation relations.

B.4.2

Eulerian Perturbation, Non-Commuting Operators

To prove the inequality of equation (B.17), recall from equation (B.4) that the advective
derivative is the same as a derivative at constant a. In light of this fact, it is easiest to consider
the Eulerian perturbation in the Lagrangian description. Expanding equation (B.14) to ﬁrstorder, the Eulerian perturbation becomesQ (a, t) = Q − Q0 + a · ∇a Q0 .

(B.20)

So, remembering that a = a (a, t)∂∂∂(Q ) =
Q−
Q0 +
∂t a
∂t a
∂t a



∂
a
∂t a



· ∇a Q0 + a ·∂(∇a Q0 ).
∂t a(B.21)

108
Applying the operations in the opposite order gives


∂Q
∂t a


=

∂Q
∂Q0
∂Q0
−
+ a · ∇a
.
∂t a
∂t a
∂t a

(B.22)

Switching the order of diﬀerentiation in the last term of either equation (B.21) or (B.22), subtracting, and changing to the notation of equation (B.4)





D
DQ
∂a
D 
, =
(Q ) −=· ∇a Q0 .
Dt
Dt
Dt
∂t a

(B.23)

In general, the Eulerian perturbation and the advective derivative do not commute.

B.4.3

Lagrangian Perturbation, Commuting Operators

Next, consider the Lagrangian perturbation δ. To prove the equality of equation (B.16),
use equation (B.4) and the form of the perturbation given in equation (B.12).
∂∂∂(δQ) =
(Q) −
(Q0 )
∂t a
∂t a
∂t a

(B.24)

and

δ

∂Q
∂t a


=

∂Q
∂Q0
−
.
∂t a
∂t a

(B.25)The right-hand sides of equations (B.24) and (B.25) are equal, proving that the Lagrangian
perturbation commutes with the advective derivative.

B.4.4

Lagrangian Perturbation, Non-Commuting Operators

The inequalities of equation (B.17) are best treated in the Eulerian description. Expanding equation (B.9) to ﬁrst-order,δQ(x, t) = Q − Q0 + δx · ∇Q0 .

(B.26)

Working ﬁrst with the ∇ operator and recalling that δx = δx(x, t),∇(δQ) = ∇Q − ∇Q0 + (δx · ∇)∇Q0 + (∇Q0 · ∇)δx + ∇Q0 ×(∇×δx).

(B.27)

109
The ﬁnal three terms result from applying a vector identity to ∇(δx · ∇Q0 ). Applied in thereverse order, the operations produceδ(∇Q) = ∇Q − ∇Q0 + (δx · ∇)∇Q0 .

(B.28)[∇, δ] = ∇(δQ) − δ(∇Q) = (∇Q0 · ∇)δx + ∇Q0 ×(∇×δx).

(B.29)

Hence,

Finally, consider the commutation properties of the Lagrangian perturbation and thepartial with respect to time at constant x.
∂∂∂(δQ) =
Q−
Q0 +
∂t x
∂t x
∂t x




∂
∂
δx · ∇Q0 + δx ·
∇Q0 ,
∂t x
∂t x

(B.30)

but

δ

∂Q
∂t x


=

∂Q
∂Q0
∂Q0
−
+ δx · ∇
∂t x
∂t x
∂t x

(B.31)

Again, by switching the order of diﬀerentiation in the last term of two preceding equation andsubtracting,





∂
∂Q
∂δx
∂
, δ = (δQ) − δ=· ∇Q0 .
∂t x
∂t
∂t
∂t x

(B.32)

Equations (B.29) and (B.32) show that the Lagrangian perturbation commutes with neither thegradient operator nor the partial derivative with respect to time at constant position.

Appendix C

Alternate Forms of the Wave Equation

In §2.1.4, several approximations were made to arrive at the ﬁnal form of the wave
equation (2.30). For this appendix, we relax two of these assumptions, constant gravity andadiabatic stratiﬁcation, and derive the resulting wave equations.
The starting point is equation (2.22), the perturbed momentum equation under Cowling’s
approximation:∂ 2 δx∇ρ= ∇(c2 χ − gξ) + (c2+ gr̂)χ.
2
∂t
ρ

(C.1)

Our previous assumptions of constant gravity allowed us to write expressions such as ∇2 (gξ) =g∇2 ξ. In the general case, we cannot make such a statement. Instead we have four parameters(the three components of δx and χ) and thus need four separate equations to reduce the equation
to one in terms of χ only. Of course, since χ = ∇ · δx the four variables are interrelated. We
proceed by applying four operators to (C.1), yielding the following equations:(a) ∇·∇2 (gξ) = ∇2 (c2 χ) −∂2χ2βχ
+ r̂ · ∇(βχ) +
,
∂t2
r

(C.2)

(b) r̂ · ∇ × ∇×
∂2 2
2βχ
∂2
2
∇
ξ
=
∇
(βχ)
+
+
(r̂ · ∇χ),
h
∂t2
r2
∂t2

(C.3)

(c) ∇h ·(horizontal component)
∂2
∂t2



∂ξ 2ξ
+
∂r
r


− ∇2h (gξ) =∂2χ− ∇2h (c2 χ),
∂t2

(C.4)

111
(d) r̂·∂2ξ
+ r̂ · ∇(gξ) = r̂ · ∇(c2 χ) + βχ.
∂t2

(C.5)

Equations (C.2) and (C.3) are the same as equations (2.24) and (2.25), respectively.
From our work in §2.1.3, we know that the angular parts of the perturbation variables
are given by spherical harmonics. Hence, we are free to make the substitution∇2h = −

( + 1)
.
r2

(C.6)

With this simpliﬁcation, equations (C.2)–(C.5) are four linear equations for the four unknowns
˙ ξ̈, and χ where a dot represents a derivative with respect to the radial coordinate r. Simple
ξ, ξ,substitution quickly reduces the equations to the form


2ġ
g̈ −
r




∂2
2( + 1)
∂2
−
g
ġ
ξ
=
B1 (χ) − gB2 (χ) − 2ġB3 (χ)
∂t2
r2
∂t2

(C.7)

and





2g ∂ 2( + 1) 2
∂2
∂4
+
ġ
−
−
g
ξ = −gB3 (χ) + 2 B4 (χ),
4
2
2
∂t
r ∂t
r
∂t

(C.8)

where B1 (χ), B2 (χ), B3 (χ), and B4 (χ) are the right-hand sides of equations (C.2)–(C.5), respectively.
In order to eliminate the partial derivatives with respect to time, we work with only one∂
Fourier component of the solution and hence make the substitution ∂t→ −iω. Although it

is possible to produce a ﬁnal partial diﬀerential equation without any further restrictions, thealgebra is quite involved. This is why the simplifying assumptions of §2.1.4 were made: theyilluminated the underlying physics without clouding the result with mathematical complications.
In that spirit we consider a looser set of assumptions, with a better physical rationale, thanthose considered in the main text. From the spherical symmetry of the reference state, we
already know that g = g(r). However, the solar mass distribution is centrally condensed —the convection zone accounts for only 2% of the total mass despite occupying the outer 30% byradius. Furthermore, p modes of large spherical degree traverse only the outermost regions of the112convection zone. Thus, gravity for the p modes under consideration is well approximated by afree-space potential. Under this assumption, there is a simple relation between the gravitational
acceleration and its radial derivatives,
ġ =

−2g
r

and

g̈ =

6g
.
r2

(C.9)

We ﬁrst consider a system with an adiabatic stratiﬁcation (in other words, β = 0), but
gravity as speciﬁed in (C.9). Equations (C.7) and (C.8) can be reduced to the following form:


 2

∂ χ
4g ∂ 2
∂4
2 2
2 2
−
+
g
∇
−
∇
(c
χ)
+
gr̂
·
∇χ
−
h
∂t4
r ∂t2
∂t2

 
 2

4g ∂ 2
6g 2
4g 2 2∂ χ
10g ∂ 22 2
− 2
− ∇h (c χ) +−∇ r̂ · ∇(c2 χ) = 0, (C.10)
r ∂t2
r
∂t2
r2 ∂t2r h∂where equation (C.6) and the correspondence between ∂tand −iω have been used to put theresult in a form similar to (2.26). To compare to our previous results, we take the plane-parallel
limit (r → ∞). Doing so yields
 2

∂ χ
2 2
−
∇
(c
χ)
+
gr̂
·
∇χ
= 0,
∂t2

(C.11)

which is equivalent to equation (2.26) in the limit of adiabatic stratiﬁcation.
Although equation (C.10) contains odd derivatives of χ, we cannot transform it into a
Helmholtz-like form, as we did with equation (2.27), due to the ﬁnal terms on the left-hand side.
Deﬁning the diﬀerential operators
L1 =∂4
4g ∂ 2−+ g 2 ∇2h ,
∂t4
r ∂t2

L2 =

4g ∂ 2
6g 2
−
, and
r ∂t2
r2

L3 =

10g ∂ 2
4g 2 2−∇ , (C.12)
r2 ∂t2r h

and making the same substitution as before, Ψ = c2 ρ1/2 χ, we ﬁnd that

L1


 2



∂2Ψ
∂ Ψ
gΨ
2 2
2
2 2
2
−
c
∇
Ψ
+
ω
Ψ
−
L
−
c
∇
Ψ
+
L
c
r̂
·
∇Ψ
+
= 0,
2
3
c
h
∂t2
∂t2
2

(C.13)

where the acoustic cutoﬀ frequency ωc is the same as deﬁned in equation (2.28). Again thisform is equivalent to equation (2.30) in the plane-parallel limit. If desired this equation can beplaced in the form of a dispersion relation analogous to equation (2.35), although the expression
is somewhat complicated.

113
Finally, for completeness, we present the governing equation for the more general casewhere β = 0 and where derivatives of the gravitational acceleration are of the form given inequation (C.9). After some algebra, and in terms of Ψ, the result is

L1

 4



 
2βχ ∂ 2 Ψ
2βgc2
∂ Ψ
2 2
2
2 2 2
−
c
∇
−
ω
−
−
c
N
∇
−
Ψ −
c
h
∂t4
r
∂t2
r2
 4

 2 


2
gΨ
∂2Ψ
∂ Ψ
∂
2 ∂
2
2
L2
−c
∇ Ψ + L3c r̂ · ∇Ψ ++ β 2 = 0. (C.14)
∂t4
∂t2 h
∂t2
2
∂t

The ﬁrst term in square brackets is identical with equation (2.27) and in the appropriate limits
this form reduces to our previous results.

Appendix D

The Computer Code

Originally written in IDL during 1995 to implement the method described in §2.4.4, thecode has undergone several revisions since then. In 1997 it was translated into FORTRAN90
when it was realized that the beneﬁts of IDL’s graphical capabilities did not outweigh FORTRAN’s greater speed and numerical precision. We ﬁrst used adiabatic switching in a fashion
similar to the current version in 1998. For structural reasons, the code has been broken intomodules, each of which is described in the table below.

115

Table D.1: Modules used in our adiabatic switching code. Files above the line are necessary
while those below the line generate test data, but are not required for use.

File (*.f90)

Functionconstants

Deﬁnes useful constants such as the polytropic index, gravitational acceleration, domain size, and the switching period.
The top level of the integration routines from Shampine & Gordon (1975).
This routine allocates storage for variables and calls subsequent programs.
Deﬁnes and evaluates the diﬀerential equations used by the ray-tracing routine— in other words, the right-hand-side of Hamilton’s equations.
The top-level program. It sets the initial conditions, writes useful data to aﬁle for later reference, and calls the integration routines.
Part of the integration routine. It supervises the integration and evaluates its
success or failure.
The heart of the integration routine, it uses the Adams PECE formulas to solve
Hamilton’s equations. The code adjusts its order and step size to control theerror.
The ﬁnal part of the integration routine. After previous parts of the code
approximate the solution with a polynomial, this subroutine ﬁnds the solution
at the desired point.
On the ﬁrst call it loads the mesh of sound speeds into memory for quicker
access. On all subsequent calls it calculates the sound speed and directional
derivatives at a point in the domain.
Deﬁnes and evaluates the form of the switching function.
On the ﬁrst call it loads the mesh of advective velocities into memory forquicker access. On all subsequent calls it calculates the velocity and directional
derivatives at a point in the domain.
Deﬁnes and evaluates the functional form of the sound speed when creating agrid of values.
Creates a grid with a given spacing (set in constants.f90) of sound speeds.
Creates a grid with a given spacing (set in constants.f90) of advective velocities.
Deﬁnes and evaluates the functional form of the advective velocity when creating a grid of values.

eureka
func
main
mod1
mod2

mod3

sound

switchingvelocityc funcmake cmake vv func