Title: Dispersion diagrams of waves in the simulated solar chromosphere
Authors and affiliations: Dove, C., Zita, E.J., Bogdan, T.J.,
What role do MHD waves play in the dynamics of the solar chromosphere? Observation, theory
and simulation address this question in distinctly different ways, so it is helpful to draw points of
contact between them. Simulation gives us a relatively omniscient view of the system we model,
but it is only as good as the theory upon which it depends. And the only way to absolutely test
theory is by observation. By generating dispersion diagrams from simulation data, we can view
the data in a form that is easier to compare directly with simulation and theory. Further, we can
bring the data analysis techniques of observers to bear on our simulation, possibly revealing
dynamics that have eluded prior analyses.
Solar physicists expect that magnetic fields contribute to chromospheric and coronal heating.
Candidates for heating include flares, currents, and MHD waves. Where field lines or flux tubes
are twisted, bent, or moved by their embedding plasma, local field nulls can become sites of
reconnection and heating. Observations of solar flares show that when oppositely-directed field
lines cross, magnetic reconnection can drive energy releases such as coronal mass ejections.
Non-potential force-free plasmas carry field-aligned currents, though solar plasmas require
enhanced resistivity (e.g. from turbulence) for currents to provide significant ohmic heating.
While acoustic waves dissipate with altitude, MHD waves may be able to transport energy up
through the chromosphere. Clearly magnetic fields are capable of heating the solar atmosphere –
but exactly how, and to what extent?
To understand coronal heating, we must understand the dynamics of MHD waves in the solar
chromosphere. Solar physicists can take three complementary approaches toward understanding
solar magnetohydrodynamics: analysis, numerical modeling, and observation. Analytically, we
can use the MHD equations to derive wave equations and dispersion relations for waves in a
solar plasma, as described in Section 3 below. However, nonlinear relationships render many
MHD systems analytically intractable. For example, in regions where the plasma gas pressure is
comparable to the magnetic pressure, waves are not pure but mixed, as discussed in Section 3
below. In such cases, we can apply numerical methods to approximate wave dynamics.
Computer simulations can offer insights into MHD wave dynamics that are difficult or
impossible to obtain through purely analytical methods. Finally, since simulations have
limitations, the ultimate test of MHD models is direct observation. This is a two-way
relationship: at their best, MHD models can also contribute to richer understanding of
Observers and modelers generally analyze different data in different ways. A theorist would like
to be able to solve the MHD equations and obtain a dispersion relation for each wave solution, or
Page 1 of 19
a relationship between frequency omega and wavenumber k. If we plot such dispersion relations
for as omega vs. k, we may see structures that correspond to characteristic waves. The observer
measuring, say, Doppler velocity and intensities, may not have the luxury of dispersion relations,
due to noise and limitations of line-of-sight observations. Two of our goals in this paper are
(1) to generate and analyze dispersion diagrams from a computer simulation (described in
Section 3 below) in order to better understand chromospheric MHD waves and
(2) to note points of contact between our simulation analyses and observational techniques, in
order to deepen understanding of these complementary analyses.
We start by describing, in Section 3, the MHD model that generates our numerical data. We
summarize earlier analyses of the same datasets by Bogdan et al. (2001, 2003). Then we
describe our methods for generating and analyzing dispersion diagrams from the simulation data.
We discuss the results of our analyses and their relation to a few observations. Finally, we
outline future work.
3) Summary of Calculations:
Our data are generated by a 2D MHD simulation code written by Nordlund & Galsgaard (1995).
The code self-consistently evolves density rho, pressure p, magnetic field B, and velocity in a
stratified (del*p= - rho*g) isothermal atmosphere, by solving the non-linear equations of ideal
While the code is capable of runs in 3D, we concentrate on data from a run limited to 2 spatial
dimensions, z and x where z is the height above the photosphere and x is the horizontal distance
along the photosphere. In the z-direction there are 294 height steps at 4.33km per step for a total
of 1.27 Mm. In the x-direction there are 500 horizontal steps at 15.8km per step for a total of
7.9 Mm. The code is evolved in time at 1.3 s per timestep with 161 steps for a total duration of
The initial magnetic field configuration consists of a large, central unipolar field flanked by two
smaller concentrations of opposite polarity field. The central concentration is 2.0 Mm wide at
the photosphere, while the two satellites are each 750 km wide. This configuration imitates a
Page 2 of 19
magnetic region of the photosphere scaled down by a factor of 10 (Johnson, Petty-Powell 2002).
Flux is concentrated near the driver in the photosphere and spreads out to form a “magnetic
canopy.” The strong magnetic field simulates a network region of the solar chromosphere (e.g.
near a sunspot), with penumbral “wings” of return flux.
Our isothermal atmosphere at a temperature of 5785 K, has constant sound speed of 8.49 km/s.
Since the density and magnetic field change independently in space, the Alfvén speed
v A = B 2 4πρ is not constant and generally increases rapidly with height (Bogdan et al., 2002).
Our driver is a piston 0.4 Mm wide located slightly left-of-center in the simulation. The piston
oscillates in the z-direction with a period of 29.3 s for a frequency of 43.2 mHz, about 10 times
faster than the familiar 5-min p-mode oscillations (Leighton 1962).
We are interested in how MHD waves are generated at the photosphere, how they travel through
the chromosphere, and how they transform into other waves, or “mix” into other “modes”. MHD
waves speeds are generally a combination of two characteristic speeds: the sound speed
cs = γ P ρ and the Alfvén speed v A = B 2 4πρ , where ρ = mass density, P = gas pressure, γ
= ratio of specific heats, and B = magnetic field strength. Plasma magnetohydrodynamics are
generally controlled by the local ratio of gas pressure to magnetic pressure, β ≡ 8π P B 2 , and
magnetic fields and plasma are “frozen” together by the medium’s high conductivity. Where β
>>1, the gas pressure exceeds the magnetic pressure, and moving plasma drags its embedded
field along. Since β = 2cs2 γ vA , acoustic waves are “fast” and Alfvénic waves are “slow” in
these “high-beta” regions. Where β <<1, the magnetic pressure dominates, and moving field
lines drag their embedded plasma along. In these “low-beta” regions, Alfvénic waves are fast
and acoustic waves are slow. In our simulation, the driver excites acoustic waves at the base of
the computational domain. These pressure modes, or p-modes, propagate upward along the
spreading magnetic field lines. The photospheric driver also can excite different types of
magnetic waves. Magnetosonic waves propagate longitudinally across field lines, compressing
field lines and their plasma. Alfvén waves propagate transversely along field lines and can be
thought of as disturbances excited by the “plucking” of magnetic field lines.
Those are the simpler cases.
We are interested in more complex cases. Mode mixing is especially rich where the two
characteristic speeds are comparable, that is, where beta approx 1. This is the case near the
“magnetic canopy” at the boundary of strongly magnetic “network” regions; note the bold
Page 3 of 19
beta=1 line in Fig.1.
Figure 1 : Fluctuations in fractional density deltaRho/Rho_0 at an elapsed time of 19.5 s. Z is height above the
photosphere in Mm. X is the position along the surface of the photosphere in Mm. Light, mostly vertical lines are
initial magnetic field. Numbered lines are plasma beta contours.
In such regions, MHD waves are generally not pure modes, but hybrids of Alfvénic, acoustic,
and magnetosonic waves, depending on the local gas pressure and field strength, and on the
orientation of oscillations with respect to the local magnetic field.
Page 4 of 19
Magnetic and acoustic waves can perturb the plasma properties and therefore the nature of waves
traveling through them. For example, magnetic compressions due to magnetosonic modes
periodically decrease local beta, and gas pressure increases due to acoustic waves periodically
increase local beta. Magnetic perturbations can yield plasma oscillations, and vice versa. The
resulting nonlinear wave dynamics can be modeled by a good MHD simulation, such as the
model of Nordlund et al. described above. By comparing dispersion plots of simulation data
with actual observations of the chromosphere, we can investigate how MHD waves change
character as they propagate, and how the waves transport energy into the solar magneto-
We aim to reveal new MHD wave dynamics by generating dispersion plots from 2D simulation
data described in Section 3. We hope that our plots will also provide a point of contact between
our simulation and observation. We wrote Interactive Data Language (IDL) codes to compute
the 2D Fast Fourier Transform of our data, transforming a signal of amplitude varying with x and
t to a power spectrum varying with frequency and wavenumber. We also used averaging
techniques and a Hanning data window to reduce spurious peaks in our spectra.
Finally, we used only a small portion of our total FFT data for a given plot, limiting our analyses
to a domain where MHD is valid. It is only valid for low frequencies and wavenumbers because
of two central assumptions:
1) Plasma is a fluid, so we ignore its kinetic properties, assuming that collisions do not
contribute significantly to plasma dynamics.
2) Linearization: we look at small perturbations from equilibrium quantities and throw out
resultant higher-order terms in our equations of motion and continuity. Because the
perturbations are small, we can assume that they are adiabatic, and we can assume that
the electrons and ions move together in the plasma. Therefore, we disregard ion/electron
Page 5 of 19
To obtain a dispersion plot for the wave objects in our simulation, we need to transform
amplitudes (e.g. of density or speed) varying in space and time to powers in frequency and
wavenumber. The Fourier transform converts spatial signals to wavenumber spectra and time-
varying signals to frequency spectra. For a continuous function, the Fourier transform H(f) of a
signal h(t) can be found from
∫ h(t )e
H( f ) = dt
∫ H ( f )e − 2πift df
h(t ) =
This integral transform depends on the assumption that any periodic function can be represented
as a sum of sinusoids, or a Fourier series (Brigham 1988). However, our data are not continuous.
To achieve analogous results to the integral Fourier transform we need to use a special
application for the discrete data. The Discrete Fourier Transform (DFT) G of a discrete signal g
can be found from
Gn ≡ ∑ g k e 2πikn / N Eq. 4
Now we can transform finite datasets with discrete steps in space or time to finite datasets with
discrete steps in frequency and wavenumber. While calculating a Discrete Fourier Transform is
computationally intensive, Tukey and Cooley developed an algorithm in 1965 to efficiently
compute the DFT. The Fast Fourier Transform (FFT) algorithm decreases the computation time
of the DFT dramatically for datasets with 2^n elements (Brigham 1986). We use the FFT
algorithm included with IDL (Interactive Data Language), which is our primary tool for the
analyses below. Raw FFT output in IDL gives data incremented by frequency and inverse
wavelengths with a spacing of 1/N where N is the number of discrete steps in the given data
series (IDL Online Help). To convert the inverse wavelengths to actual wavenumbers, we
multiply each inverse wavelength element by 2pi/Deltax where Deltax is the extent of each step
in x. To convert FFT frequencies to actual frequencies, we divide each frequency by 1/Deltat
where Deltat is the extent of each step in t.
Numerical data from the 2D MHD code are stored in a 3D array or “datacube”
approximately 750MB in size. It contains elements in z, x, and t that represent rho, p, B, and u
(Nordlund & Galsgaard 1995). We extract these data from the datacube using an IDL routine
authored by Bogdan et al. 2001. This routine returns a 2D "slab" of data in one spatial
dimension vs. time. For the bulk of our analyses, we chose slabs of data in x vs. t at a given z.
This corresponds to the horizontal position along the surface of the photosphere with respect to
Page 6 of 19
time at a given height. We then extracted slabs for several different variables, including
fractional density Deltarho/rho_0, plasma velocity perpendicular to the magnetic field u_perp,
horizontal plasma velocity u_z, and plasma velocity parallel to the magnetic field u_parallel
We wrote IDL programs using a 2D FFT function to generate plots of frequencies versus
wavenumbers from our 2D slabs of data, as in Fig 2.
Figure 2: Frequency f vs. wavenumber k plot of fractional density deltaRho/Rho_0 fluctuations in x and in t at an
altitude of 0.2165 Mm. Power, on the vertical axis, is scaled by the natural log to reveal small-scale features.
Page 7 of 19
Then we used a combination of methods to get clearer information from the raw output of these
FFTs. The raw FFTs present challenges for several reasons. Recall that an FFT generates a
representation of space/time data in frequency/wavenumber space for *periodic* data. When
non-periodic data are transformed, or the signal contains a strong frequency component that falls
between frequency bins in the FFT dataset, some spurious peaks result, as evident in Fig 3b.
This phenomenon is called "leakage" and manifests itself in spectral plots as small, regularly
spaced peaks surrounding a real peak at a given wavenumber or frequency. These
discontinuities can be caused by the discrete representation of periodic data whose periods do not
fit perfectly into the given data window, as illustrated in Fig 4.
Figure 4: Sin(x) on the vertical axis and x on the horizontal axis. Gray box shows how a rectangular data window
can introduce a discontinuity to discretely sampled data.
Thus these peaks "leak" into neighboring frequencies or wavenumbers (Press 1986). We
circumvent this problem by applying a window function to the input data. A window function
tapers a dataset to eliminate spurious edge discontinuities. We choose a Hanning window (Fig
Page 8 of 19
Figure 4a: Hanning window in two dimensions, used to taper edges of a 2D dataset, eliminating discontinuities (see
This yields omega-k diagrams with the spurious edge peaks reduced (Fig 5b).
Figure 4b: Frequency f vs. wavenumber k plot of Hanning windowed fractional density deltaRho/Rho_0 fluctuations
in x and in t at an altitude of 0.2165 Mm. Power, on the vertical axis, is scaled by the natural log to reveal small-
scale features. Note smoother peaks and clearer features.
Page 9 of 19
Another limitation of the FFT is "aliasing". Recall that the highest frequency signal
which can be resolved in data spaced by time intervals dt is f_N = 1/2dt, and it is defined as the
Nyquist critical frequency (Press 1986). Similarly, the shortest wavelength resolvable in data
spaced by intervals of dx yields the largest wavenumber, k_max = 2 pi / lambda_min = 2 pi / dx.
When the data contain frequencies and wavenumbers greater than the Nyquist critical
frequency/wavenumber, these frequencies and wavenumbers are spuriously shifted below the
Nyquist limit (Fig. 6).
Figure 5: Right: Sin(x) + Sin(2x) + Sin(5x) sampled with 32 points, giving a sampling frequency of 5.1
samples/2Pi. The Nyquist critical frequency corresponds to a sinusoid that completes a total of 16 cycles, giving
exactly 2 samples per cycle. Left: FFT of Sin(x) + Sin(2x) + Sin(5x). Note that there should be 3 clear peaks, but
the peak corresponding to Sin(5x) would lie above the Nyquist critical frequency, so it is spuriously shifted under
the Nyquist limit, which is the 16th frequency bin.
One way to circumvent this problem is to analyze bandwidth-limited data, i.e. data which contain
no high frequency or wavenumber components (Press 1986). Another way is to increase the
sampled resolution of the discrete data to include all possible frequencies (Press 1986). In our
case, we concentrate our analyses on signals with low frequency and wavenumber. Since our
driving frequency of 43.2 mHz is much smaller than the Nyquist frequency (104.5 Hz), we
expect little aliasing in frequency. Similarly, our driver, with a spatial extent of dx=0.4 Mm, has
a wavenumber k=2 pi / dx of 15.7 Mm^-1. Since the driver wavenumber is much smaller than
the Nyquist wavenumber (1570 Mm^-1), we expect little aliasing in wavenumber. We apply an
averaging algorithm (Fig. 7a) to remove the general contribution of small "wiggles" in adjacent
slabs of data. This yields smoother omega-k diagrams, as in Fig 7b, which emphasize strong,
For each altitude z above the photosphere, we analyze a set of "slabs", or data varying
with x and t. We apply the methods described above to a total of 11 slabs near a given altitude: 5
slabs above the given height, 5 below, and 1 slab at the height itself.
Then we compute the 2D FFT of each slab and average together the results, spanning a thickness
of 17.38 km in altitude. This numerical averaging is analogous to intrinsic observational
Page 10 of 19
averaging due to the finite spread in altitude of a signal of given wavelength. For example, the
91 A UV line samples plasma from z=__ to __ in the chromosphere (Judge et al 2002)
Combining the averaging technique with the Hanning window function, we are able to smooth
the numerical data substantially. By limiting our analyses to a small "corner" of our 2D FFT
data, f : 0 - 60 mHz and k : 0 - 80 Mm^-1, and plotting the log of the power, we eliminate high-
frequency and wavenumber peaks that may be aliased.
We use three types of plots to visualize our dispersion relations: contour plots, scaled
images, and 3D surfaces, as in Fig 8.
Figure 8: Averaged frequency f vs. wavenumber k plot of fractional density deltaRho/Rho_0 fluctuations in x and in
t at an altitude of 0.2165 Mm. 11 FFT’s were averaged around the given height, corresponding to a total spread of
47.6 km in altitude. Power, on the vertical axis, is scaled by the natural log to reveal small-scale features. Top:
contour plot, middle: surface plot, bottom: grayscale image
Contour plots (Fig. 12) are useful for estimating slopes and thus speeds. Recall that slopes (d
omega/dk) in omega-k space represent group velocities or the speeds at which waves propagate
in the x-direction in our simulation. We plot f vs. k instead of omega vs. k, in order to directly
Page 11 of 19
see the effects of the simulation and its harmonics, so we can convert slopes to group velocities
by applying a scale factor of 2*pi. Scaled images (Fig 8, bottom) and 3D surfaces (Fig 8,
middle) are good for spotting overall patterns.
We will first describe fundamental features common to our dispersion diagrams, and then discuss
in more detail the patterns evident in different variables at different heights. Common features
include horizontal bands centered at various frequencies, vertical bands centered at various
wavenumbers, and a plaid pattern due to their intersections; adjacent stripes with similar slopes,
and s-curves of varying slopes. Analyzing variables at different heights will reveal wave
transformations as they rise through the chromosphere.
I) Common features: One dominant feature is the horizontal band of power at a frequency of
~40 mHz , and weaker bands at harmonics. For example, in Fig 8(bottom) and Fig 9(top), the
second band is clearly visible around 80 mHz. This feature is visible in plots of rho, u_perp, and
u_z, and becomes fainter with increasing height. At ~0.21 Mm, the feature is a clear band in all
plots, but at a height of ~0.8 Mm the power has spread to neighboring frequencies and drops off
exponentially with increasing wavenumber. The ~40 mHz band and harmonics are due to our
driver, with its fundamental frequency of 43.2 mHz and spatial extent of 0.4 Mm at the
photosphere. The relatively high amplitude and sharply defined boundaries are consistent with
wavefronts excited by this radial oscillator.
Similarly, there are evenly spaced vertical bands in wavenumber, and their spacing
changes with height (Fig 9).
Page 12 of 19
Frequency vs. wavenumber for fractional density Frequency vs. wavenumber for vertical velocity u_z at z =
deltaRho/rho_o at z = 0.2165 Mm 0.2165 Mm
Frequency vs. wavenumber for fractional density Frequency vs. wavenumber for vertical velocity u_z at z =
deltaRho/rho_o at z = 1.0825 Mm 1.0825 Mm
Figure 9: k-spacing of vertical bands decreases with height in dispersion plots of fractional density deltaRho/rho_o
and vertical velocity u_z.
These bands are clearest in rho and in u_z. Near the photosphere, the bands are clear and tightly-
spaced by Del k~13 Mm^-1. This corresponds to a wavelength of 0.4 km, which is the width of
the driver. We see evenly-spaced harmonics in wavenumber as well. The spacing between
wavenumber bands decreases with height, consistent with a spreading spatial extent of the driver
as p-modes are channeled up the spreading magnetic structure (Fig 1). At a height of 0.8 Mm the
wavenumber spacing is much smaller, at 3 Mm^-1, corresponding to a larger wavelength of 2.0
Mm. We can quantify the increasing spread of the driver's effect with height:
(x = xd e z/Hd)
where we find the characteristic height for the spread of the rising p-modes Hd = 0.7 Mm (Fig
Page 13 of 19
Effective driver width (Mm) varies with height (Mm)
y = 0.407e
1 R = 0.9649
Width of driver in Mm (y)
Effective width of driver
Expon. (Effective width of
0 0.2 0.4 0.6 0.8
Height in Mm (x )
Figure 11: Effective driver width varies with height. 0.4 Mm is the actual driver width and 0.7
Mm (1/1.4) is the driver scale-height.
This scale driver scale-height depends on both the particular shape of the network region and on
the pressure scale height H_p=_(equation)_ (CBD: I don’t understand why pressure comes in
here). The increasing spread of the network region can be similarly characterized by
x = xb e z/Hb
where the characteristic height for the spread of the magnetic canopy Hb = 0.7 Mm . These scale
heights are unique to our simulation (where H_p=__#_) and may vary in real network regions
with various shapes and magnetic field strengths.
The grid-like or plaid pattern of intersecting horizontal and vertical bands in the
dispersion diagrams is therefore due to the harmonics in time and space of the driving piston at
the photosphere, as acoustic modes rise into the chromosphere and spread across their magnetic
Next consider the regions of adjacent stripes, all with similar slopes. By calculating the
slopes, or group speeds, of such features, we gain insight into the waves present at each height in
the simulation. For example, in Fig 12 we see slopes corresponding both to weak fast modes
(slopes a, b and c) and slow modes (slope d).
Page 14 of 19
Figure 12: Frequency vs. wavenumber for fractional density deltaRho/rho_o at height z = 0.65 Mm. Labeled slopes
correspond to group speeds. a = 27 km/s, b = 36 km/s, c = 74 km/s and d = 9 km/s.
II) Changes with height:
Another interesting feature in our dispersion plots is an “s-curve". At a given altitude z, the
slope in f-k appears to change smoothly. These slopes also vary with altitude in the
chromosphere. At heights near z ~ 0.4 Mm, we see vertical, repeating s-shaped curves with the
same slopes above and below the driving frequency, and with slopes approaching infinity near
the driving frequency (Fig 12a). Their spacing in k is consistent with the spreading driver effect
discussed in the paragraph above. At z ~ 0.65 Mm, the s-curve slopes increase with
wavenumber below the driving frequency (Fig 12b) while curves above the driving frequency
have slopes of roughly 9 km/s /pm dv. This is close to the sound speed cs = 8.49 km/s. Curves
below the driving frequency, labeled a, b and c respectively in Fig 13b, have supersonic slopes.
The speeds of each curve are found to be v_a = 27 km s-1, v_b = 36 km s-1, and v_c = 74 km s-1.
Supersonic waves always have a magnetosonic component, as evident from Equation 2.
We cannot readily distinguish Alfvénic from magnetosonic waves in the dispersion diagrams,
and in practice these modes are generally mixed in a real magnetic canopy region, as discussed
in Section 3 above. Since the Alfvén speed varies with position x, due to the varying magnetic
field strength, we cannot assign a unique value for this eigenmode at a given height z. We can,
however, gain insight by comparing our omega-k plots with data plotted in space (x) and time,
at the same altitudes (z) (Fig. 14a,b).
Page 15 of 19
Figure 13: Fractional density deltaRho/rho_o varies with time t and position along photosphere x at given altitudes
z = 0.43 Mm (left) and z = 0.65 Mm (right). Note that greater wavefront slopes indicate slower waves in these t vs.
Notice in Fig.13a that the dominant features are very fast waves centered over the driver.
These are flanked by slow, weak fluctuations near x = 2 Mm and x = 6 Mm whose speeds are
comparable to the sound speed of 8.49 km s-1. This may seem curious since acoustic waves
centered over the driver ought to be slow in a low-beta plasma. What we see here is a projection
effect. The dispersion diagram represents waves traveling in the x-direction, and since the
magnetic field is mostly vertical, acoustic waves travel primarily in the z-direction, especially
near the photosphere. We detect their wavefronts passing through our plane of view as strong
density fluctuations spreading out horizontally as they rise in z, therefore traveling very rapidly
in the x-direction The more planar the acoustic modes are, the faster they spread in the x-
direction at a given height. As the sound waves change shape with increasing height, channeled
by bending field lines, they become less planar.
The supersonic slopes of curves in f vs. k of density fluctuations at a given height
correspond to fast hybrid waves. This is evident when we compare plots of f vs.k and t vs. x at z
= 0.4 Mm and z = 0.65 Mm and consider the hybrid wave speeds in Eqn 2 above.
At z = 0.4 Mm, our magnetic field has a substantial horizontal component (Fig 13a).
Since we limit our analysis to waves traveling in the x-direction, we should see waves in our f vs.
k plots with hybrid speeds in the neighborhood of v+ with phi=0 (Eqn 2a). When we examine
similar plots at z = 0.65 Mm, we see that the magnetic field is mostly vertical, so we see waves
traveling perpendicular to the magnetic field. The fast hybrid mode travels fastest perpendicular
to field lines ( at v+ with phi=0 (Eqn 2a)) due to its magnetosonic component. So, at greater
heights we see more magnetosonic than Alfvén modes traveling in the x-direction, and at lower
heights we see a mix of both. This makes sense, since the magnetic field lines, along which
Alfvenic waves travel, are more nearly vertical at higher altitudes. The only waves available (in
2D) to cross vertical field lines are magnetosonic modes. At lower altitudes, e.g. below the
magnetic canopy and in the return flux wings, oblique field lines can support components of both
magnetosonic and Alvenic waves in the x-direction.
At lower heights, we see interference between waves that are mode-converting near beta
= 1. This makes sense if we refer to the initial magnetic field configuration (Fig 1) and recall
Page 16 of 19
that, at low heights, the magnetic field bends into penumbral return-flux loops, changing from
vertical to horizontal. In general, weak fast modes should get faster with increasing altitude
where the magnetic field is everywhere vertical. Comparing dispersion plots …
At greater heights, we see two distinct dark bands with negative slope, as in Fig.14 for
Figure 14: Grayscale dispersion plot of perpendicular velocity u_perp at z = 0.87 Mm. Aspect
ratio adjusted to match previous plots. Note dark right-corner and two clear, dark bands with
decreasing negative slopes.
They become apparent in our plots of both u_perp and rho at a height of z ~ 0.65 Mm. These
dark bands represent troughs between the ridges of high signal in perpendicular velocity
fluctuations. Recall that these fluctuations are primarily due to magnetic waves, as acoustic
waves channeled along magnetic field lines have parallel velocity fluctuations. What causes
Another clear feature in both variables is the suppression of short-wavelength (high-k),
low-frequency waves with increasing height. This is evident in the dark lower-right corner of
grayscale plots such as Figs. 14 and 15
Page 17 of 19
Figure 15: Dispersion plot of fractional density deltaRho/rho_o at z = 0.87 Mm. Dark lower-right corner of
grayscale plot shows suppression of short-wavelength, low-frequency waves. Also note the two dark bands with
decreasing negative slope as discussed in the paragraph above and shown in Fig. 14.
Note that this corner is brighter, that is the high-k, low-f signal is stronger, nearer to the
photosphere, as in Fig._. Two possible explanations for this could be that:
1) There are more short-wavelength features at low altitude because more slow mode
wavefronts have had time to mode-convert to magnetic waves and interfere with one
2) High-frequency ripples travel faster than low-frequency ripples, so at greater heights we
see fewer ripples.
Further analysis will be necessary to determine what is responsible for this effect.
Page 18 of 19
6) Future Work
For future work, we have further refinements to our method that we may implement. To aid the
comparison of our data with observational data, we may weight our averaging algorithm (see
Section 4 above) with a Gaussian curve and adjust its span in altitude to more closely
approximate the signal spread that observers see due to Doppler broadening of the signal (Judge
When we obtain output from our FFT, we focus on the quadrant of the f-k diagram near
the zero-order frequencies and wavenumbers. However, focusing here only reveals information
about left-traveling waves. The information about right-traveling waves appears near the
maximum-order frequency and zero-order wavenumber, or vice-versa. Because our driver is
offset to the left (Fig 1), our data is asymmetrical. So there may be new right-traveling waves
that have no left-traveling analogues, and these could be revealed by performing the same
analyses on FFT output in the area mentioned above.
As our MHD model is refined, we will gain access to new simulation data with, for
example, better boundary conditions that limit reflections along the upper computational
boundary. This will allow us to examine longer timeseries and greater heights, such as the
“ghost zone” region above 1 Mm (Bogdan et al 2003).
We will also examine smaller regions in space and/or in time by “zooming” in on areas of
interest in our slabs of data. Currently, the resolution of our data limits renders this impractical.
But the MHD model actually outputs data at a far finer resolution than we have employed in our
analyses. Our choice of 15.8km steps in height and 1.3s steps in time is purely practical and we
are not limited to these dimensions.
Neither are we limited to our choice of examining slabs of data varying with x and with t.
We can also look at slices of data varying in z and in t, which would reveal far more information
about waves traveling in the z-direction than do our slabs of data.
Page 19 of 19