THE HYDROMAGNETIC INTERIOR OF A SOLAR QUIESCENT PROMINENCE. II. MAGNETIC DISCONTINUITIES AND CROSS-FIELD MASS TRANSPORT

We adopt the following one-fluid hydromagnetic description of the corona and prominence:

$$\begin{equation} \rho \left[{\partial {\bf v} \over \partial t} + \left({\bf v} \cdot \nabla \right) {\bf v} \right] = {1 \over 4 \pi }(\nabla \times {\bf B}) \times {\bf B} -\nabla p -\rho g {\hat{z}}, \end{equation} \tag{ 1 }$$

$$\begin{equation} {\partial {\bf B} \over \partial t} - \eta \nabla ^2 {\bf B} = \nabla \times ({\bf v} \times {\bf B} ), \end{equation} \tag{ 2 }$$

$$\begin{equation} {\partial \rho \over \partial t} + \nabla \cdot (\rho {\bf v} ) = 0, \end{equation} \tag{ 3 }$$

where v, B, ρ, p, and g denote, respectively, the velocity, magnetic field, density, pressure, and uniform gravitational acceleration in the negative Cartesian z-direction. We also need an equation of state and another for energy transport, the latter involving anisotropic thermal conduction, radiative processes, and volumetric heating operating in the corona. If these energy processes are specified explicitly in terms of the fluid and field variables, Equations (1)–(3) together with the scalar equations of state and energy transport are a complete set to determine these variables, including the temperature T.

For simplicity in our analysis of basic effects, we neglect viscosity in momentum equation (1) and take electrical conductivity to be isotropic, characterized by a constant σ related to the resistivity η = c²/4πσ, where c is the speed of light. The length scale l₀ ≈ 7 × 10⁷ cm is of the order of spatial resolutions achievable in solar observations from the ground. In the corona and prominence interior, the resistive diffusion of a magnetic field of that length scale has a timescale τ_D = l²₀/η of tens of years and longer if σ is estimated with the formula of Spitzer (1962) in the temperature range of 10⁴–10⁶ K for a fully ionized gas. Therefore, for dynamical events in the corona and prominences observed over timescales of hours at spatial resolutions not drastically smaller than l₀, the frozen-in condition is an excellent approximation with an important qualification. There can be processes, unresolved in such observations, creating field structures of sufficiently small scales over which resistive dissipation is significant. The observable large-scale consequences of this dissipation are not possible to be described under the frozen-in condition. For example, at the large scales of resolution-limited observations, a significant larger "effective" resistivity may appear to be operating. This general physical point important for understanding quiescent-prominence observations is a central motivation of this series of papers.

Consider the relaxation of a fluid to static equilibrium under the frozen-in condition, setting η = 0. This equilibrium is described by

$$\begin{equation} {1 \over 4 \pi }(\nabla \times {\bf B}) \times {\bf B} -\nabla p -\rho g {\hat{z}} = 0, \end{equation} \tag{ 4 }$$

$$\begin{equation} \nabla \cdot {\bf B} = 0, \end{equation} \tag{ 5 }$$

$$\begin{equation} \nabla \cdot \left[ \kappa { \left({\bf B} \cdot \nabla T \right) {\bf B} \over |{\bf B}|^2 } \right] - {\mathcal L} + {\mathcal H} = 0, \end{equation} \tag{ 6 }$$

$$\begin{equation} p = {k_B \over m_0} \rho T. \end{equation} \tag{ 7 }$$

In addition to force balance, we have thermal balance among a radiative-loss sink ${\mathcal L}$, a heating source ${\mathcal H}$, and a field-aligned thermal flux with conductivity κ, all these quantities assumed to be expressible as known functions of the fluid and field. We assume the ideal gas law, k_B and m₀ being the Boltzmann constant and the mean molecular mass, respectively.

Consider a 2D Cartesian system with the solenoidal field

$$\begin{equation} {\bf B} = \left[ Q, {\partial A \over \partial z}, - {\partial A \over \partial y} \right], \end{equation} \tag{ 8 }$$

in terms of the flux function A constant along each field line and an x-component Q, x being an ignorable coordinate. The electric current density J = (c/4π)∇ × B is given by

$$\begin{equation} \nabla \times {\bf B} = \left[ -\nabla ^2 A, {\partial Q \over \partial z}, - {\partial Q \over \partial y} \right]. \end{equation} \tag{ 9 }$$

The equilibrium Lorentz force cannot have an x-component, requiring Q(y, z) = Q[A(y, z)] so that Equation (4) reduces to the two components in the y − z plane,

$$\begin{equation} \nabla ^2 A + Q(A)Q^{\prime }(A) + 4 \pi {\partial p(A, z) \over \partial A} = 0, \end{equation} \tag{ 10 }$$

$$\begin{equation} {\partial p(A, z) \over \partial z} + \rho g = 0, \end{equation} \tag{ 11 }$$

for force balance across and along the field, respectively (Low 1975).

The condition Q = Q(A), together with Equation (10), implies that in equilibrium the current density takes the form

$$\begin{equation} \nabla \times {\bf B} = 4 \pi {\partial p(A, z) \over \partial A} {\hat{x}} + Q^{\prime }(A) {\bf B}, \end{equation} \tag{ 12 }$$

a field-aligned component with a proportionality Q'(A) constant on constant-A field lines and a cross-field component in the x-direction that generates the Lorentz force to balance the forces of pressure and weight. If the latter forces are absent, setting p ≡ 0, we have a force-free field.

In a stratified atmosphere, we express p(y, z) = p[A(y, z), z], to describe its hydrostatic variation with height z along a constant-A field line. Each point (y, z) on the y − z plane is identified by the value A = A₀, a constant, at that point and a height z-varying along the curve A = A₀, the way to visualize the partial derivatives of p with respect to (A, z) as independent variables. Henceforth, the dependence of a scalar function on (A, z) shall be indicated when we mean partial derivatives with respect to these variables. The hydrostatic profile of p is subject to the balance of forces across the field lines, the Grad–Shrafanov equation (10).

Denote the field-aligned thermal flux in Equation (6) by ${\bf F}={\mathcal F}{\bf B}$, where

$$\begin{equation} {\mathcal F}= {\kappa \over \left(|\nabla A|^2 + Q^2 \right)^{1/2}} {\partial (T, A) \over \partial (y, z)}, \end{equation} \tag{ 13 }$$

introducing the Jacobian operator. Then energy-balance equation (6) takes the form

$$\begin{equation} {\partial ({\mathcal F}, A) \over \partial (y, z)} - {\mathcal L} + {\mathcal H} = 0. \end{equation} \tag{ 14 }$$

Equations (7), (10), (11), and (14) constitute a complete system of one algebraic and three partial differential equations (PDEs) determining the four scalar unknowns A, p, ρ, T.

To see the structure of this system of equations, introduce a constant normalization T₀, setting T = T₀θ and defining the constant hydrostatic scale height Λ₀ = k_BT₀/m₀g, using the notations in Paper I. First eliminate ρ between Equations (7) and (11) and then integrate with respect to z, treating (A, z) as the independent variables, to obtain

$$\begin{equation} p(A, z) = p_0(A) \exp \left[ - \int _0^ z {dz^{\prime } \over \Lambda _0 \theta (A, z^{\prime })} \right], \end{equation} \tag{ 15 }$$

where p₀(A) arises as an arbitrary "constant" in that integration with respect to z. In the representation θ(y, z) = θ(A, z), Equations (13) and (14) take the forms

$$\begin{equation} {\mathcal F}= {\kappa T_0 \over (|\nabla A|^2 + Q^2 )^{1/2}} {\partial \theta (A, z) \over \partial z}, \end{equation} \tag{ 16 }$$

$$\begin{equation} {\partial {\mathcal F}(A, z) \over \partial z} - {\mathcal L} + {\mathcal H} = 0, \end{equation} \tag{ 17 }$$

the latter effectively being an ordinary differential equation balancing a field-aligned thermal conduction against energy loss and gain along each constant-A field line. By using Equations (7) and (15) to express p and ρ in terms of θ, the problem is then reduced to solving Equations (10) and (17) for A and θ as the two coupled unknowns subject to suitable boundary conditions. This is a 2D extension of the 1D KS slab in Paper I.

Taken as a boundary-value problem in PDEs, there are two types of free functions to prescribe in order to render Equations (10) and (17) specific. We must prescribe the thermal conductivity κ, energy sink ${\mathcal L}$, and heat source ${\mathcal H}$ as explicit functions of the fluid and field variables, to identify a specific fluid by its properties. The thermally balanced equilibrium states of this fluid are then generated by all the admissible forms of Q(A) and p₀(A), usually called the generating functions. For each pair of analytic [Q(A), p₀(A)], the boundary-value problem in the tradition of classical analysis gives either no solution, because the governing PDEs are nonlinear, or an everywhere-analytic solution (Courant & Hilbert 1962). The principal result of Paper I is that the set of all everywhere-analytic solutions is a subset of measure zero of the physically admissible thermally balanced equilibrium states. So the same result may be expected in the 2D extension of the model.

If the electrical resistivity η ≠ 0, there is no steady state because the current density would then decay in time. If η = 0, the fluid surface identified as a flux surface of B at any one time remains a flux surface for all time under continuous displacement, the frozen-in condition expressed by the induction equation

$$\begin{equation} {\partial {\bf B} \over \partial t} = \nabla \times \left({\bf v} \times {\bf B} \right). \end{equation} \tag{ 18 }$$

This equation implies, in particular,

$$\begin{equation} {\partial A \over \partial t} + {\bf v} \cdot \nabla A = 0, \end{equation} \tag{ 19 }$$

for our 2D system. It is easy to show that as long as v is continuous, setting A(y, z, t) = A₀, a fixed constant, describes an evolving fluid surface. Therefore, the total mass M(A₁, A₂) distributed with x as an ignorable coordinate between any two flux surfaces A = A₁ and A = A₂, A₁ and A₂ being constants, is a constant in time. More generally, induction equation (18) states that the net magnetic flux across any fluid surface area is conserved as it deforms continuously. Therefore, the total magnetic flux Φ(A₁, A₂) contributed by B_x = Q across the area in the y − z plane between two field lines, A = A₁ and A = A₂, is also a constant in time. When given one of the above 2D thermally balanced equilibrium states, the system must be thought of physically as having arrived at this state with a specific M(A₁, A₂) and Φ(A₁, A₂) fixed for the system for all time. If we specify M(A₁, A₂) and Φ(A₁, A₂), we are asking for the equilibrium state for a completely identified fluid in terms of its mass and x-directed flux distributed over the constant-A fluid surfaces in some specific manner. For that equilibrium, the functions p₀(A) and Q(A) must be given forms that recover M(A₁, A₂) and Φ(A₁, A₂), respectively, in the solution. That is, p₀(A) and Q(A) are unknowns in this physical formulation of the equilibrium expressed by integral equations that relate them to M(A₁, A₂) and Φ(A₁, A₂), respectively; see the discussion of this integro-PDE system by Wu & Low (1987).

Conceptually, M(A₁, A₂) and Φ(A₁, A₂) are physical properties that can be freely prescribed. It is thus interesting that, in general, not all prescribed pairs yield an everywhere-analytic solution to their integro-PDE system. The unavoidable singularities in the non-analytic solutions are the weak solutions of the boundary-value problem, meaning that these singularities are not arbitrary but subject to the same physical laws as the integral equivalence to the PDEs governing the continuous part of the global solutions (Courant & Hilbert 1962). From this analysis, an exhaustive search of everywhere-analytic solutions to the boundary-value problems posed by PDEs (7), (10), (11), and (14) with an a priori prescription of continuous p₀(A) and Q(A) would miss out on the singular solutions that populate the total set of physically admissible equilibrium states.

Coupling the Grad–Shrafanov equation (10) to the energy-transport equation (17) is an analytically intractable problem. The basic property of interest here is in equilibrium states with a temperature T that is discontinuous across a magnetic flux surface. This effect can be investigated with the simpler equilibrium states in which θ(y, z) = θ(A) with ${\mathcal L} = {\mathcal H} = 0$. Then ${\mathcal F} = 0$ by Equation (16), and energy equation (17) is trivially satisfied. Equation (15) gives

$$\begin{equation} p(A, z) = p_0(A) \exp \left[ - {z \over \Lambda _0 \theta (A)} \right], \end{equation} \tag{ 20 }$$

reducing the problem to one of solving the Grad–Shrafanov PDE

$$\begin{equation} \nabla ^2 A + Q(A)Q^{\prime }(A) + 4 \pi {\partial \over \partial A} \left(p_0(A) \exp \left[ - {z \over \Lambda _0 \theta (A)} \right]\right) = 0, \end{equation} \tag{ 21 }$$

for A with prescribed functions p₀(A), Q(A), and θ(A), subject to boundary conditions on A. Since thermal conduction is absent across the field lines, θ(A) is not required to be a continuous function of A. Moreover, the pressure is generally not continuous when the temperature is discontinuous, that is, the function p₀(A) may also be discontinuous. The derivatives of p₀(A) and θ(A) on the right-hand side of Equation (21) contribute δ-functions to ∇²A and, hence, to the current density given by Equation (9).

All the above magnetostatic problems, as either a purely PDE system or an integro-PDE system, are generally intractable. So we shall make our physical points in Section 4 by specialized solutions to give insights into the general solutions. Before proceeding, we need to address in the next section an important assumption in our analysis. In a real plasma, not only is the electrical conductivity not infinite, that is, there is a small but nonzero resistivity, but its cross-field thermal conductivity is also small but nonzero. The non-dimensional ratio of the resistivity to that thermal diffusion of heat determines whether the breakdown of the frozen-in condition as an approximation precedes the breakdown of the assumption of zero cross-field thermal conduction.

We adopt the isotropic, temperature-dependent Spitzer (1962) resistivity,

$$\begin{equation} \eta = 5 \times 10^{12} T^{-3/2}\, \mathrm{cm}^2\,\mathrm{s}^{-1}, \end{equation} \tag{ 22 }$$

neglecting for simplicity a weak anisotropy owing to the presence of a magnetic field. We also neglect other complications like the thermoelectric effect on the resistivity arising from an extreme thermal gradient across the field.

We obtain the values η = 5 × 10⁶ cm² s⁻¹ at T = 10⁴ K in a prominence and η = 5 × 10³ cm² s⁻¹ at T = 10⁶ K in the corona. The resistive diffusion time τ_d = L²/η of a magnetic structure of size L = 200 km, representing a lower-limit length scale observable with Hinode/SOT, for example, is of the order of 10⁸ s and 10¹¹ s for the prominence and corona, respectively. These times are too long to be relevant at the timescale of 100 s characteristic of the prominence's vertical threads. Turning the question around, the length scales for which τ_d ≈ 100 s are L_P = 2 × 10⁻¹ km in the prominence and L_C = 10⁻² km in the corona, both three orders of magnitude, at least, smaller than their respective hydrostatic scale heights of Λ_P = 300 km and Λ_C = 3 × 10⁴ km.

The mean free paths of the electrons and protons for like-particle collisions have the same dependence on T²/N, where N is the number density of the species. Let us simplify by assuming the same temperature for the two species. Then, this common mean free path is of the order of l_P = 5 cm in a prominence plasma of density N ≈ 10¹¹ cm⁻³ and T = 10⁴ K and l_C = 50 km in a coronal plasma of density N ≈ 10⁹ cm⁻³ and T = 10⁶ K. Therefore, when resistive diffusion in the prominence becomes relevant at a scale of L_P = 2 × 10⁻¹ km, this scale is still much larger than the mean free path of l_P = 5 cm. On the other hand, because of its high electrical conductivity, resistive diffusion becomes relevant in the corona for scales of L_C = 10⁻² km or smaller, well below its mean free path of l_C = 50 km. In other words, as a field gradient steepens toward a tangential discontinuity under the frozen-in condition, it proceeds until the fluid description breaks down, locally at where the tangential discontinuity is developing without arriving at a resistive fluid regime. The situation is actually more optimistic because in the presence of a 10 G magnetic field, the ions are tied to the field with a cyclotron radius of about 90 cm for protons at T = 10⁶ K, so that the fluid picture may be said to break down only along the field. A proper treatment of this breakdown is outside the scope of the present paper, but it seems reasonable to assume that such a developing tangential discontinuity will excite a host of plasma-kinetic effects to dissipate its growing current density while the fluid description and the frozen-in condition break down. Note that the fluid picture is approximately good along the field to some relatively small scales in the low corona, considering that the mean free path in the corona is of the order of 10⁻² of the coronal hydrostatic scale height Λ_C.

We adopt the Spitzer thermal conductivity and thermometric diffusivity

$$\begin{eqnarray} \kappa _{\parallel }(T)&=&2 \times 10^{-6} T^{5/2} \,\mathrm{erg} \,\mathrm{cm}^{-1} \,\mathrm{s}^{-1} \,\mathrm{K}^{-1},\nonumber \\ K_{\parallel }(T)&=& \kappa (T)/(5Nk_B) \nonumber \\ &=&3 \times 10^9 T^{5/2} N^{-1} \,\mathrm{cm}^2 \,\mathrm{s}^{-1}, \end{eqnarray} \tag{ 23 }$$

respectively, taken from Parker (1994), using a subscript to indicate thermal conduction along the magnetic field. The electrons as the carriers of thermally conducted heat are tightly tied to the magnetic field lines. Thermal conduction across the field has the coefficients

$$\begin{equation} \kappa _{\perp }=\kappa _{\parallel }/(\omega _e \tau _e)^2, \end{equation} \tag{ 24 }$$

$$\begin{equation} K_{\perp }=K_{\parallel }/(\omega _e \tau _e)^2, \end{equation} \tag{ 25 }$$

greatly reduced from the respective parallel coefficients by a factor of (ω_eτ_e)⁻² ≪ 1, where

$$\begin{equation} \omega _e=eB/m_e c \ \mathrm{rad} \ \mathrm{s}^{-1}, \end{equation} \tag{ 26 }$$

$$\begin{equation} \tau _e={11.4 m_e^{1/2} \over {\rm log} \Lambda } {T^{3/2} \over N} \ \mathrm{s} \end{equation} \tag{ 27 }$$

are the electron cyclotron frequency and the time for an electron to be significantly deflected by Coulomb interaction, respectively, in standard notation. For our purpose we set the Coulomb factor log Λ ≈ 10 for application. For example, if the field strength B = 100 G and the electron density is N ≈ 10¹⁰ cm⁻³ at T = 10⁶ K, the cross-field thermal conduction is reduced by a factor of (ω_eτ)⁻² ≈ 10⁻¹⁰ (Parker 1994). Taking the field to be 10 G, closer to the fields threading quiescent prominences, that reduction factor is still exceedingly small at 10⁻⁸.

The thermometric diffusivity, K_⊥(T), describing the diffusion of temperature, has the same dimension as resistivity η, and we have an interest in their dimensionless ratio. Direct evaluation gives

$$\begin{equation} \omega _e=1.6 \times 10^7 B, \end{equation} \tag{ 28 }$$

$$\begin{equation} \tau _e=2.3 \times 10^{-2} {T^{3/2} \over N}. \end{equation} \tag{ 29 }$$

Then, we have

$$\begin{equation} \omega _e \tau _e\approx 3.73 \times 10^5 {B T^{3/2} \over N}, \end{equation} \tag{ 30 }$$

$$\begin{equation} K_{\perp }=2.2 \times 10^{-2}{N \over B^2 T^{1/2} }. \end{equation} \tag{ 31 }$$

Evaluate the ratio

$$\begin{eqnarray} \epsilon & = & {K_{\perp } \over \eta } \nonumber \\ & = & 4.4 \times 10^{-15} {N T \over B^2} \nonumber \\ & = & 32 {Nk_BT \over B^2}. \end{eqnarray} \tag{ 32 }$$

The quantity p = 2Nk_BT is just the pressure of the fully ionized gas, which defines the plasma beta

$$\begin{equation} \beta = {8 \pi p \over B^2}. \end{equation} \tag{ 33 }$$

Hence, we have the ratio

$$\begin{equation} \epsilon = 0.64 \beta. \end{equation} \tag{ 34 }$$

The field intensity in a prominence is about 10–60 G, from polarimetric observations (Casini et al. 2003; Lopez Ariste & Casini 2003). If we take N = 10¹¹ cm⁻³, T = 10⁴ K, and B = 30 G, we get p = 0.28 erg cm⁻³, β = 7.7 × 10⁻³, and = 5 × 10⁻³.

With ≪ 1, the thinning of a steep gradient layer to zero will cross the threshold thickness for resistive dissipation of the growing current density in the layer, well before cross-field thermal conduction becomes important. If ≫ 1, the thinning steep gradient will excite cross-field thermal diffusion before resistive diffusion. In this case, the thermal balances in the two adjacent tubes are perturbed and a time-dependent evolution is initiated without breaking the frozen-in condition. It is possible that this evolution on its own develops steep gradient layers even thinner to also excite resistive dissipation. This is just as interesting an outcome, but the case of ≪ 1 is simpler, a clear case of the magnetic thermal insulation directly causing a magnetic tangential discontinuity to form at a flux surface, a novel case of the Parker spontaneous current sheets. Note that β ≪ 1 does not mean that the field is so rigid that it cannot change when the plasma pressure changes. Rather, the field being strong can change by just a small amount to accommodate any plasma-pressure change, including an inevitable jump in the plasma pressure.

Consider equilibrium in the presence of a spatially uniform temperature. Set θ ≡ 1, that is, T = T₀, the uniform temperature, and the Grad–Shafranov equation (21) takes the form

$$\begin{equation} \nabla ^2 A + Q(A)Q^{\prime }(A) + 4 \pi {dp_0(A) \over dA} \exp (-k_0 z) = 0, \end{equation} \tag{ 35 }$$

introducing the inverse scale height k₀ = 1/Λ₀. Let us prescribe

$$\begin{equation} p_0(A)= p_1 \exp (\lambda A ), \end{equation} \tag{ 36 }$$

$$\begin{equation} Q(A) = B_1, \end{equation} \tag{ 37 }$$

where p₁, λ, and B₁ are constants. The pressure in Equation (20) takes the form

$$\begin{equation} p(A, z) = p_1 \exp \left[ - k_0 z + \lambda A \right], \end{equation} \tag{ 38 }$$

and the Grad–Shrafanov equation (35) reduces to

$$\begin{equation} \nabla ^2 A + 4 \pi p_1 \lambda \exp \left[ - k_0 z+ \lambda A \right] = 0. \end{equation} \tag{ 39 }$$

The substitution,

$$\begin{equation} A = {k_0 \over \lambda } z + \phi (y, z), \end{equation} \tag{ 40 }$$

transforms the Grad–Shrafanov equation into the Liouville PDE

$$\begin{equation} \nabla ^2 \phi + 4 \pi p_1 \lambda \exp \left(\lambda \phi \right) = 0, \end{equation} \tag{ 41 }$$

with the complete integral

$$\begin{equation} \exp (\lambda \phi ) = {2 (\nabla u )^2 \over \pi p_1 \lambda ^2 (1 + u^2 + v^2 )^2 }, \end{equation} \tag{ 42 }$$

introducing u and iv as the real and imaginary parts of an arbitrary analytic complex function F(ω) of the complex variable ω = y + iz, that is,

$$\begin{equation} u(y, z) + {i} v(y,z) = F (\omega). \end{equation} \tag{ 43 }$$

Each prescription of F(ω) yields a Liouville solution ϕ defining a static equilibrium. The field, with B_x = B₁, is given by the flux function (40), frozen into an isothermal density distribution

$$\begin{equation} \rho g = k_0 p_1 \exp \left(\lambda \phi \right), \end{equation} \tag{ 44 }$$

using the ideal gas law and the pressure (38). The field and density are expressed as functions of space through ϕ(y, z).

For our purpose in this paper, we concentrate on the family of solutions generated by

$$\begin{eqnarray} F(\omega) & = & \exp \left(k \omega \right) + a, \nonumber \\ u(y, z) & = & \exp \left(k y \right) \cos kz + a, \\ v(y, z) & = & \exp \left(k y \right) \sin kz,\nonumber \end{eqnarray} \tag{ 45 }$$

where k and a are two arbitrary constants. Equation (42) gives the solution

$$\begin{equation} \exp (\lambda \phi ) = {2 k^2 \exp (2 k y ) \over \pi p_1 \lambda ^2 (1 + a^2 + 2 a \exp (k y ) \cos kz + \exp (2 k y ) )^2 }. \end{equation} \tag{ 46 }$$

For a given a, define $\exp \left(k y_1 \right) = \sqrt{1 + a^2}$. Then a transformation y → y + y₁ redefines the function ϕ(y, z):

$$\begin{equation} \exp (\lambda \phi ) = {2 k^2 \over \pi p_1 \lambda ^2 (\sqrt{1 + a^2} \cosh ky + a \cos kz )^2 }. \end{equation} \tag{ 47 }$$

This transformation has moved the origin along the y-axis to a location about which the redefined function ϕ(y, z) is symmetric in y. The free parameter k is the wavenumber of this function's periodicity in z, as well as its exponential decline with y in either direction.

To relate to the physics we are interested in, we prescribe the uniform temperature T₀, fixing the inverse scale height k₀. Take the free parameter λ to have the dimension of the reciprocal of the field intensity times length. We may then set λ = k₀/B₂, introducing a constant field intensity B₂ to replace λ as a free parameter. We now introduce a third field intensity B₃ = (k/k₀)B₂ to represent the free parameter k as may be convenient to do depending on the physical discussion. Note that λ = k/B₃ by definition. Then the magnetic field is given by

$$\begin{equation} A=B_2 z - 2\left({B_3 \over k}\right) \log (\sqrt{1 + a^2}\cosh ky + a \cos kz ), \end{equation} \tag{ 48 }$$

$$\begin{eqnarray} {\bf B}&=&\left[ B_1, B_2 + {2 a B_3 \sin kz \over \sqrt{1 + a^2} \cosh ky + a \cos kz},\right.\nonumber\\ &&\left. {2 B_3 \sqrt{1 + a^2} \sinh ky \over \sqrt{1 + a^2} \cosh ky + a\cos kz} \right], \end{eqnarray} \tag{ 49 }$$

in which (B₁, B₂, B₃) are three field intensities we may arbitrarily prescribe with k = (B₃/B₂)k₀. This is a four-parameter family of everywhere continuous solutions, the fourth being the dimensionless constant a.

Each of these fields is in equilibrium with the isothermal density

$$\begin{equation} \rho g = {B_2 B_3 k \over 2 \pi (\sqrt{1 + a^2} \cosh ky + a\cos kz )^2}, \end{equation} \tag{ 50 }$$

which is periodic in z, confined about the central plane y = 0 with an exponential decay over a scale k⁻¹ to zero as y → ±∞. The density is independent of B_x = B₁. From force-balance equation (4) with x as an ignorable coordinate, the uniform B_x produces only a uniform magnetic pressure that is immaterial to the force balance. This component plays a more explicit role through the anisotropic influence of the field on thermal conduction in energy equation (6).

4.1.1. The 1D Isothermal KS Slab

Set a = 0 and we have the 1D-slab isothermal solution of Kippenhahn and Schlüter (1957),

$$\begin{equation} A=B_2 z - 2\left({B_3 \over k}\right) \log (\cosh ky ), \end{equation} \tag{ 51 }$$

$$\begin{equation} {\bf B}=[ B_1, B_2, 2 B_3 \tanh k y ], \end{equation} \tag{ 52 }$$

$$\begin{equation} \rho g = {B_2 B_3 k \over 2 \pi }{\rm sech}^2 k y. \end{equation} \tag{ 53 }$$

Thus, setting a ≠ 0 generalizes this strictly y-dependent solution to vary with both y and z, keeping x as an ignorable coordinate. In Section 4.3 we will be examining a different generalization to 2D, allowing variations in y and x while keeping z as an ignorable coordinate. Therefore, it is useful to see clearly the following specific property of this 1D KS solution as a reference state for these two generalizations.

Equation (52) shows that the equilibrium field is the result of weighing down a uniform background field ${\bf B}_{{\rm bg}}=B_1{\hat{x}}+B_2{\hat{y}}$ into a sheared, bowed shape. By the 1D nature of the equilibrium, a measure of the total weight is given by

$$\begin{eqnarray} W &=& \int _{-\infty }^{\infty } \rho g dy \nonumber \\ &=& {B_2 B_3 \over \pi }. \end{eqnarray} \tag{ 54 }$$

This is the weight in an infinite column of unit cross-sectional area aligned in the y-direction. In the y − z plane, W may be called the integrated weight per unit length in the z-direction, which we henceforth refer to as the weight for brevity. The greater the weight, the greater is the sag defined by B_z = ±2B₃ at y → ±∞ for a fixed B_y = B₂. If we hold B_bg fixed, it is useful to replace B₃ = πW/B₂ with W as a free parameter. Then the wavenumber k = k₀B₃/B₂ can be expressed as k = πk₀W/B²₂, and Equations (51)–(53) read

$$\begin{equation} A=B_2 \left(z - 2\Lambda _0 \log \left[ \cosh \left({\pi W \over B_2^2}{y \over \Lambda _0} \right) \right]\right), \end{equation} \tag{ 55 }$$

$$\begin{equation} {\bf B}= B_1 {\hat{x}} + B_2 \left[ 0, 1, 2 {\pi W \over B_2^2} \tanh \left({\pi W \over B_2^2}{y \over \Lambda _0} \right)\right], \end{equation} \tag{ 56 }$$

$$\begin{equation} \rho g = \frac{\pi W^2}{2 \Lambda _0 B_2^2} {\rm sech}^2 \left({\pi W \over B_2^2}{y \over \Lambda _0} \right), \end{equation} \tag{ 57 }$$

where we have used k₀ = 1/Λ₀.

In 3D space, we have an infinite, vertical slab represented by a density ρ(y) symmetrically concentrated about the vertical plane y = 0. The characteristic slab width is B²₂Λ₀/πW. Figures 1 and 2, with a fixed B_y = B₂, suffice for illustrating the parametric dependence of the solution on the isothermal temperature and the total mass W/g in the frozen-in condition.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 1 shows the shapes of the bowed field lines projected on the y − z plane for three slab solutions of increasing values of W, all having the same scale height Λ₀. The greater W is, the more inclined the field is below the horizontal indicated by |B_z/B_y| = B₃/B₂ = πW/B²₂ at y → ±∞. The density distributions show the corresponding parametric trend of a narrowing of slab width and an increasing peak value of the density at y = 0. When the total weight W is greater, it weighs more on the background field into a deeper bow such that the slab in equilibrium has a smaller width.

Figure 2 shows the shapes of the bowed field lines projected on the y − z plane for three slab solutions of increasing values of the isothermal temperature T₀ and its scale height Λ₀, all having the same weight W. Increasing the scale height implies a more gradual spatial decline of stratified density with z along the inclined field. This means that the same total weight W in the three cases is distributed over increasingly thicker slab widths as shown. In all three cases, the field lines tend to the same inclination below the horizontal at y → ±∞ because that inclination is determined by W alone.

4.1.2. The 2D Isothermal KS Slab

If a ≠ 0, we have a 2D KS slab with a slab width k⁻¹ and a periodic vertical structure of wavenumber k. Figures 3 and 4 show four examples with the same values of B₂, B₃ but with increasing values of a as a measure of departure from one-dimensionality. In each figure, the thin lines are constant-A field lines projected on the y − z plane and the thick lines are contours of constant density. The constant field component B_x = B₁ may be taken to be the same for the three examples. This field component does not affect the equilibrium in the y − z plane, but it gives the field, in 3D space, a sheared configuration for the field lines extending to y → ±∞ and a twisted flux-rope configuration for the magnetic islands, or regions of closed field lines, found in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For a < a₀ = B₂/2B₃, all the field lines extend to infinity and are unevenly bowed because the loaded mass per unit magnetic flux varies across the field, the cases in Figure 3. For a > a₀, periodic magnetic islands appear, lined up along the z-axis, the two illustrative cases in Figure 4. In each magnetic island, the closed field lines nest around an O neutral point. The island boundary is a part of a separatrix line, drawn as a dashed line in Figure 4, associated with an X neutral point. This separatrix line comes in from y → −∞, passes through the X neutral point to go around the magnetic island, passes through that neutral point again, and then goes off to y → ∞. As a increases past a₀, the special case of a = a₀ shown on the right in Figure 3 corresponds to the parametric simultaneous appearances of the O and X neutral points in coincidence, manifesting as a cusp on the dashed line shown.

The density ρ(y, z) concentrates around the vertical plane y = 0 in a vertical series of alternating maxima and saddle points. Along y = 0, the saddle points are minima in density alternating with the maxima, having the respective values

$$\begin{equation} \rho _{{\rm min}} = {B_2 B_3 k \over 2 \pi }(\sqrt{1 + a^2} + a )^{-2}, \end{equation} \tag{ 58 }$$

$$\begin{equation} \rho _{{\rm max}} = {B_2 B_3 k \over 2 \pi }(\sqrt{1 + a^2} + a )^2. \end{equation} \tag{ 59 }$$

Each of these equilibrium fluids when perturbed to evolve under the frozen-in condition moves with its flux surfaces of constant-A deforming as fluid surfaces, described by induction equation (19). The mass sandwiched between any two flux surfaces on which A = A₁, A₂, two constants,

$$\begin{equation} {\mathcal M}(A_1, A_2)=\int _{A_1<A(y, z, t)<A_2} \rho (y, z, t) dy dz, \end{equation} \tag{ 60 }$$

is a constant in time, where ρ(y, z, t) is the evolving density and A(y, z, t) satisfies Equation (19). To clarify our notation, with no loss of generality, we can fix A₁ to define a reference flux surface. Then ${\mathcal M}(A_1, A)$ is a function of a single variable A such that setting this variable A = A₂, a constant, picks out the second flux surface to define the invariant total mass measured from the reference surface. We call ${\mathcal M}(A_1, A)$ the (invariant) mass function. Its derivative $d{\mathcal M}(A_1, A)/dA$ is the total mass per unit A-flux, which we refer to as mass per unit flux for brevity. We can calculate ${\mathcal M}(A_1, A)$ from the field and density in the a-family of equilibrium states. By this calculation we can attribute the two-dimensionality of the equilibrium state to the fact that the mass per unit flux varies with A when a ≠ 0, that is, flux tubes of a fixed unit axial flux are loaded with unequal total masses. The field is thus unevenly bowed by the weight varying from flux tube to flux tube. The special case of a = 0 corresponds to a uniform loading of mass for each flux tube of a fixed axial flux. These flux tubes are identical, thus compatible with equilibrium varying only in the y dimension.

The total weight averaged over a full periodicity in the z-direction given by

$$\begin{eqnarray} W &=& {k \over 2 \pi }\int _0^{2\pi /k} \int _{-\infty }^{\infty } \rho g dy dz \nonumber \\ &=& {B_2 B_3 \over \pi } \end{eqnarray} \tag{ 61 }$$

is independent of parameter a, exactly the weight W given by Equation (54) for the a = 0 1D isothermal KS slab. We recover the formula k = πk₀W/B²₂ showing that the product of the total average weight W with the inverse isothermal scale height determines the width of the 2D slab, as well as the wavenumber of the vertical periodicity. This property characterizes the a-family of solutions in terms of their different mass functions (${\mathcal M}$) all having the same total mass 2πW/kg per periodic cell, e.g., −∞ < y < ∞, 0 < z < 2π/k. The four states of different values of a in Figures 3 and 4 thus have the significance of having the same total mass 2πW/kg being distributed in different manners over the A-flux.

4.1.3. The 0 < a < a₀ States

4.1.4. The a > a₀ States

4.1.5. Endless Intermittently Resistive Evolution

There is a family of analytic solutions to the isothermal Grad–Shrafanov equation (35) especially suited for illustrating 2D equilibrium states in which the mass per unit flux, $d{\mathcal M}(A_1, A)/dA$, is a discontinuous function (Low & Manchester 2000; Manchester & Low 2000). These solutions satisfy the Grad–Shrafanov equation (35) with p₀(A) being a discontinuous function of A.

Prescribe the constant-A field lines to be contours of constant ϕ(y, z), where

$$\begin{equation} \phi = \exp (-k_0z) + \cos (k_0 y). \end{equation} \tag{ 62 }$$

This simply means that A(y, z) = A(ϕ) and Equation (35) reduces to the following ordinary differential equations:

$$\begin{equation} (1-\phi ^2){d \over d\phi }\left({dA \over d\phi }\right)^2- 2 \phi \left({dA \over d\phi }\right)^2+ {d \over d\phi }Q^2 = 0, \end{equation} \tag{ 63 }$$

$$\begin{equation} \phi {d \over d\phi }\left({dA \over d\phi }\right)^2+2\left({dA \over d\phi }\right)^2+ 4\pi {d p_0\over d\phi } = 0, \end{equation} \tag{ 64 }$$

relating A(ϕ), Q(ϕ), and p₀(ϕ). Figure 5 displays the contours of constant ϕ(y, z), comprising horizontally periodic cells of U-shaped field lines resting on a bed of horizontally oriented field lines. The field line ϕ(y, z) = 1 separates the latter set of lines from the periodic cells, running up to z → ∞ and back as the boundary of each successive cell. To build an equilibrium of this family of solutions, we may prescribe p₀(ϕ) arbitrarily, that is, specifying the isothermal pressure distribution p(y, z) = p₀[ϕ(y, z)]exp (− k₀z). This specification is equivalent to prescribing $d{\mathcal M}(A_1, A)/dA$ on the field lines of a fixed shape. Equations (63) and (64) then determine A(ϕ) and Q(ϕ) in terms of the prescribed p₀(ϕ) and yield the equilibrium magnetic field given by Equation (8). Note that in Figure 5 the upper periodic U-shape fields can be loaded with vertically elongated fluids with their lower ends terminating at a height beneath which the field is of a contrastingly different topology, the latter suggestive of the observed quasi-steady macroscopic voids at the base of prominences (Berger et al. 2010, 2011).

Zoom In Zoom Out Reset image size

In the original studies (Low & Manchester 2000; Manchester & Low 2000), interest was concentrated on A(ϕ) having an everywhere continuous derivative, so that both the field and plasma are everywhere continuous in space. This family of solutions includes the case of p₀(ϕ) prescribed as a discontinuous function corresponding to a discontinuous $d{\mathcal M}(A_1, A)/dA$. Suppose that p₀(ϕ) is discontinuous across a flux surface ϕ = ϕ₁, a constant; then Equations (63) and (64) demand that both (dA/dϕ)² and Q² are also discontinuous across ϕ = ϕ₁. The latter implies that the magnetic pressure given by

$$\begin{equation} B^2=\left({dA \over d\phi }\right)^2 | \nabla \phi |^2 + Q^2 \end{equation} \tag{ 65 }$$

is discontinuous across ϕ = ϕ₁. Independent of whether or not p₀(ϕ) is continuous, the total pressure is

$$\begin{equation} p + {B^2 \over 8 \pi } = \left[P_0 - {1 \over 4 \pi } \int \left({dA \over d\phi }\right)^2 d\phi \right] \exp (-k_0z) + {\lambda _0 \over 8 \pi }, \end{equation} \tag{ 66 }$$

derived from the force-balance equation (4). Thus, we are assured that the discontinuities in fluid pressure and magnetic pressures in the above construction automatically satisfy the continuity of the total pressure as a requirement for the surface of discontinuity dA/dϕ to be in force balance. We leave the reader to explore the properties of this family of solutions.

The sets of continuous and discontinuous equilibrium solutions are infinite, spanned by all forms of $d{\mathcal M}(A_1, A)/dA$, all sharing the same field lines displayed in Figure 5. Clearly, the set of continuous equilibria is a subset of measure zero of the set of all realizable equilibria.

The nonlinear Grad–Shrafanov equation (21) for a prescribed temperature θ(A), constant on each flux surface but allowed to be of different values on different flux surfaces, is generally analytically intractable and non-trivial to solve numerically. We turn to the 2D extension of the 1D isothermal KS slab, in an entirely different geometry, by Low & Petrie (2005) for an explicit illustration of inevitable discrete currents arising from distributions of temperature and total frozen-in mass that are discontinuous across magnetic flux surfaces.

The 1D KS slab solution described by Equations (51)–(53) has three free parameters, (B₁, B₂, B₃). Let us set B₁ = 0, limiting our attention to fields that lie in the y − z plane, and redefine B₂ = B₀cos Φ, 2B₃ = B₀sin Φ in terms of two equivalent free parameters B₀ and Φ. With k = k₀B₃/B₂ = (1/2Λ₀)tan Φ, we obtain the solution

$$\begin{equation} A=B_0 \cos \Phi \left(z - 2 \Lambda _0 \log \left[ \cosh \left(\frac{1}{2 \Lambda _0}\tan \Phi (y-y_0) \right) \right]\right), \end{equation} \tag{ 67 }$$

$$\begin{equation} {\bf B}=B_0 \left[ 0, \cos \Phi , \sin \Phi \tanh \left(\frac{1}{2 \Lambda _0} \tan \Phi (y-y_0) \right) \right], \end{equation} \tag{ 68 }$$

$$\begin{equation} \rho g = {B_0^2 \over 8 \pi \Lambda _0} \sin ^2 \Phi {\rm sech}^2 \left(\frac{1}{2 \Lambda _0}\tan \Phi (y-y_0) \right), \end{equation} \tag{ 69 }$$

$$\begin{equation} p=\Lambda _0 \rho g, \end{equation} \tag{ 70 }$$

where we have introduced a constant translation displacement y₀ in the y-direction with no loss of generality. The weight defined by Equation (54) is now expressed as

$$\begin{equation} W=\frac{1}{4 \pi }B_0^2 \sin 2 \Phi. \end{equation} \tag{ 71 }$$

We have done no more than rewrite the 1D solution in another equivalent form.

What is interesting is that if we now formally take B₀ a constant in space but let Φ, Λ₀, y₀ be arbitrary functions of x, we still have a solution to the force-balance equation (4), which, for our purpose here, we rewrite as

$$\begin{equation} {1 \over 4 \pi } \left({\bf B} \cdot \nabla \right) {\bf B} - \nabla \left(p + {B^2 \over 8 \pi } \right) -\rho g {\hat{z}} = 0. \end{equation} \tag{ 72 }$$

First, the variation with x introduced retains the solenoidal condition on B. Direct evaluation shows that p + B²/8π = B²₀/8π, uniform in 3D space, and (B · ∇)B gives only a vertical component that balances the weight ρg exactly. Thus, the force balance is satisfied for all prescribed Φ(x), Λ₀(x), and y₀(x). This was shown by Low & Petrie (2005), except that they had not realized that the scale height Λ₀, too, may vary freely with x, allowing for a corresponding temperature T = T₀(x), for a solution of this kind. Finally, an arbitrary function y₀(x) displaces the solutions on different x-planes arbitrarily, creating quite complex structures (Low & Petrie 2005; Petrie & Low 2005).

Each of these 2D equilibrium states is a bowed magnetic field lying in planes of constant x, wherein B_y = B₀cos Φ(x) and the temperature T = T₀(x) are uniform, supporting the weight W(x), given by Equation (71), frozen into the field. The 1D KS slab corresponds to Φ, Λ₀, and y₀ being the same constants on all x-planes as magnetic flux surfaces. Arbitrary loadings of weight W(x) and different temperatures T₀(x) on these flux surfaces are admissible, giving rise to field patterns varying across the flux surfaces in terms of the sag of the bowed field, density, and field intensity as functions of x. The original study of Low & Petrie (2005) was interested in continuous variations.

The study in the present paper introduces a new interest, namely, the equilibrium states of fields endowed by their resistive past with a distribution of total mass in magnetic flux tubes that is not continuous across flux tube boundaries. Once the discontinuous distribution has been acquired and the frozen-in condition restored as an excellent approximation, the field will relax into an equilibrium state that is expected to be discontinuous in fluid and field properties across the flux surfaces. Taking the final equilibrium state to be isothermal on flux surfaces as a simplifying assumption, we have two independent origins for the discontinuity, one being the different constant temperatures on isothermal flux surfaces and the other being the different total mass per unit magnetic flux in flux tubes on either side of a flux surface. We present examples of these two properties.

Figure 1 shows three solutions of different total weight W at the same temperature T₀ and scale height Λ₀. A continuous variation with x through these three solutions gives an everywhere continuous equilibrium state varying with x and y. A prescription of W(x) that jumps discontinuously from two different states in Figure 1 across a flux surface x = x₀, a constant, would have the superimposed field lines shown on the left in Figure 6. Everywhere on x = x₀ we see a discontinuous B across this flux surface, a magnetic tangential discontinuity carrying a discrete current flowing in that flux surface. Note that independent of the choice of W(x) and T₀(x), the total pressure p + B²/8π = B²₀/8π is uniform in space and therefore continuous across x = x₀, satisfying the condition for x = x₀ to be in force balance. In this case, not only is the direction of B discontinuous across x = x₀, but so is B², from which follows that the density and pressure are both discontinuous across that surface. As long as electrical conductivity is perfect, this discontinuous equilibrium is admissible as a steady state, but our interest in obtaining this solution is to make the point that electrical conductivity in a real solar plasma is extremely high but finite, so that the formation of this singularity leads to resistive dissipation and magnetic reconnection.

Zoom In Zoom Out Reset image size

As pointed out by Low & Petrie (2005), the far fields on different x-planes have different inclinations to the horizontal because the sets of field lines in the x-planes are loaded with different weight W. It is simple to show that B in the far region y → ±∞, ρ → 0 and the sheared field is in a force-free state ∇ × B = α(x)B with α = ±dΦ(x)/dx. This field feeds into or takes out field-aligned currents from both sides of the 2D slab in order to regulate the cross-field current flowing in the x-direction within the slab. This current density distribution self-consistently accounts for just the right amount of cross-field current in the slab so that Lorentz force is nonzero and exactly supports the weight of the slab in the different x-planes.

The solution on the right in Figure 6 corresponds to the three solutions in Figure 2 with different temperature T₀ on different x-planes but the same weight W. Using these solutions to construct a 2D solution varying with x and y, let us suppose that T = T₀(x) jumps discontinuously from one state to another in Figure 2 across x = x₀. Then, the fields on the two sides of x = x₀ seen projected on x = x₀ are different, as shown on the right in Figure 6, producing a discrete current flowing on that surface. Again, the field, density, and pressure are all discontinuous across x = x₀. Since the weight W is the same for all x, that is, Φ is a constant, it follows that α = 0, and the far field in y → ±∞ does not vary with x, evident in the right panel of Figure 6. The current density generating the x-varying supporting Lorentz force in the 2D slab is not strictly in the x-direction, meandering around in the x−y plane but confined around the central plane y = 0 of the slab. Again, the discrete current in x = x₀ must dissipate at the small resistivity of a real plasma, leading to reconnection and changes in mass distribution over the magnetic flux.