simulation of magnetohydrodynamic shock wave ... - IOPscience

The Astrophysical Journal, 708:268–287, 2010 January 1  C 2010.

doi:10.1088/0004-637X/708/1/268

The American Astronomical Society. All rights reserved. Printed in the U.S.A.

SIMULATION OF MAGNETOHYDRODYNAMIC SHOCK WAVE GENERATION, PROPAGATION, AND HEATING IN THE PHOTOSPHERE AND CHROMOSPHERE USING A COMPLETE ELECTRICAL CONDUCTIVITY TENSOR Michael L. Goodman1 and Farzad Kazeminezhad Advanced Technologies Group, West Virginia High Technology Consortium Foundation, 1000 Galliher Drive, Fairmont, WV 26554, USA; [email protected], [email protected] Received 2009 February 6; accepted 2009 November 9; published 2009 December 10

ABSTRACT An electrical conductivity tensor is used in a 1.5D magnetohydrodynamic (MHD) simulation to describe how MHD shock waves may form, propagate, and heat the photosphere and chromosphere by compression and resistive dissipation. The spatial resolution is 1 km. A train of six shock waves is generated by a sinusoidal magnetic field driver in the photosphere with a period T = 30 s, mean of 500 G, and variation of 250 G. The duration of the simulation is 200 s. Waves generated in the photosphere evolve into shock waves at a height z ∼ 375 km above the photosphere. The transition of the atmosphere from weakly to strongly magnetized with increasing height causes the Pedersen resistivity ηP to increase to ∼2000 times the Spitzer resistivity. This transition occurs over a height range of a few hundred kilometers near the temperature minimum of the initial state at z ∼ 500 km. The initial state is a model atmosphere derived by Fontenla et al., plus a background magnetic field. The increase in ηP is associated with an increase in the resistive heating rate Q. Shock layer thicknesses are ∼10–20 km. They are nonzero due to the presence of resistive dissipation, so magnetization-induced resistivity plays a role in determining shock structure, and hence the compressive heating rate Qc . At t = 200 s the solution has the following properties. Within shock layers, Qmaximum ∼ 1.4–7 erg cm−3 s−1 , and Qc,maximum ∼ 10–103 Qmaximum . Between shock waves, and at some points within shock layers, Qc < 0, indicating cooling by rarefaction. The integrals of Q and Qc over the shock wave train are F ∼ 4.6 × 106 erg cm−2 s−1 and Fc ∼ 1.24 × 109 erg cm−2 s−1 . A method based on the thermal, mechanical, and electromagnetic energy conservation equations is presented for checking the accuracy of the numerical solution, and gaining insight into energy flow and transformation. The method can be applied to higher dimensional simulations. It is suggested that observations be performed to map out the transition region across which the transition from weakly ionized, weakly magnetized plasma to weakly ionized, strongly magnetized plasma occurs, and to correlate it with net radiative loss. Key words: MHD – shock waves – stars: chromospheres – Sun: chromosphere – Sun: magnetic fields Online-only material: color figures

the electromagnetic, thermal, and center-of-mass (CM) kinetic energy reservoirs in these regions. Except near the temperature minimum, where most of the free electrons come from singly ionized heavy ions, protons are the dominant ion species (FAL), and Mi ∼ Mp , where Mp is the proton magnetization. The chromosphere is commonly defined as the layer of gas between the photosphere and corona in which H is mostly neutral, and that has a lower boundary at the height above the photosphere at which the gas temperature first begins to increase with height. A different definition is proposed here. Chromospheric plasma is defined here as H plasma that is weakly ionized and strongly magnetized. Chromospheric heating may occur in shock waves or in the absence of shock waves. Previous modeling of chromospheric heating in the absence of shock waves leads to the conclusion that emission identified with the chromosphere does not exist if the plasma is not strongly magnetized and weakly ionized, but does exist if the plasma has these properties, and if one or more sufficiently strong electric current generating processes are present in the plasma (Goodman 2000, 2001, 2004a, 2004b, 2008; Kazeminezhad & Goodman 2006, henceforth KG06). These processes drive proton Pedersen currents that heat the plasma through proton–H i collisions, generating observed chromospheric emission. Electron and heavy ion current dissipation make a relatively small contribution to this heating, and do so only near the base of the

1. INTRODUCTION For collision dominated plasmas, such as the photosphere and chromosphere, the magnetization of a charged particle is the product of its cyclotron frequency and total collision time, which is the inverse of the sum of its collision frequencies with neutral and unlike charged particles. Consider a three-species plasma of electrons, singly charged ions, and neutral atoms. Let Me and Mi be the electron and ion magnetizations. For the purpose of this paper, define the plasma to be weakly (strongly) magnetized if Me Mi  1( 1). The reason the product of the magnetizations is used in this definition is that it determines the effective resistivity of the plasma, described in Section 2. Let ρ and ρn be the total and neutral mass densities. Define the plasma to be weakly (strongly) ionized if ρn /ρ ∼ 1( 1). For the temperatures and densities characteristic of the solar atmosphere (Vernazza et al. 1981; Fontenla et al. 2002, henceforth FAL; Avrett & Loeser 2008), and for measured or estimated magnetic field strengths in the atmosphere, the photosphere is weakly ionized and weakly magnetized, the chromosphere is weakly ionized and strongly magnetized, and the corona is strongly ionized and strongly magnetized. The large differences in the degrees of ionization and magnetization between these regions imply large differences in how energy flows between 1

Author to whom any correspondence should be addressed.

268

No. 1, 2010

MHD SHOCK WAVE HEATING OF THE CHROMOSPHERE

chromosphere where Me Mi  1. The results of the magnetohydrodynamic (MHD) shock wave modeling presented here suggest that magnetization-enhanced resistivity is important in determining shock structure, and that the compressive heating rate in shock layers is orders of magnitude greater than that due to Pedersen current dissipation. This result places emphasis on the long-standing question of to what degree is chromospheric heating due to processes in shock waves, and due to processes that are independent of shock waves. All of the mass and energy that drive coronal dynamics flows from the photosphere through the chromosphere into the corona. The chromosphere merges with the corona through the transition region (TR). The TR has a structure strongly determined by chromospheric dynamics. It might be created and maintained by the diffusion of H i from the interior of magnetic structures containing chromospheric plasma at a temperature T ∼ 104 K into surrounding corona plasma at T ∼ 106 K where it is collisionally excited and ionized (Judge 2008; Judge & Centeno 2008). Spicules are examples of such magnetic structures that drive mass from the chromosphere into the TR. The chromosphere and TR form a boundary layer between the photosphere and corona. Coronal models must use a lower boundary that represents the merging of chromospheric plasma into coronal plasma. Boundary conditions strongly determine solutions to models. Chromospheric dynamics is expected to strongly determine coronal dynamics. Although the approximate thickness of the chromosphere is a few thousand kilometers, its net radiative loss (NRL) of ∼107 erg cm−2 s−1 in network regions, and ∼(4–5) × 106 erg cm−2 s−1 in internetwork regions is ∼10–100 times greater than that of the corona. Only a few percent of the fraction of the total energy flux driven from the photosphere into the overlying atmosphere that is thermalized needs to reach the corona to balance its thermal energy loss by radiation and by thermal energy flow to the lower atmosphere, and to provide the thermal component of the energy needed to accelerate the solar wind. The remaining 90%–99% of this fraction balances the NRL from the chromosphere. The chromosphere is highly structured by laminar and turbulent cellular convection and magnetic fields (Cauzzi et al. 2008). Compared with the chromospheric internetwork, the chromospheric network, which outlines the boundary regions of the supergranulation, is the site of relatively high concentrations of kilogauss strength magnetic field structures. These structures are strongly correlated with enhanced emission, and hence heating (Solanki 1993; Solanki et al. 1994; Schrijver & Zwaan 2000; Akasofu & de Jager 2001; Tu et al. 2005). Although it is generally believed that MHD processes dominate heating of the chromospheric network, debate continues as to the importance of acoustic wave dissipation in heating the chromospheric internetwork (Fossum & Carlsson 2005a, 2005b, 2006; Carlsson et al. 2007; Cuntz et al. 2007). The question most recently addressed is whether acoustic waves with frequencies ∼5–40 mHz are generated with sufficient power in the convection zone to substantially heat the internetwork chromosphere, noting that acoustic waves with frequencies 50 mHz are strongly radiatively dissipated in the photosphere. Observations suggest the power in acoustic waves with frequencies of 5–40 mHz is a factor of ∼5 too small to balance the NRL from the internetwork chromosphere. Identifying the mechanisms in the solar atmosphere that generate plasma kinetic energy and the resulting radiation requires

269

using models that include sufficiently accurate descriptions of the processes that determine the flow of energy between the electromagnetic and particle kinetic energy reservoirs. One such class of models is MHD models that include effects of transport processes. These processes are determined by the local single particle distribution functions. Transport processes are represented in MHD models by tensors that are anisotropic and inhomogeneous. The anisotropy results from the preferred direction defined by the magnetic field when it is strong enough to cause the particle cyclotron frequencies to be much greater than the particle collision frequencies due to short range, binary collisions, or due to the longer range, collective interaction of many particles. These two basic types of scattering processes occur on lengthscales less than the Debye length λD , and on lengthscales  λD , respectively. Examples of transport processes are electrical and thermal conduction, and viscous and thermoelectric effects. A proposed atmospheric heating mechanism can be tested in a meaningful way using an MHD model only if the model includes a sufficiently accurate description of the relevant transport tensors, and is solved with sufficient spatial and temporal resolution in the parameter ranges of density, temperature, velocity, and magnetic field strength that characterize the region and processes of interest. The definition of “sufficient(ly)” depends on the problem. The fact that the chromosphere, TR, and corona are strongly magnetized plasmas implies that the relevant transport tensors, whether describing transport due to short-range binary, or long-range many particle interactions, are anisotropic. Modeling electrical conduction in the chromosphere requires using an anisotropic electrical conductivity tensor. Temperature gradients in the TR and in current sheets are expected to be sufficiently large to require using anisotropic tensors to describe thermoelectric effects (Goodman 1998, 2005), and effects of electron and ion thermal conduction and viscosity. Ideal MHD models do not include transport processes. These models are not useful for estimating the importance of proposed atmospheric heating mechanisms since they cannot predict resistive or viscous heating rates, or electron and ion thermal energy fluxes due to thermal conduction and thermoelectric effects. Poynting’s theorem shows that ideal MHD allows the electromagnetic and CM kinetic energy reservoirs to exchange energy, but does not allow the electromagnetic and thermal energy reservoirs to exchange energy. Here an MHD model that contains a complete electrical conductivity tensor for a three-specie plasma of electrons, singly charged ions, and neutral atoms is used to simulate a process by which driven, smooth MHD waves are generated in the photosphere, propagate upward, steepen into shock waves near the FAL temperature minimum, and continue to propagate upward and heat the atmosphere within the shock wave. The height integral of the resistive heating rate per unit volume over the simulated shock wave train yields a heating flux sufficient to balance the NRL of the internetwork chromosphere. The corresponding integral of the compressive heating rate per unit volume is ∼102 times larger. The results presented here are complementary to the onedimensional, non-LTE (NLTE) hydrodynamic simulations of Carlsson & Stein (1992, 1994, 1995, 1997, 2002) of acoustic shock wave heating of the internetwork chromosphere. Those simulations reproduce the main features of the Ca ii H and K bright points, but produce an NRL rate that is too low by at least one order of magnitude to balance the chromospheric NRL. The collective results of these simulations and of observations of

270

GOODMAN & KAZEMINEZHAD

chromospheric emission suggest that acoustic wave dissipation is not a significant source of heating for the network or internetwork chromosphere, or the TR (Carlsson et al. 1997, 2007; Judge et al. 1997, 2003; Judge & Carpenter 1998; Fossum & Carlsson 2005a, 2005b, 2006). The resistive heating mechanism that can by itself provide a heating rate sufficient to create and maintain the chromosphere is outlined in Section 1.1. This mechanism involves a magnetization-induced, several orders of magnitude increase in the Pedersen resistivity with increasing height above the photosphere. This resistivity is present whether or not shock waves form, although it is modified by conditions in a shock layer since it depends on magnetic field strength, temperature, and charged and neutral particle densities. As a diffusive process, the associated resistive dissipation plays a role in determining shock wave structure, and hence the compressive heating rate. Shock structure is determined by a balance between the ideal process of nonlinear steepening, present in ideal MHD, and diffusive or effectively diffusive processes such as resistive and viscous dissipation, thermal conduction and thermoelectric effects, and radiative cooling. 1.1. The Basic Chromospheric Resistive Heating Mechanism Define Γ ≡ (ρn /ρ)2 Me Mi . In the photosphere and chromosphere, Γ ∼ Me Mi . The onset of the heating mechanism is triggered in the weakly ionized region of the atmosphere when Γ makes the transition from being 1 to being 1. The quantity Γ increases with height from values 1 at the photosphere, reaches unity near the height of the local temperature minimum, increases with height to values 1 over a height range of a few hundred kilometers, and continues to increase with height as long as the atmosphere is weakly ionized. During this transition M i → Mp . Let E, B, and V be the electric and magnetic fields, and the CM velocity of the plasma. The CM electric field ECM = E + (V × B)/c provides the electromagnetic energy for resistive heating by doing work on the charges in the current density J. In MHD, the rate at which electromagnetic energy is transformed into thermal energy is J·ECM , which may be positive or negative. Let σP be the Pedersen conductivity for current flow along ECM,⊥ , defined as the component of ECM perpendicular to B. As indicated in Section 2, σP = σ (1 + Γ)/[(1 + Γ)2 + Me2 ]. Here σ is the conductivity for current flow parallel to B. The transition from a weakly magnetized to a strongly magnetized, but still weakly ionized atmosphere with increasing height from the photosphere causes σP /σ to decrease with height. It is ∼1 in weakly magnetized regions, and orders of magnitude 0. The initial profiles for p and ρ are the FAL profiles. The initial velocity is zero. It follows from the momentum equation that vx (z, t) = vy (z, t) = 0 for t > 0. The initial magnetic field is    z dη . (1) Bx (z, 0) = B0 exp −α 0 2L(η) Here L(z) = kB T (z)/ m(z) g is the ideal gas pressure scale height for a multi-specie plasma. The average particle mass m(z) ≡ ρ(z, 0)/n(z). Here kB , T (z), n(z), and g are Boltzmann’s constant, the FAL temperature and total number density, and the gravitational acceleration at the photosphere. The constant α = 0.7326 is chosen so that Bx (z, 0) decreases by a factor of 100 from z = 0 to z = 1.904 × 104 km. The choice of a non-equilibrium, non-force-free initial state resulting from a non-uniform, horizontal magnetic field Bx (z, 0) combined with an FAL hydrostatic equilibrium state is made for the following reasons. (1) The resistive heating rate is determined by ηP (z, t) ≡ (1 + Γ)η . Since Γ ∝ B 2 , and since Γ  1 throughout almost the entire chromosphere, ηP is essentially ∝ B 2 in the chromosphere. Then it is important to choose B(z, 0) in as realistic a way as possible within the limitations of the model to obtain a best estimate of ηP (z, 0) since doing so is important for obtaining a best estimate of ηP and the associated resistive heating rate for t > 0. The only height-dependent magnetic field allowed by the model is a horizontal field. It was decided to choose the simplest heightdependent magnetic field that has some basis in reality, and that decreases exponentially fast with increasing height with the local ideal gas pressure scale height to no more than a few tens of Gauss at the base of the TR. Orozco Su´arez et al. (2007) and Lites et al. (2008) present Hinode observations at a spatial resolution of 230 km that indicate that the magnetic field in the internetwork photosphere is mainly horizontal with magnitudes up to several hundred Gauss, a mean value of ∼100 G, and a filling factor of ∼20%, with the strongest fields concentrated in inter-granular lanes. In addition, the magnetic field in sunspot penumbrae is mainly horizontal. (2) Wave drivers in the real photosphere and chromosphere do not operate in an equilibrium or forcefree background atmosphere. They operate in, and couple to a non-equilibrium, non-force-free environment of convection and magnetic fields. Since the wave driver and initial state are nonlinearly coupled, the evolution of the atmosphere is not the linear superposition of the evolution of the initial state without the driver, and of a perturbation due to the driver. For example, the resistive heating rate cannot be decomposed into the sum of two resistive heating rates, one due to the un-driven evolution of the initial state, and one due to the driver. Although the dynamics cannot be separated into a sum of driven and un-driven components, an approximate way to check that the driver is the main cause of the resistive and compressive heating rates predicted by the simulation is as follows. The simulation with the same initial state is run without any wave driver, and the resistive and compressive heating rates per unit volume are computed at t = 200 s, which

Vol. 708

is the duration of the simulation. Integrating these heating rates over a given height range gives the corresponding heating fluxes. These fluxes are compared with those for the driven simulation. For the driven case, a train of six shock waves with magnetic field jumps ∼10–102 G is generated. This wave train extends from ∼375 km to the TR. The properties of this wave train are discussed in detail in Section 3.1. For the un-driven case, three wave-like structures appear in the magnetic field, two of which appear to be largely formed shock waves. The magnetic field changes by ∼1–2 G across these waves, so they are much weaker than the shocks generated in the driven case. This three wave train extends from ∼1100 km to 1600 km. First consider the resistive heating rates. For the driven case, the resistive heating fluxes generated by and below the shock wave train are 4.6 × 106 and 1.49 × 107 erg cm−2 s−1 . For the un-driven case, the resistive heating fluxes generated by and below the three wave train are 1.26 × 103 and 1.66 × 105 erg cm−2 s−1 . Then the resistive heating fluxes for the driven case are 102 –103 times greater than those for the un-driven case. Next consider the compressive heating rates. The resistive heating rate is always positive. The compressive heating rate at a point may be positive or negative depending on whether the plasma is undergoing compression or expansion. For an oscillatory driver, as is used here, the compressive heating rate is expected to oscillate between positive and negative values except in the shock waves, which compress the plasma. For the driven case, the compressive heating fluxes generated by and below the shock wave train are 1.24×109 and −5.57×109 erg cm−2 s−1 . For the un-driven case, the compressive heating fluxes generated by and below the three wave train are 2.16 × 106 and 4.39 × 108 erg cm−2 s−1 . Then the compressive heating flux generated by the waves for the driven case is ∼103 times greater than that for the un-driven case, and the compressive heating flux generated below the wave train for the driven case is ∼ − 10 times that for the un-driven case. The conclusion is that the resistive and compressive heating rates of the driven system are not significantly influenced by the heating rates associated with the un-driven initial state, and instead are determined by the evolution of the initial state due to its nonlinear coupling to the driver. 2.2.2. Basic Equations

The shock wave solutions are obtained by solving the following simplified form of the KG06 equations. The mass and momentum conservation equations are

and

∂ρ + (ρVz ) = 0 ∂t

(2)

  B2 ρGM  ∂(ρVz ) . + ρVz2 + p + x = − ∂t 8π (R + z)2

(3)

Here G, M , and R are the gravitational constant, solar mass, and solar radius. The prime denotes differentiation with respect to z. Since the model is applied at heights such that z  R , GM  /(R + z)2 ∼ g. Faraday’s law is c2 ∂Bx + (Vz Bx ) = ∂t 4π



(1 + Γ)  Bx σ

 .

(4)

The derivation of this particular form of Faraday’s law is outlined below in the discussion of Equation (15).

No. 1, 2010

MHD SHOCK WAVE HEATING OF THE CHROMOSPHERE

The energy equation is   ∂ ρVz2 3p Bx2 ρGM + + − (R + z) ∂t 2 2 8π    2 ρGM 1 Bx 5p 2 ρV + + − Vz = + (R + z) 2 z 4π 2   c 2  1 + Γ Bx Bx . 4π σ

(5)

The equation of state is p = ρkB T /mp . Here mp is the proton mass. The specification of the equations is completed by defining how Γ and σ are computed. This is done in the following section in the context of a general discussion of the conductivity tensor σ¯ used in Ohm’s law J = σ¯ · ECM . 2.2.3. Ohm’s Law

Ohm’s law may be written as ˆ J = σ E + σP ECM,⊥ + σH ECM × B.

(6)

Here Bˆ is a unit vector parallel to B. Mitchner & Kruger (1973) derive an electrical conductivity tensor for a three-component plasma of electrons (e), singly charged ions (i), and neutral particles (n). The degree of ionization of the plasma is arbitrary. The components of this tensor are given by σ = σP = σH =  Γ=

ρn ρ

ne e2 , me (νei + νen ) σ (1 + Γ) (1 + Γ)2 + Me2 −Me σ (1 + Γ)2 + Me2

2 Me Mi , Me =

, ,

ωe ωi , Mi = ∗ , (νei + νen ) νi

eB eB ∗ mn , ωi = ,ν = ωe = νin . me c mi c i mi + mn

(7)

profiles for the temperature and particle densities are used to determine the collision frequencies that enter σ¯ . The only part of σ¯ that is computed self-consistently as a function of height and time is the magnetic field strength. Explicit expressions for the collision frequencies as functions of temperature and particle densities are given in Goodman (2004a). Including a radiative cooling term in the energy equation is necessary for predicting observed emission, and determining radiative damping effects on MHD shock waves. A selfconsistent calculation of ionization fractions is necessary for a self-consistent determination of the Hall, Pedersen, and parallel conductivities since they are functions of the charged and neutral particle densities, in addition to depending on temperature and magnetic field strength. Conversely, the ionization fraction of a species is expected to be a sensitive function of the conductivities for the following reason. The ionization fraction is a sensitive function of temperature through the Planck function; the temperature depends on the resistive heating rate through the energy equation, and the resistive heating rate depends on the conductivities. If the resistive heating rate is small compared with the compressive heating rate, as might be the case in shock waves, it still has an effect on shock structure, which determines the compressive heating rate, and hence the temperature. Important effects of transport processes may be understood by using MHD models that do not include complications of self-consistent treatments of radiative processes. That is the approach taken here. Coupling this or other MHD models to the statistical equilibrium radiative transfer code developed by Uitenbroek (2001), which is based on the method of Rybicki & Hummer (1991, 1992), is a tractable next step for including quasi-realistic radiative transfer and ionization effects. Only a model that includes these effects can determine the fraction of the thermal energy generation rate due to compression and resistive and viscous dissipation that is converted into radiation that escapes to infinity. Ohm’s law may also be written in terms of the resistivity tensor η. ¯ It is the inverse σ¯ −1 of the conductivity tensor. The corresponding form of the Ohm’s law is ˆ ECM = η J + ηP J⊥ + ηH J × B.

(8)

(9)

Here e is the electron charge, ne is the electron number density, mn is the neutral particle mass, mi is the ion mass, me is the electron mass, νei is the electron–ion collision frequency, νen is the electron–neutral collision frequency, and νin is the ion– neutral collision frequency. Goodman (2004a) used an averaging procedure to represent the solar atmosphere as a three-species plasma of electrons, singly charged ions, and neutral particles for the purpose of computing the electrical conductivity tensor using Equations (7)–(9). The properties of the singly ionized species are derived from the masses and FAL number densities of protons, He ii, He iii, C, Si, Al, Mg, Fe, Na, and Ca. The properties of the single neutral species are derived from the masses and FAL number densities of H i and He i. The conductivity tensor also depends on B through Me and Mi . For the numerical results presented in Section 3, B = |Bx (z, t)|, and is computed self-consistently during the simulation. The model in KG06 does not include radiative cooling in the energy equation, or a calculation of ionization fractions. FAL

273

(10)

Here η = 1/σ , ηH = η Me = B/(ecne ), and ηP = (1 + Γ)η are the parallel, Hall, and Pedersen resistivities. In the chromosphere ηP  η . In the photosphere and corona ηP ∼ η . The Spitzer resistivity ηS ≡ η . The elements of the resistivity and conductivity tensors are connected by the relations σ =

1 −ηH ηP , σH = 2 , σP = 2 . 2 2 η ηP + ηH ηP + ηH

(11)

Since ηH does not explicitly involve collision frequencies it is sometimes called a collisionless resistivity, although it does not contribute directly to the resistive heating rate Q = Q + Q ⊥

(12)

2 = σ E 2 + σP ECM,⊥

(13)

= η J 2 + ηP J⊥2 .

(14)

274

GOODMAN & KAZEMINEZHAD

Although σH and ηH do not appear explicitly in Q, they determine Q indirectly. The ratio σH /σP = −ηH /ηP determines the angle that J⊥ makes with ECM,⊥ , and hence plays a role in determining Q⊥ = J⊥ · ECM,⊥ (Goodman 2004a). Either form of the Ohm’s law given by Equations (6) and (10) shows that electrical conduction, and the associated resistive dissipation, is highly anisotropic when the electrons or ions are magnetized. The conductivities and resistivities differ and vary by orders of magnitude throughout the solar atmosphere. It is only in weakly magnetized regions where σ¯ reduces to the Spitzer conductivity σ , since in such regions σH → 0 and σP → σ . Plots of typical variations of the conductivities, electron, proton, and heavy ion magnetizations, and relative collision frequencies from the photosphere to the lower corona are given in Goodman (2004a). Returning to the special case considered in this paper, for which the only spatial dependence is on z, and Bx (z, t) is the only nonzero component of B, the Ohm’s law Equation (10) reduces to

ECM = ηP yˆ − ηH zˆ Jy . (15) Computing ∇ × E = ∇ × (ECM − V × B/c), and substituting the result into Faraday’s law yields Equation (4). The Pedersen resistivity ηP = (1 + Γ)/σ may be written as σP /(σP2 + σH2 ). This shows that the Hall and Pedersen conductivities combine to yield the effective resistivity of the plasma that occurs in the Faraday law and energy Equations (4) and (5). In general, the Hall conductivity couples orthogonal components of B in Faraday’s law. This allows a component of B to generate components orthogonal to it. This effect does not occur for the special case considered here since only Bx is nonzero. It remains to specify the boundary conditions at the upper and lower boundaries of the computational domain to determine the solution to the model. 2.2.4. Boundary Conditions

The pressure and mass density at z = 0 are fixed at their FAL values. The magnetic field boundary condition at z = 0 defines the driver of the simulation. This boundary condition is Bx (0, t) = B0 + δBx sin ωt.

(16)

Here B0 is the value of the initial magnetic field strength at the photosphere given by Equation (1). The simulation results presented in Section 3 are generated using B0 = 500 G, δBx = 250 G, and ω = 2π/T , with T = 30 s. The upper boundary of the simulation domain is in the lower corona at z = 1.7 × 104 km. The pressure and density at this boundary are fixed at their FAL values. The values of Vz (z, t) at the upper and lower boundaries, and of Bx (z, t) at the upper boundary are not constrained. 3. NUMERICAL RESULTS The shock wave train generated by the simulation is analyzed in this section. Section 4 presents an analysis of electromagnetic, thermal, and CM kinetic energy flow and transformation based on substituting the numerical solution into the three corresponding energy equations. These equations are also used in Section 4 to obtain estimates of numerical error.

Vol. 708

The duration of the simulation is 200 s. The simulation becomes unstable for run-times much greater than 200 s.3 The simulation runs long enough to suggest that shock wave trains can extend over the height range of the chromosphere and be the sites of resistive and compressive heating rates sufficient to balance the chromospheric NRL. Since the model is one dimensional, the horizontal extent, or equivalently the filling factor of the shock-forming region cannot be modeled. Resistive and viscous dissipation, and mechanical compression are the local mechanisms that generate thermal energy by converting electromagnetic and CM kinetic energy into thermal energy. They appear in the general form of the thermal energy equation, given by ∂ + ∇ · ( V) = Q + Qc + Qvis − Qrad − ∇ · q. ∂t

(17)

Here = 3p/2 is the thermal energy per unit volume, and Qc ≡ −p∇ · V is the compressive heating rate per unit volume. This expression for Qc may be derived as follows. The rate per unit mass at which work is done compressing a fluid element is −p dρ −1 /dt = (p/ρ 2 )dρ/dt = −(p/ρ)∇ · V = Qc /ρ. The next to last equality follows from the mass conservation equation. When ∇ · V < 0 (>0) at a point, the fluid element flowing through that point undergoes heating (cooling) by compression (rarefaction). During compression, CM kinetic energy is converted into thermal energy, and conversely during rarefaction. The terms Qvis , Qrad , and q in Equation (17) are the viscous heating and radiative cooling rates per unit volume, and the total diffusive thermal energy flux. These three terms are omitted in the model. It is well known that radiative cooling is important in determining shock structure and damping in the photosphere and lower chromosphere (Withbroe & Noyes 1977, and references therein). It is expected that in a more accurate model of shock structure that includes these omitted terms, gradients of temperature and velocity are so large in the shock layer that viscous heating and thermal energy diffusion become important. Khodachenko et al. (2004, 2006) provide estimates showing that the effects of viscosity and thermal conduction on linear MHD waves in the chromosphere are minor compared with the effect of anisotropic electrical conduction in determining the wave damping rate, and state that preliminary rough estimates suggest that thermoelectric current drive might be as important as electrical conductivity in determining the damping rate. Then heating by linear wave dissipation can probably be accurately described by models that include only the transport process of anisotropic electrical conduction, and possibly thermoelectric current drive. By contrast, an accurate description of thermal energy generation in shock layers probably requires including the full tensor transport effects of electrical and thermal conduction, viscosity, and thermoelectric effects in the magnetized chromosphere. Including all these tensor effects in an MHD 3

The instability is due to the steepening of one or more shock fronts beyond the point at which gradients at the fronts can be resolved with the spatial resolution of 1 km. The gradients are largest at a shock front, which is trailed by a shock layer in which the gradients rapidly decrease in magnitude with increasing distance from the shock front. A Godunov viscosity is included in the code to damp out numerical oscillations, increasing code stability. Gradients eventually become too large for this mechanism to stabilize the code. Numerical oscillations grow, the time step → 0 as determined by the Courant condition, and the simulation breaks down. Increasing the Godunov viscosity, which is a form of artificial diffusion, increases code stability and run time at the expense of increasing the amount of unphysical dissipation. Simulations that use artificial diffusivities orders of magnitude greater than what corresponds to solar conditions incur large errors in predicted heating rates.

No. 1, 2010

MHD SHOCK WAVE HEATING OF THE CHROMOSPHERE

275

1

3

10

10

0

10

t=200 s t=0.38 s

−1

10

2

10

−2

x

B (G)

|Jy|(A−m−2)

10

−3

10

−4

10 1

10

−5

10

t=200 s t=0.38 s

−6

10

−7

10

0

10

0

250

500

750

1000

1250 1500 z(km)

1750

2000

2250

0

250

500

750

1000

2500

Figure 1. Magnetic field vs. height at the beginning and end of the simulation. (A color version of this figure is available in the online journal.)

1250 1500 z(km)

1750

2000

2250

2500

Figure 3. Magnitude of the current density vs. height at the beginning and end of the simulation. (A color version of this figure is available in the online journal.)

−9

10

2

10 −10

10

0

10 −11

−2

10 −1

−12

10

−3

ergs−cm −sec

Resistivity (sec)

10

−13

10

ηP at t=200 s

−14

10

ηP at t=0.38 s

−6

10

−8

10

−10

η||=ηP/(1+Γ)

−15

10

−4

10

10

Q at t=200 s Q /Q = Γ/(1+Γ) at t=200 s

−12

Γ

10

Q at t=0.38 s

−16

10

0

250

500

750

1000

1250 1500 z(km)

1750

2000

2250

2500

Figure 2. Total resistivity ηP vs. height at the beginning and end of the simulation, and parallel resistivity η vs. height. η is independent of time since FAL data are used to compute it. (A color version of this figure is available in the online journal.)

shock wave simulation has never been done, and would have to be done at a resolution of 1 km to accurately resolve the transport physics. This is a difficult task for the future. This paper focuses on compressive and anisotropic resistive heating in shock waves at a resolution of 1 km in the simple magnetic field geometry of a height-dependent horizontal magnetic field as a step toward a more comprehensive description of transport processes, especially magnetization-induced transport processes in shock waves in the photosphere and chromosphere. For the case considered here, Q = Q⊥ . Then Q = ηP J⊥2 = (1 + Γ)η J⊥2 ≡ QS + QΓ . Here QS = η J⊥2 is the Spitzer heating rate, and QΓ = ΓQS is the heating rate due to the magnetizationinduced resistivity Γη . 3.1. Shock Wave Train Figures 1–14 show properties of the shock wave train. The TR in these figures lies between 2150 and 2200 km. Figure 1 shows the magnetic field near t = 0 and at t = 200 s. The initial field gives rise to a Lorentz force −Bx Bx /4π > 0 that tends to accelerate plasma upward. There are six shock waves.

−14

10

0

250

500

750

1000

1250 1500 z(km)

1750

2000

2250

2500

Figure 4. Total resistive heating rate per unit volume (Q) vs. height at the beginning and end of the simulation, and QΓ /Q vs. height at the end of the simulation. QΓ ∝ Me Mi is the contribution to Q that explicitly depends on particle magnetizations. (A color version of this figure is available in the online journal.)

The leading shock wave reaches z ∼ 2000 km. Magnification of individual shocks shows they have thicknesses that increase from ∼10 km to ∼20 km with increasing height. The shock strength, defined as the jump in Bx across the shock front, divided by the field strength immediately ahead of the shock front, is ∼0.51, 0.84, 0.74, 0.66, 0.66, and 1.14 for the six shocks in order of increasing height. The shock wave near z = 103 km is distorted by an instability that is probably part of the cause of the eventual breakdown of the simulation. Magnification shows there is a clearly defined shock front near z = 103 km. The shock waves form near z = 375 km. This is the height at which ηP ∼ 1.1η . The parallel and Pedersen resistivities are shown in Figure 2. The increase of Me Mp with height causes ηP to exceed η by a factor of ∼2 × 103 in the upper chromosphere. Outside of the shocks, ηP is mostly smaller than the initial resistivity. At t = 200 s, Γ first reaches unity for z ∼ 500 km. In the underlying region where Γ  1, the conductivity tensor is isotropic and reduces to σ , and J is mainly an electron current density. In the overlying region where Γ  1, the conductivity

276

GOODMAN & KAZEMINEZHAD 5

Vol. 708

18

10

10

4

− p∇⋅ V < 0 (Rarefaction) − p∇⋅ V > 0 (Compression)

10

3

t=200 s t=0.38 s

16

10

10

2

1

14

10

10 n(cm−3)

−3

ergs−cm −sec

−1

10

0

10

−1

12

10

10

−2

10

−3

10

10

10

−4

10

−5

10

8

0

250

500

750

1000

1250 1500 z(km)

1750

2000

2250

10

2500

Figure 5. Compressive heating rate per unit volume (Qc ) vs. height at the end of the simulation. (A color version of this figure is available in the online journal.) 12

2.5

10

2

10

8

1.5

6

1

10

−1

10−2 g−cm−2−sec−1

−1

10

4

10

2

10

0

10

500

750

1000

1250 1500 z(km)

1750

2000

2250

2500

F (0,t) m

F (z m

,t) × 106

max

0.5 0 −0.5 −1

−2

10

Qm at t=200 s

−4

10

−1.5

Qm at t=0.38 s

−6

−2

10

−2.5

−8

10

250

Figure 8. Number density vs. height at the beginning and end of the simulation. (A color version of this figure is available in the online journal.)

10

ergs−g −sec

0

0

250

500

750

1000

1250 1500 z(km)

1750

2000

2250

2500

Figure 6. Total resistive heating rate per unit mass (Qm ) vs. height at the beginning and end of the simulation. (A color version of this figure is available in the online journal.)

0

20

40

60

80

100 120 Time(sec)

140

160

180

200

Figure 9. Mass flux through the lower (photosphere) and upper boundaries vs. time. (A color version of this figure is available in the online journal.)

13

10

4 2

12

10

0 11

−2 Vz (km−sec−1)

ergs−g−1−sec−1

10

10

10

9

10

8

−4 −6 −8 −10

10

−12 7

− p∇⋅ V/ρ < 0 (Rarefaction) − p∇⋅ V/ρ > 0 (Compression)

10

6

10

0

250

500

750

1000

1250 1500 z(km)

1750

t=200 s t=0.38 s

−14 2000

2250

2500

Figure 7. Compressive heating rate per unit mass (Qcm ) vs. height at the end of the simulation. (A color version of this figure is available in the online journal.)

−16

0

250

500

750

1000

1250 1500 z(km)

1750

2000

2250

2500

Figure 10. Velocity vs. height at the beginning and end of the simulation. (A color version of this figure is available in the online journal.)

No. 1, 2010

MHD SHOCK WAVE HEATING OF THE CHROMOSPHERE

6

277

2

10

10

1

10

5

10

0

t=200 s t=0.38 s

4

10

t=200 seconds

10

−2

p(dynes−cm )

−1

10 3

10

−2

10 2

10

−3

10

−4

1

10

10

−5

10

0

10

−1

|Ey+VzBx/c|(V−m )

−6

10

−1

10

0

250

500

750

1000

1250 1500 z(km)

1750

2000

2250

2500

Figure 11. Pressure vs. height at the beginning and end of the simulation. (A color version of this figure is available in the online journal.) 6

−1

|E |(V−m ) y

−7

10

0

250

500

750

1000

1250 1500 z(km)

1750

2000

2250

2500

Figure 14. Magnitudes of Ey and ECM,y vs. height at t = 200 s. (A color version of this figure is available in the online journal.)

10

5

10

t=200 s t=0.38 s

4

T(K)

10

3

10

2

10

1

10

0

250

500

750

1000

1250 1500 z(km)

1750

2000

2250

2500

Figure 12. Temperature vs. height at the beginning and end of the simulation. (A color version of this figure is available in the online journal.) 45 35 25 15 5 −5 −15 −25 −35 −45 −55 −65

t=200 s

−75

(E +V B /c)(10−2V−m−1)

−85

y

−95 Minimum value ~ −227

z x −1

E (V−m ) y

−105 −115

0

250

500

750

1000

1250 1500 z(km)

1750

2000

2250

2500

Figure 13. Electric field Ey , and CM electric field ECM,y vs. height at t = 200 s. (A color version of this figure is available in the online journal.)

tensor is anisotropic, J is an electron current density, and J⊥ is a proton Pedersen current density that flows parallel to ECM,⊥ .

The smooth photospheric waves that drive the simulation and the resulting shock waves are perpendicular propagating waves in that B is perpendicular to the direction of propagation. The smooth waves can steepen into shocks only if their characteristic speed, determined by their interaction with the medium through which they propagate, decreases with height. In general, a continuous wave steepens into a stable shock if and only if the propagation velocity of the wave decreases in the direction of wave propagation (Whitham 1974, Section 2.6). Although the simulation is nonlinear and includes an electrical conductivity tensor and the effect of gravity, the driver waves are analogous to fast, linear magnetoacoustic waves in ideal MHD that steepen into shocks. A detailed analysis of the linear dispersion relations for parallel and perpendicular propagating waves that include the effects of gravity and the electrical conductivity tensor in an FAL background atmosphere is presented in Sections 3.2 and 3.3 of KG06. The characteristic wave speed of the smooth waves is close to the fast magnetoacoustic wave speed Vm ≡ (VA2 + 2 )1/2 . Here VA = B/(4πρ)1/2 and Vsound = (5p/3ρ)1/2 Vsound are the Alfv´en and sound speeds (Freidberg 1982). Then Vm must decrease with height in order that the smooth-driven waves steepen into shock waves with increasing height. At the start of the simulation, dVm /dz < 0 for z  250 km, and dVm /dz > 0 for z > 250 km. This allows the smooth-driven waves to begin +/− to steepen into shocks in the region z  250 km. Let Vm be the values of Vm immediately in front of and behind a shock wave, respectively. As time increases, the forming shock waves modify the surrounding plasma such that Vm− − Vm+  1 km s−1 for each shock. The shock waves become fully formed and continue to propagate upward, modifying their own and each others environment so that from the photosphere to the upper chromosphere the shocks always propagate into a medium in which Vm− > Vm+ . Figure 3 shows the magnitude of the current density. It increases by 2–3 orders of magnitude in the shock waves relative to its external value. The shocks are propagating layers of relatively intense current. Outside of the shocks the current density is mostly smaller than the initial current density. Figure 4 shows the resistive heating rate per unit volume Q, and the fraction QΓ /Q ∼ Me Mi /(1 + Me Mi ) of the heating rate due to magnetization effects. Q increases by 4–6 orders of magnitude in the shock waves relative to its external value.

278

GOODMAN & KAZEMINEZHAD

QΓ /Q is essentially unity for z  750 km. For z  500 km, Q is mainly due to dissipation of electron currents. For 500  z  600 km, Q is mainly due to dissipation of proton Pedersen currents. For z  600 km, Q is entirely due to dissipation of proton Pedersen currents. ECM,⊥ provides the electromagnetic energy for resistive heating. It is generated by a combination of induction effects according to Faraday’s law when ∂B/∂t = 0, and by convection effects due to CM flow orthogonal to B, which generates a convection electric field. The chromospheric heating flux F due to resistive heating is obtained by integrating Q at t = 200 s over the height range of the shock wave train, which is 375 km  z  2100 km. The result is F = 4.6 × 106 erg cm−2 s−1 . Integrating Q over the height range 0 km  z  2100 km gives the total resistive heating flux Ftotal = 1.95 × 107 erg cm−2 s−1 . Then ∼76.4% of the resistive heating occurs below the shock wave train, in the region where the wave generated by the driver is still a smooth wave. Figure 4 shows that almost all of this heating occurs close to the driver at z = 0. At z = 375 km, Me Mi ∼ 0.1, so the plasma is not yet magnetized, and resistive heating is mainly due to dissipation of electron currents. At z = 500 km, Me Mi ∼ 1, and heating by proton Pedersen current dissipation begins to dominate the resistive heating rate. There is only one shock wave below z = 500 km. It has the smallest resistive heating rate of the six shock waves. About 90% of F is due to proton Pedersen current dissipation in the shocks above z = 500 km. Figure 5 shows the compressive heating rate per unit volume Qc . It is ∼10–103 Q within shock layers at points where Qc > 0. Between shock waves, and at some points within shock layers, Qc < 0, indicating cooling by rarefaction. The chromospheric heating flux Fc due to compressive heating is obtained by integrating Qc over the height range of the shock wave train. The integral begins at z = 370 km, which is the first point at which Qc > 0 at the base of the lowest altitude shock wave. Integrating over points at which Qc > 0 ( 100 s. The reason for this is not clear, although only three of the six shock waves are formed by t = 100 s. This change in sign of the time average of Fm (0, t) might indicate the emergence of an oscillation period longer than the driver period due to nonlinear dispersive effects. The simulation must be run for much longer than 200 s to determine if this is the case. In Figure 9, the amplitude of Fm (0, t) ∼ 10−2 g cm−2 s−1 . This corresponds to a velocity amplitude of ∼0.1–1 km s−1 , which is reasonable. Then the predicted values of Fm (0, t) are reasonable if the FAL density is a reasonable value for the actual photospheric density. However, the amplitude of Fm (0, t) is ∼109 times larger than the observed solar wind mass flux into interplanetary space. The related reasons for this difference are as follows. (1) The mass flux into the upper atmosphere is Fm (zmax , t). Its value determined by the shock wave train is not known until the shock wave train reaches zmax , and until the average of Fm (zmax , t) over some appropriately chosen time interval T, discussed below, becomes constant. This state is defined as the dynamic steady state of the simulation. The simulation is not in this state by t = 200 s since the shock wave train does not extend to zmax by this time. Even if the simulation can be run long enough to reach a dynamic steady state, the requirement that some time average of Fm (zmax , t) approach the solar wind mass flux imposes a strong restriction on Fm (0, t). This is shown as follows. Choose times t2 > t1 where t1 is after

279

the shock wave train reaches zmax . Let T = t2 − t1 . The spatial integral of the time average of the mass conservation equation gives (M(t1 + T ) − M(t1 ))/T = Fm (0, t) T − Fm (zmax , t) T . Here T denotes the time average over t1  t  t1 + T , and M(t) is the total mass in the computational domain. Assume M is bounded. Then the left-hand side of this equation is zero if M(t1 + T ) = M(t1 ), as is the case if M(t) becomes periodic with period T. Even if M does not become periodic, the left-hand side must approach zero for fixed t1 , and T → ∞ since M is bounded. The point is that after a sufficiently long time the condition Fm (0, t) T = Fm (zmax , t) T must hold. Assuming Fm (0, t) ∼ 10−2 Vz (km s−1 ) g cm−2 s−1 , and requiring Fm (zmax , t) T ∼ 10−11 g cm−2 s−1 , which is a typical solar wind flux from regions with a closed magnetic field, it follows that Vz (0, t) T ∼ 10−9 km s−1 , which is essentially zero. Then in a dynamic steady state Vz (0, t) must be almost exactly periodic with period T. This condition is satisfied if Vz (0, t) has any combination of periods Ti such that T /Ti is any positive integer. For a monochromatic driver, it is expected that the appropriate choice for T is the driver period. Even if ρ(z = 0) is fixed, and the system is driven by a photospheric magnetic field driver with period T, Vz (0, t) might not become periodic with period T due to nonlinear dispersive effects. One way to enforce the generation of a realistic solar wind flux is to fix ρ(z = 0), and drive the system with a periodic, multiply periodic, or random vertical velocity such that the statistics of the mass flux through the photosphere, such as its mean, standard deviation, and fluctuation timescales are comparable to those of the solar wind. This is an unrealistic way to enforce a realistic mass flux in a one-dimensional model since the mass flux through a given point in the real photosphere is not expected to have a time average close to the solar wind mass flux given the relatively enormous amplitude of photospheric mass fluxes, as estimated above. This leads to the next point. (2) The restriction of the model to one dimension prevents mass from returning to the photosphere by a combination of horizontal and downward flow. The time average of the mass flux through a point in the photosphere can, and is expected to be, much larger than the solar wind mass flux provided most of the mass driven upward returns to the photosphere. Since the predicted values of Fm (0, t) are reasonable, mass fluxes through the photosphere can be many orders of magnitude larger than the solar wind flux, but almost all of this mass returns to the photosphere. (3) The radial expansion of the magnetic field with increasing height is a multi-dimensional effect that is important in determining the solar wind mass flux. These reasons imply that a multidimensional model that uses spherical coordinates, that has an upper boundary at least one solar radius above the photosphere, and that is run until a dynamic steady state is reached must be used to obtain a relevant prediction of the rate of mass loss into interplanetary space. Similar remarks, based on the thermal energy Equation (17), apply to the rate of energy flow into interplanetary space. The shocks have infinitesimal thickness on the scale of the chromosphere and the local pressure scale height. The MHD jump condition corresponding to the Faraday law Equation (4) may then be applied on either side of the shock fronts where conditions are relatively uniform to obtain an estimate of the shock speeds. This jump condition is  c2 (1 + Γ)Bx Vshock [Bx ] = Vz Bx − . 4π σ 

(18)

280

GOODMAN & KAZEMINEZHAD

Vol. 708

14

4 Q at t=200 s QΓ/Q = Γ/(1+Γ) at t=200 s

3.5

Qm at t=200 s 12

Q x 103 at t=0.2 s

3

m

Q x 10 at t=0.2 s 3

109 ergs−g−1−sec−1

−3

ergs−cm −sec

−1

10 2.5 2 1.5

6

4

1 0.5

2 900

910

920

Figure 15. Resistive heating rate per unit volume in the shock wave near z = 885 km at the beginning and end of the simulation. (A color version of this figure is available in the online journal.)

Here Vshock is the shock velocity, [f ] ≡ f (z+ ) − f (z− ) is the jump in any quantity f across the shock, and z± denotes points immediately ahead of and behind the shock wave front. Since the shocks have thicknesses ∼10–20 km, z± are visually approximated as positions immediately outside the region where f is changing rapidly. The term in Equation (18) that involves σ is insignificant outside the shocks compared with Vz Bx . This is because Bx is orders of magnitude larger inside the shock than outside the shock. Figure 10 shows [Vz ] across the shocks. Equation (18) predicts that for the shocks from the lowest to the greatest height, Vshock (km s−1 ) ∼5.76, 7.25, 6.32, 4.10, 4.67, and 3.62 at t = 200 s. These estimates of Vshock are consistent with estimates obtained directly from the simulation. Figures 11 and 12 show the pressure and temperature. They are increased by the passage of a shock wave, along with the density, flow velocity, and magnetic field. Above z ∼ 1000 km, the temperature is far below the FAL temperature. This is mainly due to the equation of state T = mp p/kB ρ, and the relatively large values of ρ generated by the simulation. Figure 13 shows the electric field Ey , and the CM electric field ECM,y = Ey + Vz Bx /c. The spatial variation of Ey is determined by the Faraday law c ∂Ey /∂z = ∂Bx /∂t. Across a shock wave, this equation yields the jump condition [Ey ] + Vshock [Bx ]/c = 0.

(19)

This predicts the jumps in Ey shown in Figure 13. This jump condition may be written as [Ey + Vshock Bx /c] = 0. This states that in the frame of reference moving with a given shock wave, the jump in ECM,y across the shock wave is zero. The field ECM,y 2 determines the resistive heating rate Q = ECM,y /ηP . Figure 13 shows that |ECM,y | has sharply peaked maxima in the shock waves, where it reaches values ∼0.1–1 V m−1 . This is seen more clearly in Figure 14. Figure 14 shows the magnitudes of Ey and ECM,y . The |Ey | decreases by factors of ∼5–100 in the shock waves relative to its value outside the shock waves. The |Ey | is a local minimum in the shock waves. Conversely, |ECM,y | is a local maximum in the shock layers. It increases by orders of magnitude in the layers relative to its value outside the layers. It is only in the shock layers where |ECM,y | is large enough to drive a significant resistive heating rate.

0 860

870

880

890 z(km)

900

910

920

Figure 16. Resistive heating rate per unit mass in the shock wave near z = 885 km at the beginning and end of the simulation. (A color version of this figure is available in the online journal.) 6

−1

890 z(km)

5

Q at t=200 s

4

Q × 103 at t=0.38 s

c c

3 2

−3

880

1 0

2

870

10 ergs−cm −sec

0 860

8

−1 −2 −3 −4 860

870

880

890 z(km)

900

910

920

Figure 17. Compressive heating rate per unit volume in the shock wave near z = 885 km at the beginning and end of the simulation. (A color version of this figure is available in the online journal.)

Figures 4 and 6 show that Q and Qm decrease to insignificant values through the TR and into the lower corona. The reason is that Γ → 0 due to (ρn /ρ)2 → 0, although Me Mp continues to increase through the TR and into the corona. Then Q → η J⊥2 , and J becomes an electron current density. Electron Pedersen current dissipation may play a significant role in coronal heating. This is seen by noting that when Γ = 0, Q may be written as 2 σ E 2 + σP ECM,⊥ . This explicitly shows the dependence of Q on the Pedersen conductivity, which is determined by electron– proton collisions in the TR and corona. 3.2. Internal Structure of a Shock Wave Figures 15–25 show the internal structure of the shock wave near z = 885 km at t = 200 s. Figures 15 and 16 show that the resistive heating rates per unit volume and mass are confined to a layer with a thickness of ∼10 km. QΓ /Q has a minimum value outside the shock layer of 0.9374, and a maximum value in the shock layer of 0.981, corresponding to values of Γ ∼ 15 and 51.6, and resistivities

No. 1, 2010

MHD SHOCK WAVE HEATING OF THE CHROMOSPHERE 65

2 Qcm at t=200 s 1.5

Q

cm

t=200 s t=0.2 s

60

× 103 at t=0.38 s

55

1

1012 ergs−g−1−sec−1

281

50 Bx(G)

0.5

0

45 40

−0.5

35 −1

−1.5 860

30

870

880

890 z(km)

900

910

25 860

920

Figure 18. Compressive heating rate per unit mass in the shock wave near z = 885 km at the beginning and end of the simulation. (A color version of this figure is available in the online journal.)

870

880

890 z(km)

900

910

920

Figure 20. Magnetic field in the shock wave near z = 885 km at the beginning and end of the simulation. (A color version of this figure is available in the online journal.)

1

25 t=200 s t=0.2 s

0.5

20 15 10

−2

Jy(A−m )

5 0

0

−5 −10 −15

−0.5

−20 −25 −30

−1

−35

Ey(V−m−1)

−40

E

−1.5 860

870

880

890 z(km)

900

910

920

Figure 19. Current density in the shock wave near z = 885 km at the beginning and end of the simulation. (A color version of this figure is available in the online journal.)

ηP = (1 + Γ)η ∼ 16 η and 53 η . Integrating Q(z) over this shock wave gives a heating flux of 7 × 105 erg cm−2 s−1 . Figures 17 and 18 show the compressive heating rates per unit volume and mass. The regions of compressive heating and cooling are confined to a layer with a thickness of ∼20 km. Integrating Qc over this shock wave gives a total compressive heating flux of 4.2 × 107 erg cm−2 s−1 . The regions of heating and cooling contribute fluxes of 1.4 × 108 and −9.8 × 107 erg cm−2 s−1 to the total flux. Figures 19 shows that the current density is confined to a layer with a thickness of ∼15 km, consistent with the variation of Bx shown in Figure 20. Figure 21 shows the variation of Ey , ECM,y , ηP , and the convection electric field Econvection = Vz Bx /c across the shock layer. The maximum values of |ECM,y |, Q, and ηP occur at z = 887 km. There ECM,y = −0.5666 V m−1 , Ey = −0.8370 V m−1 , Econvection = 0.2704 V m−1 , ηP = 9.6455 × 10−11 s, and Q = 3.6982 erg cm−3 s−1 , as shown in Figure 15. A way to simultaneously estimate the shock velocity Vshock and test the accuracy of the simulation is as follows. Faraday’s law implies the jump condition [Ey ] + Vshock [Bx ]/c = 0,

=(V B /c)(V−m−1)

convection

−45

z x −2

−1

−50

(Ey+VzBx/c)(10 V−m )

−55

η (10

−60 860

−11

P

865

870

875

880

885

890 895 z(km)

900

sec) 905

910

915

920

Figure 21. Electric fields and resistivity in the shock wave near z = 885 km at t = 200 s. (A color version of this figure is available in the online journal.)

discussed in the context of Figure 13. The mass conservation Equation (2) implies the jump condition [nVz ] = Vshock [n].

(20)

The momentum conservation Equation (3) implies the jump condition  Vshock [ρVz ] = ρVz2 + p + Bx2 /8π . (21) The energy conservation Equation (5) implies the jump condition  Vshock ρVz2 /2 + 3p/2 + Bx2 /8π − ρgR 

= Vz ρVz2 /2 + 5p/2 + Bx2 /4π − ρgR − cηP Bx Jy /4π . (22) Figures 20–24 allow the jumps in Bx , Ey , ρ, Vz , and p to be computed once the corresponding points z± immediately ahead

282

GOODMAN & KAZEMINEZHAD 220

2.6

200

2.4 t=200 s t=0.2 s

2.2

160 p(dynes−cm−2)

1.8

14

−3

n(10 cm )

t=200 s t=0.2 s

180

2

1.6

140 120

1.4

100

1.2

80

1

60

0.8 860

Vol. 708

870

880

890 z(km)

900

910

40 860

920

Figure 22. Number density in the shock wave near z = 885 km at the beginning and end of the simulation. (A color version of this figure is available in the online journal.)

870

880

890 z(km)

900

910

920

Figure 24. Pressure in the shock wave near z = 885 km at the beginning and end of the simulation. (A color version of this figure is available in the online journal.) 6.5

1

6

0 −1

5.5 t=200 s t=0.2 s

−3

5

3

−1

T(10 K)

−2 Vz (km−sec )

t=200 s t=0.2 s

−4

4.5

−5

4 −6

3.5

−7 −8 −9 860

3 860 870

880

890 z(km)

900

910

920

Figure 23. Velocity in the shock wave near z = 885 km at the beginning and end of the simulation. (A color version of this figure is available in the online journal.)

of and behind the shock wave are defined. It follows from Figures 20–24 that an appropriate choice for z+ is 890 km. Since the shock has finite thickness, and merges smoothly into the downstream plasma, the appropriate choice for z− is not as clear. If the simulation is accurate there is a single value of z− at which the jump conditions in Equations (19)–(22) predict essentially the same shock velocity. Such a value of z− exists, and is determined as follows. Define Vs1 , Vs2 , Vs3 , and Vs4 to be the values of Vshock computed from the jump conditions in Equations (19)–(22). Define d = (Σ4i,j =1 (Vsi − Vsj )2 )1/2 . Evaluate d for 860 km  z−  889 km, assuming z+ = 890 km. It is found that d decreases monotonically with increasing z− from 17.45 km s−1 at z− = 860 km to a minimum value of 0.963 km s−1 at z− = 887 km. For z−  888 km, some of the jump conditions predict negative shock velocities, which is not physical. Then z− = 887 km, and is identified as being immediately behind the leading edge of the shock wave. The values of Vshock at z− predicted by the jump conditions are Vs1 = 4.3922, Vs2 = 4.4034, Vs3 = 4.9557, and Vs4 =

870

880

890 z(km)

900

910

920

Figure 25. Temperature in the shock wave near z = 885 km at the beginning and end of the simulation. (A color version of this figure is available in the online journal.)

4.4035 km s−1 . Their average value is 4.5387 km s−1 . Q is a maximum at z = z− . The sound, Alfv´en, and fast magnetoacoustic wave speeds at z = z+ are 6.8215, 5.4384, and 8.7241 km s−1 . The flow velocity Vz at this height is −8.7183 km s−1 . Then the shock wave velocity relative to the medium immediately ahead of the shock wave is Vshock − Vz ∼ 4.5387 + 8.7183 = 13.257 km s−1 . This is ∼1.52 times the fast magnetoacoustic wave speed. The sound, Alfv´en, and fast magnetoacoustic wave speeds at z = z− are 8.4782, 7.0607, and 11.033 km s−1 . At this height, the flow velocity is −3.5449 km s−1 , and Vshock − Vz ∼ 8.0836 km s−1 . Then the shock velocity relative to the flow velocity is supermagnetoacoustic at z+ , and sub-magnetoacoustic but superAlfv´enic at z− . This is consistent with the interpretation of the shock waves as fast magnetoacoustic shock waves. At z = 876 km, the four shock speeds computed from the jump conditions are 5.2276, 3.6214, 6.8532, and 3.6214 km s−1 . They are within a factor of 2 of one another, and close to their mean of 4.8309 km s−1 , which is close to the mean at z = 887 km, computed above. For decreasing values of z, the

No. 1, 2010

MHD SHOCK WAVE HEATING OF THE CHROMOSPHERE

Energy flow and transformation in the solar atmosphere is a complex process. Fluxes of thermal, CM kinetic, and electromagnetic energy flow in both directions through the photosphere, and are transformed into one another in the overlying atmosphere. Transport processes play a major role in determining the fluxes and transformation rates of these forms of energy. The thermal energy conservation equation is Equation (17), omitting Qvis , Qrad , and ∇ · q. The CM kinetic and electromagnetic energy conservation equations are ∂ ∂t



ρV 2 2



   ρV 2 = −∇ · p+ V − Qc + V · Fm + ρV · g 2 (23)

and ∂ ∂t



B2 8π

−3

4. ENERGY FLOW AND TRANSFORMATION, AND ESTIMATES OF NUMERICAL ERROR

ergs−cm −sec

−1

four predicted values of the shock speed steadily diverge. This suggests that the shock extends from z ∼ 876 km to z ∼ 887 km, consistent with visual inspection of the figures. The temperature is shown in Figure 25. It is computed from the ideal gas equation of state.

11 10 9 8 7 6 5 4 3 2 1 0 −1 −2 −3 −4 −5 −6 −7 −8

283 〈∂ε/∂ t 〉 z 〈 −∇⋅(ε V)−p∇⋅ V+Q 〉 z

0

20

40

60

80

100 t (sec)

120

140

160

180

200

Figure 26. Numerical error as measured by the degree to which the numerical solution satisfies the spatial average of the thermal energy equation vs. time. The difference between the two curves is due to numerical error. (A color version of this figure is available in the online journal.) 0.8



2

0.7

〈 ∂(ρ V /2)/∂ t 〉 z

0.6

〈 − ∇⋅(p V+ VρV /2) + p∇⋅ V+ V⋅ Fm+ρ V⋅ g 〉 z

2

= −∇ · S − V · Fm − Q.

(24)

0.5 ergs−cm−3−sec−1

0.4

Here Fm = (J × B)/c is the magnetic Lorentz force per unit volume, and S is the Poynting flux. The thermal, CM kinetic, and electromagnetic energy conservation equations are each coupled to the other two by the pairs of terms (Qc , Q), (−Qc , V · Fm ), and (−V · Fm , −Q), respectively. These terms sum to zero, and describe the flow of energy between the three energy reservoirs. The energy conservation Equations (17), (23), and (24) are not part of the model in that they are not among the equations that are solved to determine the state of the plasma. However, these energy equations may be mathematically derived from the equations of the model. Then, if there is no numerical error, the numerical solution of the model satisfies Equations (17), (23), and (24) exactly. This is never the case since there is numerical error in every simulation. If this error is sufficiently small, then using the numerical solution to compute the terms in Equations (17), (23), and (24) provides insight into how energy flows through the plasma, and is transformed from one form to another. The predictions of any simulation must be understood in the context of numerical error, and of the degree to which the model on which the simulation is based contains an accurate description of the relevant transport processes. It is emphasized that the degree to which the numerical solution of the model equations satisfies the energy equations does not say anything about the physical content of the model equations. It is only a measure of how accurately the model equations are solved numerically. The model equations used here omit important effects, but as shown in Section 4.1 their numerical solution satisfies the energy equations derived from the model equations reasonably well. It is necessary to make the distinction between solving a set of model equations accurately, and the physical content of the model equations. In addition to the energy conservation test just discussed, the fact that the MHD Rankine–Hugoniot relations predict shock speeds in Sections 3.1 and 3.2 that are close to one another is another indication that the numerical solution is reasonably accurate.

0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 −0.5

0

20

40

60

80

100 t (sec)

120

140

160

180

200

Figure 27. Numerical error as measured by the degree to which the numerical solution satisfies the spatial average of the CM kinetic energy equation vs. time. The difference between the two curves is due to numerical error. (A color version of this figure is available in the online journal.)

4.1. Estimates of Numerical Error Figures 26–28 show the spatial averages, denoted by z , of Equations (17), (23), and (24) over the computational domain as functions of time. The time derivatives of the thermal, CM kinetic, and electromagnetic energy densities are plotted separately from the sum of the other terms in the equations. The difference between the two curves in each plot is a measure of numerical error. The overall error is small to moderate. Most of the error occurs for t > 100 s. At t = 100 s two shock waves are fully formed, a third is forming, and the leading shock wave is at 1040 km. The error is believed to arise mainly from errors in the calculation of spatial derivatives in the shock layers, which are not fully resolved by the spatial resolution of 1 km. Numerical error is expected to increase with the order of the spatial derivatives that are computed. The model involves second-order spatial derivatives through terms involving Bx in Equations (4) and (5).

284

GOODMAN & KAZEMINEZHAD 1

Vol. 708

10

0.75

〈∂ε /∂ t 〉

〈 ∂(B2/8π)/∂ t 〉 z

9

〈 −(∇⋅ S+ Fm⋅ V+Q)〉 z

8

〈 −∇⋅(ε V)〉

7

〈 −p∇⋅ V 〉 z

6

10 〈 Q 〉

0.5

z z

2

z

0.25 −1

4 ergs−cm −sec

0

−3

ergs−cm−3−sec−1

5

−0.25 −0.5

3 2 1 0 −1 −2 −3

−0.75

−4 −5

−1

0

20

40

60

80

100 t (sec)

120

140

160

180

200

Figure 28. Numerical error as measured by the degree to which the numerical solution satisfies the spatial average of Poynting’s equation vs. time. The difference between the two curves is due to numerical error. (A color version of this figure is available in the online journal.)

Among Equations (17), (23), and (24), the only second-order spatial derivative occurs in ∇ · S, which involves Bx . The method for estimating numerical error presented in this section can and should be used to estimate error in multidimensional simulations. Using the numerical solution to compute combinations of time and space averages of the thermal, CM kinetic, and electromagnetic energy conservation equations in selected regions of the domain of a multi-dimensional simulation, especially in regions of large spatial or time derivatives, provides an important estimate of numerical error. Although Figures 26–28 provide an important measure of numerical error, they do not indicate how this error is distributed among the various terms in these equations. Methods for determining this error distribution are needed. Better understanding of how to quantify the physical effects of numerical error and the use of unphysical transport processes such as artificial diffusion in simulations is needed. Artificial diffusion stabilizes a code, but may cause large spurious heating and heat flux, and possibly other unphysical transport effects. Increasing demand for simulations that predict spectra for comparison with observations makes understanding these effects increasingly important. There is a difference between predicting observed spectra, and correctly describing the physical mechanisms that cause the spectra. Even if model-predicted spectra agree in some sense with a particular set of observations, this does not imply the model correctly describes the mechanisms that produce the spectra. It is possible that a model with specified, unrealistic values of parameters that determine various transport processes predicts spectra that agree with observations. Such a model might also omit important transport processes. Better analysis of these issues is needed to properly interpret model predictions. 4.2. Energy Flow and Transformation Figures 29–31 provide measures of the contributions of energy fluxes, and the energy conversion rates Q, Qc , and Fm · V to ∂ /∂t z , ∂(ρV 2 /2)/∂t z , and ∂(B 2 /8π )/∂t z . The flux terms have the form − ∇ · F z , for any flux F, and may be written as (F (0, t) − F (zmax , t))/zmax . Aside from the factor of 1/zmax , this is the difference between the fluxes into and

−6 −7 −8

0

20

40

60

80

100 t (sec)

120

140

160

180

200

Figure 29. Spatial average of the terms in the thermal energy equation vs. time. This shows the contribution of the thermal energy flux, and resistive and compressional heating to the rate of change of the thermal energy density. (A color version of this figure is available in the online journal.)

out of the computational domain at the photosphere and upper boundary. The fluxes at zmax are orders of magnitude smaller than those at z = 0. Figures 29 and 31 show that resistive heating plays a minor role in energy balance. Figure 29 shows that the average of Qc , and the net thermal energy flux − ∇·( V) z = (0, t)Vz (0, t)− (zmax , t)Vz (zmax , t) ∼ (0, t)Vz (0, t), which is the flux through the photosphere, largely change sign as t passes through 100 s. The figure also shows that the main contribution to the maxima (minima) of ∂ /∂t z comes from thermal energy flux (compressive heating) for t < 100 s, and that the converse is true for t > 100 s due to the change of sign. Figure 30 shows that the flux of CM kinetic energy plus that of p = 2 /3 into the system is largely cancelled by the sum of the source terms that describe the transformation of thermal, electromagnetic, and gravitational potential energy into CM kinetic energy. The total flux term and the total source term change phase by ∼π as t passes through 100 s. Figure 31 shows that for t < 100 s the contributions of the Poynting flux and the work done by the magnetic Lorentz force to the rate of change of the magnetic energy are positive and negative, respectively. Write the averaged Poynting equation as ∇ · S ∼ ∂(B 2 /8π )/∂t + Fm · V z , ignoring the relatively small resistive heating rate Q, and recall that −∇ ·S z ∼ S(0, t)/zmax . Then for t < 100 s the figure shows that as electromagnetic energy flows into the plasma through the photosphere, most of it appears as an increase in magnetic energy, and the remainder is converted into CM kinetic energy. For t > 100 s, −Fm · V z is mostly positive, indicating that CM kinetic energy is being converted into magnetic energy. The Poynting flux through the photosphere is mainly positive for t > 100 s. 5. A RELATION BETWEEN Me Mi AND β Additional insight into the meaning of the magnetization factor Me Mi is obtained by using FAL data to show that Me Mi ∼ nH 0 /(βnH i ) for 0 km  z  103 km. Here

No. 1, 2010

MHD SHOCK WAVE HEATING OF THE CHROMOSPHERE

285

7 〈 ∂(ρ V2/2)/∂ t 〉 z

6

2

〈 − ∇⋅(p V+ VρV /2) 〉

z

〈p∇⋅ V+ V⋅ F +ρ V⋅ g 〉

5

m

z

4 3

1

−3

ergs−cm −sec

−1

2

0 −1 −2 −3 −4 −5 −6

0

20

40

60

80

100 t (sec)

120

140

160

180

200

Figure 30. Spatial average of the terms in the CM kinetic energy equation vs. time. This shows the contribution of the modified flux, which is the sum of pV and the CM kinetic energy flux, and the rate at which work is done on a fluid element by compressional, magnetic, and gravitational forces to the rate of change of the CM kinetic energy density. The term ∇ · (pV) may be combined with the compressional work term (−Qc ) to yield the pressure gradient work term −V · ∇p. The terms are separated to show that −Qc balances the similar term in the thermal energy equation. (A color version of this figure is available in the online journal.) 2.5 2

〈 ∂(B /8π)/∂ t 〉

2 1.5

ergs−cm−3−sec−1

z

〈 −∇⋅ S 〉 z 〈 − Fm⋅ V 〉 z 102〈 −Q 〉

1

z

0.5 0 −0.5 −1 −1.5 −2

0

20

40

60

80

100 t (sec)

120

140

160

180

200

Figure 31. Spatial average of the terms in Poynting’s equation vs. time. This shows the contribution of the Poynting flux, resistive heating, and the work done by the CM flow on the electromagnetic field to the rate of change of the magnetic energy density. (A color version of this figure is available in the online journal.)

nH 0 = 2.35 × 1015 cm−3 , and β = 8πp/B 2 . The density nH 0 is the FAL density at z = 520 km, which is the height of the FAL temperature minimum. This relation holds as long as the electron and ion collision frequencies are determined almost entirely by collisions with H i. Then the variation of Me Mi , and hence of the magnetization-induced resistivity is determined by the variation of β and nH i . It is a property of the solar atmosphere that both these quantities decrease with height. Rosenthal et al. (2002) and Bogdan et al. (2003) use a twodimensional ideal MHD model to demonstrate that ideal MHD waves undergo mode conversion as they pass through the β = 1

surface below which β > 1, and above which β < 1. Since there are no transport effects in ideal MHD, that model isolates the effect of the transition between high and low β regimes in determining the relative importance of pressure gradient and magnetic forces in accelerating the CM flow, as described by the momentum equation. For the initial state used here, Me Mi = 1 and β = 1.5 at z ∼ 575 km, and β = 1 and Me Mi = 5 at z ∼ 695 km. This suggests that Me Mi and β pass through unity at nearly the same height. If this is true, magnetic effects begin to dominate force balance and current dissipation at approximately the same height near the base of the chromosphere. Then a combination of mode conversion of MHD wave-like disturbances, as described by Rosenthal et al. (2002) and Bogdan et al. (2003), and the onset of magnetization-induced resistivity might play an important role in chromospheric heating. If, as the FAL temperature profile indicates, the temperature varies by no more than a factor of ∼2 throughout the photosphere and chromosphere, it is fairly accurate to write Me Mi ∼ (nH 0 /nH i )2 (B/B0 )2 for 0 km  z  103 km. Here B0 = 241.7 G. Both B and nH i are expected to vary by orders of magnitude with height. Their ratio is the main factor that determines the variation of the magnetization-induced resistivity, and resistive heating rate for z  103 km. For z > 103 km, Me Mi = Me Mp ∼ (nH 0 /nH i )(np0 /np ) (B/B0 )2 . Here np is the proton density, and np0 = 1.2 × 1012 cm−3 . FAL predicts np ∼ 1011 cm−3 from z ∼ 103 km up to the base of the TR. Then Me Mp ∼ 12(nH 0 /nH i )(B/B0 )2 in this region. 6. CONCLUSIONS MHD shock waves may play an important role in chromospheric heating through the mechanisms of compressive and resistive heating. The compressive heating rate might exceed

286

GOODMAN & KAZEMINEZHAD

the resistive heating rate in a shock layer by up to two orders of magnitude. Although the resistive heating rate is much smaller than the compressive heating rate in the shock wave simulations presented here, magnetization-induced resistivity plays a role in determining shock layer structure, and hence the compressive heating rate through its action as a diffusive process. Ideal MHD shock waves have zero thickness due to the absence of any diffusive mechanism to balance nonlinear steepening. The resistivity tensor, or equivalently the conductivity tensor, represents diffusive effects in the model used here, and causes the shock waves to have nonzero thickness. The form of Equation (4) suggests that the characteristic spatial scale over which changes occur in a shock layer is roughly ∝ ηP = (1 + Γ)η . This follows from the fact that if this equation is made dimensionless, the characteristic lengthscale Lc = c2 ηP c /(4π Vc ) emerges. Here Vc is a characteristic speed to be specified, and ηP c is a characteristic value of ηP in the shock layer near the shock front where it is a maximum. Choosing ηP c = 5 × 10−11 s, consistent with Figure 2, gives Lc (km) = 0.36/Vc (km s−1 ). Choosing Vc ∼ 4 km s−1 gives Lc ∼ 0.1 km. This length, together with visual inspection of the figures, suggests that a spatial resolution of 0.1 km is necessary for full resolution of the shock layers. Kinetic lengthscales as represented through transport coefficients determine the structure of shock layers, and hence determine the thermal energy generation rates that occur in these layers. The shock waves evolve from smooth waves generated in the photosphere that propagate upward and steepen into shocks a few hundred kilometers above the photosphere. Within this height range Me Mi < 1, so ηP ∼ η . Figure 2 shows that at the height of shock wave formation ηP ∼ 10−11 s. The figure also shows that as the atmosphere becomes magnetized with increasing height, magnetization effects maintain ηP in a range of ∼10−11 –10−10 s throughout the chromosphere. By comparison, η decreases with increasing height to values 1–3 orders of magnitude below this range. Then in the absence of magnetization-induced resistivity, Lc ∼ 0.1–10 m. The structure of a shock wave is expected to vary significantly between weakly and strongly magnetized regions, reflecting the strong variation in the resistivity and other transport tensors between these regions. Important questions are whether magnetization effects contribute to the stability of a shock wave and extend the range over which it propagates and provides significant heating, and whether weakly magnetized shock waves such as acoustic shock waves are dissipated at significantly lower altitudes than strongly magnetized shock waves. When a shock wave first forms, it has a large compressive heating rate although it forms in a weakly magnetized region of the atmosphere. This is indicated by Figures 2 and 5. Then magnetization-induced resistivity is not necessary for the formation of shock waves near the photosphere that have large compressive heating rates. This is consistent with otherwise ideal, radiative-hydrodynamic simulations that show the formation and radiative damping of acoustic shock waves near optical depths of order unity (Ulmschneider et al. 1977; Ulmschneider & Kalkofen 1977). These shock waves are generated by a vertical velocity driver in the photosphere with a period of ∼25–45 s, and form within one ideal gas pressure scale height below the FAL temperature minimum. As discussed in Section 1, it appears there is insufficient power in acoustic waves generated in the photosphere with periods of ∼25–200 s to balance the chromospheric NRL, and acoustic waves with shorter periods

Vol. 708

are radiatively damped out below the temperature minimum. The MHD simulation presented here, which uses a magnetic field driver with a period of 30 s, similarly predicts the height of shock formation to be ∼150 km below the FAL temperature minimum. This is one ideal gas pressure scale height for T ∼ 5000 K. There might be two basic magnetization-induced heating mechanisms that create and maintain the chromosphere. The first is the resistive heating mechanism outlined in Section 1.1. It operates whether or not shock waves are present, and under steady state or time-dependent conditions. According to two-dimensional models, it is only effective in magnetic structures that are horizontally localized in regions with characteristic diameters 102 km. The second mechanism is compressive heating in MHD shock waves with a structure determined by magnetization-induced resistivity, and probably magnetizationinduced viscosity and thermal diffusion (Braginskii 1965; Chapman & Cowling 1970) that become effective above the temperature minimum. A question is whether one of these heating mechanisms dominates in the network, while the other dominates in the internetwork. It is proposed that there is no chromosphere in regions where Γ  1. This condition occurs where Me Mi  1, as in the photosphere, or where (ρn /ρ)2 is sufficiently small, as in the corona. A chromosphere exists in weakly ionized regions where Me Mp  1, and where there is a sufficiently strong driver of electric current. The surface above the photosphere on which Me Mp = 1 defines the base of the chromosphere. The height of this surface varies with time and position over the photosphere. The chromospheric temperature inversion occurs near this surface. The temperature inversion and associated chromospheric heating and emission occur at lower heights in regions of a stronger magnetic field since Me Mp reaches unity at lower heights in regions of a stronger magnetic field. This assumes the presence of a sufficiently strong current generating process. Such processes are expected to occur almost everywhere on the Sun due to continual photospheric dynamics. A chromosphere exists over almost the entire solar surface (Athay & Dere 1990). Just as there is a chromosphere–corona TR between the weakly ionized, strongly magnetized chromosphere and the strongly ionized, strongly magnetized corona, there is a photosphere–chromosphere transition region (PCTR) between the weakly ionized, weakly magnetized photosphere, and the weakly ionized, strongly magnetized chromosphere. The PCTR lies in the height range 400 km  z  103 km. It is within the PCTR that the weakly ionized plasma makes the transition from weakly to strongly magnetized, and from isotropic resistive heating by electron current dissipation to anisotropic resistive heating by proton Pedersen current dissipation. The transition of weakly ionized regions of the atmosphere from weakly to strongly magnetized is the key to the onset and persistence of chromospheric heating. Tritschler et al. (2008) present observational evidence for heating in a current sheet in the chromosphere above a sunspot umbra. The heating is observed to occur in a height range with upper and lower bounds of 400 km and 800 km above the photosphere. This is in the height range of the PCTR. Those authors point out that the observations are consistent with the basic chromospheric resistive heating mechanism described in Goodman (2004a), outlined here in Section 1.1. These observations may be an example of the current generating

No. 1, 2010

MHD SHOCK WAVE HEATING OF THE CHROMOSPHERE

process in a current sheet heating the chromosphere through this mechanism. The magnetization factor Me Mi varies with time due to the time dependence of the magnetic field strength, particle densities, and temperature. If it decreases from 1 to 1, the heating rate due to Pedersen current dissipation is expected to decrease by orders of magnitude. This may give rise to regions with an intermittent chromosphere. Such a region varies with time between being strongly magnetized and strongly heated, and being weakly magnetized and weakly heated. Observations  at 1 resolution indicate the presence of relatively cool regions in the chromospheric internetwork (Rezaei et al. 2008). It is conjectured that Me Mi  1 in these regions. The surface on which Me Mi = 1 is the interface between the lower atmosphere where transport processes are isotropic and the upper atmosphere where transport processes are anisotropic. Determining if this surface exists in the weakly ionized region of the atmosphere, and mapping it out in space and time should be a high priority for observations given the important implications it has for chromospheric heating. Simultaneous measurements of magnetic field strength, temperature, particle densities, and NRL in the photosphere and chromosphere allow this surface to be detected, and correlated with NRL. These measurements allow the simultaneous mapping out of the β = 1 surface which is also expected to play an important role in atmospheric heating. The authors thank the three referees for important and insightful remarks that improved the manuscript. This work was supported by grant ATM-0650443 from the Solar–Terrestrial Physics Program of the National Science Foundation to the West Virginia High Technology Consortium Foundation. This work made use of NASA’s Astrophysics Data System. REFERENCES Akasofu, S.-I., & de Jager, C. (ed.) 2001, Space Sci. Rev., 95 (Special issue: Challenge to Long-standing Unsolved Space Physics Problems in the 20th Century) Anderson, L. S., & Athay, R. G. 1989, ApJ, 346, 1010 Athay, R. G., & Dere, D. P. 1990, ApJ, 358, 710 Avrett, E. H., & Loeser, R. 2008, ApJS, 175, 229 Bogdan, T. J., et al. 2003, ApJ, 599, 626 Braginskii, S. I. 1965, in Transport Processes in a Plasma, Reviews of Plasma Physics, Vol. 1, ed. M. A. Leontovich (New York: Consultants Bureau), 205 Carlsson, M., Judge, P. G., & Wilhelm, K. 1997, ApJ, 486, L63 Carlsson, M., & Stein, R. F. 1992, ApJ, 397, L59 Carlsson, M., & Stein, R. F. 1994, in Chromospheric Dynamics, ed. M. Carlsson (Norway: Institute of Theoretical Astrophysics, Univ. Oslo), 47 Carlsson, M., & Stein, R. F. 1995, ApJ, 440, L29 Carlsson, M., & Stein, R. F. 1997, ApJ, 481, 500 Carlsson, M., & Stein, R. F. 2002, ApJ, 572, 626 Carlsson, M., et al. 2007, PASJ, 59, S663 Cauzzi, G., et al. 2008, A&A, 480, 515 Chapman, S., & Cowling, T. J. 1970, The Mathematical Theory of Non-uniform Gases (Cambridge: Cambridge Univ. Press) Cowling, T. G. 1957, Magnetohydrodynamics (New York: Interscience) Cuntz, M., Rammacher, W., & Musielak, Z. E. 2007, ApJ, 657, L57 de Wign, A. G., Stenflo, J. O., Solanki, S. K., & Tsuneta, S. 2009, Space Sci. Rev., 144, 275 Dom´ınguez Cerde˜na, I., Kneer, F., & S´anchez Almeida, J. 2003, ApJ, 582, L55 Fontenla, J. M., Avrett, E. H., & Loeser, R. 2002, ApJ, 572, 636 (FAL) Fossum, A., & Carlsson, M. 2005a, ApJ, 625, 556

287

Fossum, A., & Carlsson, M. 2005b, Nature, 435, 919 Fossum, A., & Carlsson, M. 2006, ApJ, 646, 579 Freidberg, J. P. 1982, Rev. Mod. Phys., 54, 801 Goodman, M. L. 1998, ApJ, 503, 938 Goodman, M. L. 2000, ApJ, 533, 501 Goodman, M. L. 2001, Space Sci. Rev., 95, 79 Goodman, M. L. 2004a, A&A, 416, 1159 Goodman, M. L. 2004b, A&A, 424, 691 Goodman, M. L. 2005, ApJ, 632, 1168 Goodman, M. L. 2008, ‘MHD Model Estimates of the Contribution of Driven, Linear, Non-plane Wave Dissipation to Chromospheric Heating Using a Complete Electrical Conductivity Tensor, American Geophysical Union Meeting, December 2008, Session SH51C, Abstract No. 3323. Jeffery, A. 1966, Magnetohydrodynamics (Edinburgh: Oliver & Boyd) Judge, P. G. 2008, ApJ, 683, L87 Judge, P. G., Carlsson, M., & Stein, R. F. 2003, ApJ, 597, 1158 Judge, P. G., Carlsson, M., & Wilhelm, K. 1997, ApJ, 490, L195 Judge, P. G., & Carpenter, K. G. 1998, ApJ, 494, 828 Judge, P. G., & Centeno, R. 2008, ApJ, 687, 1388 Kazeminezhad, F., & Goodman, M. L. 2006, ApJS, 166, 613 (KG06) Khodachenko, M. L., Arber, T. D., Rucker, H. O., & Hanslmeier, A. 2004, A&A, 422, 1073 Khodachenko, M. L., Arber, T. D., Rucker, H. O., & Hanslmeier, A. 2006, Adv. Space Res., 37, 447 Leake, J. E., Arber, T. D., & Khodachenko, M. L. 2005, A&A, 442, 1091 Lites, B. W., et al. 2008, ApJ, 672, 1237 Mitchner, M., & Kruger, C. H. 1973, Partially Ionized Gases (New York: Wiley) Nagata, S., et al. 2008, ApJ, 677, L145 Orozco Su´arez, D., et al. 2007, ApJ, 670, L61 Penn, M. J., & Livingston, W. 2006, ApJ, 649, L45 Rezaei, R., Bruls, J. H. M. J., Schmidt, W., Beck, C., Kalkofen, W., & Schlichenmaier, R. 2008, A&A, 484, 503 Rosenthal, Bogdan, T. J., Carlsson, M., Dorch, S. B. F., Hansteen, V., McIntosh, S. W., McMurray, A., Nordlund, A., & Stein, R. F. 2002, ApJ, 564, 508 Rybicki, G. B., & Hummer, D. G. 1991, A&A, 245, 171 Rybicki, G. B., & Hummer, D. G. 1992, A&A, 262, 209 S´anchez Almeida, J. 1997, ApJ, 491, 993 S´anchez Almeida, J. 2000, ApJ, 544, 1135 S´anchez Almeida, J. 2005, A&A, 438, 727 S´anchez Almeida, J. 2006, A&A, 450, 1199 S´anchez Almeida, J., Dom´ınguez Cerde˜na, I., & Kneer, F. 2003, ApJ, 597, L177 S´anchez Almeida, J., & Landi degl’Innocenti, E. 1996, Sol. Phys., 164, 203 S´anchez Almeida, J., Landi degl’Innocenti, E., Martinez Pillet, V., & Lites, B. W. 1996, ApJ, 466, 537 S´anchez Almeida, J., & Lites, B. W. 2000, ApJ, 532, 1215 S´anchez Almeida, J., M´arquez, I., Bonet, J. A., Dom´ınguez Cerde˜na, I., & Muller, R. 2004, ApJ, 609, L91 S´anchez Almeida, J., Ramos, A. Asensio, Bueno, J. Trujillo, & Cernicharo, J. 2001, ApJ, 555, 978 S´anchez Almeida, J., Viticchi´e, B., Landi Degl’Innocenti, E., & Berrilli, F. 2008, ApJ, 675, 906 Schrijver, C. J., & Zwaan, C. 2000, Solar and Stellar Magnetic Activity (Cambridge: Cambridge Univ. Press) Socas-Navarro, H. 2005, ApJ, 633, L57 Socas-Navarro, H., & Lites, B. W. 2004, ApJ, 616, 587 Socas-Navarro, H., & S´anchez Almeida, J. 2002, ApJ, 565, 1323 Socas-Navarro, H., & S´anchez Almeida, J. 2003, ApJ, 593, 581 Solanki, S. K. 1993, Space Sci. Rev., 63, 1 Solanki, S. K., Bruls, J. H. M. J., Steiner, O., Ayres, T., Livingston, W., & Uitenbroek, H. 1994, in NATO ASI Series C433, Solar Surface Magnetism, ed. R. J. Rutten & C. J. Schrijver (Dordrecht: Kluwer), 91 Tritschler, A., Uitenbroek, H., & Reardon, K. 2008, ApJ, 686, L45 Tu, C-Y., Zhou, C., Marsch, E., Wilhelm, K., Zhao, L., Xia, L-D., & Wang, J-X. 2005, ApJ, 624, L133 Uitenbroek, H. 2001, ApJ, 557, 389 Ulmschneider, P., & Kalkofen, W. 1977, A&A, 57, 199 Ulmschneider, P., Kalkofen, W., Nowak, T., & Bohn, U. 1977, A&A, 54, 61 Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635 Whitham, G. B. 1974, Linear and Nonlinear Waves (New York: Wiley) Withbroe, G. L., & Noyes, R. W. 1977, ARA&A, 15, 363