- Email: [email protected]
Numerical simulation of the stationary state of periodically driven coronal loops
Computer Physics Communications 59 (1990) 75—84 North-Holland
75
NUMERICAL SIMULATION OF THE STATIONARY STATE OF PERIODICALLY DRiVEN CORONAL LOOPS Stefaan POEDTS
~,
Marcel GOOSSENS
Astronomisch Instituut, Katholieke Universiteit Leuven, Celestynenlaan 200 B, 3030 Heverlee, Belgium
and Wolfgang KERNER2 Max-Planck-Institutfür Plasmaphysik, Euratom Association, Boltzmannstral3e 2~D-8046 Garching be! München, Fed. Rep. Germany
The heating of solar coronal loops by resonant absorption of Alfvén waves is studied in the framework of linearized, compressible, resistive MHD by means of numerical simulations in which the ioops are approximated by straight cylindrical, axisymmetric plasma columns with equilibrium quantities varying only in the radial direction. The incident waves that excite the loops are modelled by a periodic external source. The stationary state of this driven system is determined numerically with a finite element code. The finite element technique is extremely suitable to compute the nearly-singular solutions and yields very accurate results. The efficiency of the heating mechanism and the energy deposition profile in this stationary state strongly depend on the characteristics of both the external driver and the equilibrium. A numerical survey of the relevant parameter space shows that resonant absorption is very efficient for typical coronal parameter values and appears to be a viable candidate heating mechanism for solar loops.
1. Introduction The search for the non-radiative heating mechanism(s) active in the solar corona is going on for several decades by now and its importance and significance has still grown in the last decade with the observations of hot, solar-like, coronae in many classes of stars [1—3].The aim of the investigations presented in this paper is to clarify the role resonant absorption of MHD waves can play in the heating of the dense magnetic loops that are observed in the solar corona. The ideal MHD spectrum of oscillation frequencies of an inhomogeneous plasma contains two continuous parts: an
2
Present address: Max-Planck-Institut fur Plasmaphysik, Euratom Association, Boltzmannstral3e 2, D-8046 Garching bei Munchen Fed. Rep. Germany. Present address: JET Joint Undertaking, Theory i~, Abingdon, Oxfordshire 0X14 3EA, England.
0010-4655/90/$03.50 © 1990
-
Alfvén continuum and a slow magnetosonic continuum, with corresponding modes that possess non-square-integrable spatial singularities. Therefore, when such a plasma is driven periodically by an external source with a frequency within the range of the continuous spectrum, the energy would accumulate unbounded (in ideal MHD) in an ever diminishing plasma layer around the singular point, as time progresses. However, due to dissipative effects, this unbounded growth will be prevented and after a finite time the driven system will attain a stationary state in which all physical quantities oscillate in time with a constant amplitude and with the frequency of the external driver. In this stationary state the energy dissipation rate in the plasma exactly balances the power input from the external source. Resonant absorption of Alfvén waves was first studied in the framework of controlled thermonuclear fusion research as a supplementary heat-
Elsevier Science Publishers B.V. (North-Holland)
76
S. Poedts et a!.
/ Numerical simulation
ing mechanism to bring laboratory plasmas into the ignition regime. After lonson [4] proposed resonant absorption of MHD surface waves as a means of heating solar coronal loops, this heating mechanism has been studied by many authors in this context (see the introduction of Hollweg and Yang [5] for a summary review). In the present paper, the heating of solar coronal ioops by resonant absorption of MHD waves is studied in the framework of linearized, resistive MHD. The stationary state of a cylindrical coronal loop model, excited by a periodic external source, is determined numerically and the efficiency of the heating mechanism is expressed in terms of the fraction of the power supplied by the external source that is actually dissipated and converted into heat. We concentrate on Alfvén wave heating because the resonant absorption of slow magnetosonic waves is much less efficient: for the relevant parameter range the fractional absorption never exceeds 1% in the slow continuum. The efficiency of resonant absorption of Alfvén waves and the corresponding energy deposition profile strongly depend on the frequency and the mode numbers of the external driving source and on the characteristics of the plasma. For typical coronal parameter values the resonant absorption process is very efficient in the stationary state and the heat deposition is localized in the outer plasma layers of the ioops. The physical model that is used to study the phenomenon of resonant absorption of Alfvén waves is presented in the next section. The numerical techniques that are applied to solve the equations of linearized resistive MHD and the accuracy of the obtained results are discussed in section 3. Section 4 contains an overview of the results of our numerical simulations. A discussion of these results and the conclusions of our investigations are presented in section 5.
2. Physical model In this section we first describe and motivate the cylindrical configuration and the class of equilibria we consider to model a solar coronal loop that is excited by waves that pitch in on it. Next,
of the stationary state of coronal loops
the equations of linear, resistive MHD are presented and briefly discussed. 2.1. Configuration Poedts et a!. [6] developed a computer code to simulate resonant absorption of Alfvén waves in laboratorium devices (tokamaks) in the context of controlled thermonuclear fusion research. They considered a configuration that consists of a cylindrical, axisymmetric plasma with radius ,~,, surrounded by a vacuum and by a perfectly conducting wall at radius r~.The plasma column is assumed to be excited by a radiating external “antenna” (an external coil) situated in the vacuum region, with radius r,~(r~ r~ r~).The boundary conditions applied at the plasma—vacuum interface relate the fields at the interface to those of the external antenna such that the current in the antenna coil scales the whole solution. Clearly, solar coronal loops are not surrounded by a vacuum region and a perfectly conducting wall and they are not excited by a current in an external coil either. Therefore, for the coronal loop case, the problem has been reformulated to simulate the resonant absorption of a wave with particular wavenumbers and frequency and with unit amplitude that is incident on the ioop. The role of the vacuum region has been reduced by placing the driver at the plasma surface (ra i~,,)and the wall has been taken away from the plasma column (r~ oo). In addition, the solution is renormalized such that the wave velocity field perturbation is of unit amplitude at the plasma surface, i.e. v( r.,,) II = 1. The cylindrical geometry closely approximates that of the loops observed in the solar corona which have a large aspect ratio, i.e. e L/~u~,,>> 1, with L the length of the loop. For the calculations presented in this paper we took ~ = 20, a typical value for coronal loops. The renormalization of the solution reduces the calculation to what Grossmann and Smith [7] call the “intrinsic absorption”, which is a measure for the efficiency of the heating mechanism in the particular loop considered. The determination of the actual heating of a solar coronal loop by resonant absorption requires the knowledge of the power spectrum of the external source and con—
—~
S. Poedts et al.
/ Numerical simulation of the stationary state of coronal loops
sists of taking the convolution of the “intrinsic” energy dissipation spectrum of the coronal loop with the power spectrum of the external source [7]. Since the power spectrum of the waves that are incident on the coronal loops is not known, only the efficiency of the heating coronal loopsheat by resonant absorption and the of corresponding deposition profiles can be determined. The procedure described above has the advantage to yield the source impedance which is a very useful quantity in a system that is excited by an external source, since it provides information about the fractions of the supplied energy that are dissipated and circulating in the system. Following Grossmann and Smith [7], the efficiency of the resonant absorption process is then expressed in terms of the fractional absorption, i.e. the fraction of the power supplied to the system by the external source that is actually dissipated and converted into heat. 2.2. Equilibrium model The heating of solar coronal loops by the resonant absorption of Alfvén waves is studied within the framework of linear, resistive MHD. The resistive MHD equations are linearized around an ideal static equilibrium. We consider a cylindrical plasma with direction equilibrium varying onlycolumn in the radial (thequantities main direction of inhomogeneity in solar loops). In the usual cylindrical coordinates (r, 0, z) the equilibrium force balance then reads dr
z
dr
d r ~ (rB 9)
=
o,
(i)
with P(r) and B(r) respectively the equilibrium plasma pressure and magnetic field (magnetic induction). Two of the three profiles P(r), B9(r), and B2(r) may be chosen arbitrarily; the remaining profile is determined by eq. (1). At the plasma—vacuum interface the total pressure (plasma pressure + magnetic pressure) and B~have to be continuous (the latter condition is automatically satisfied since B,. 0). The equilibrium mass density p(r) does not appear in the equilibrium equation and can be chosen arbitrarily. We con-
77
sider a class of equilibria given by the profiles B2 = 1, (2a) J =j~(i r2)~’, (2b) 2, (2c) pin=dimensionless 1 (1 d ) r units. The distance is normalized to the radius of the loop (re), and B 2 and p are normalized to their respective values at the axis, —
—
—
i.e. B2(0) and p(O). B9(r) can be computed from the equilibrium current density given by eq. (2b). The plasma pressure P then follows from the substitution of B2 and B9 in the equilibrium force balance (1). We consider equilibria with vanishing plasma pressure at the surface (P(1) = 0), such that there are no surface currents in the equilibrium model. This is the most realistic choice for the equilibrium. The constants d, J0 and v are free parameters. d denotes the density at the plasma boundary, d = p(l), and J0 is the current density at the axis (J0 = J(0)), with J0 = 2k/q0 (where k is the inverse aspect ratio of the loop and q0 the safety factor on the axis, q(r) = rkB2/ B9). The parameter iJ determines the profile of the electric current density and hence the shear of the equilibrium magnetic field since q(1)/q(O) = v + 1, for the considered class of equilibria. Notice that these equilibria2,are force-free. Theisplasma is anot function of r and of the beta, 2 P/B order fiof=1% for the parameter range considered in this paper. 2.3. Resistive MHD equations
The displacements resistive MHDabout equations that static governequithe linear this ideal librium are 0v p-~ = =
ab -~-
=
—
v’p
V
+
V/P
(t” x B)
+
(v x b) x B,
yP’17v,
x (v x B)
t,~’
xb
—
ç” x (~Vx b),
(3) (4) (5)
with v, p and b, the Eulerian variation of, respectively, the velocity, the plasma pressure and the magnetic field. ~ is the electric resistivity and y
78
5. Poedts et al.
/ Numericalsimulation of the stationary state of coronal loops
the adiabatic index (we assume y = ~). In this paper q is assumed to be constant over the plasma column. Equations (3)—(5) are expressed in dimensionless units where the time is normalized to the Alfvén-transit-time ta = with J’~= B2(0)/
‘i,/’c~
i~ö~). Equation (3)
is the momentum equation for a non-viscous plasma. Equation (4) is the ideal MHD equation for the variation of the internal energy (adiabatic law), while eq. (5), the induction equation, includes the Ohmic term (due to the finite electric conductivity of the plasma). Notice that there is no restriction to incompressible plasmas. Equations (3)—(5) form a system of 7 linear partial differential equations which is to be completed with appropriate boundary conditions.
the aspect ratio of the loop). n is the longitudinal mode number and m is the mode number in the 9-direction. When the dissipative system is excited by an external periodic source, it reaches a stationary state i© which all physical quantities vary periodically in time with the frequency of the external driver. This stationary state and the energy dissipation rate in it can be computed by substituting the temporal dependence
f ( r; t)
=
f(r) e~pt
(7)
in eqs. (3)—(5), where to1, is the frequency of the external source. With the substitutions (6) and (7) the system (3)—(5) reduces to a set of ordinary differential equations of order 6 in the radial coordinate r, which can be written in the form
3. Numerical procedure
Lu =0,
This section gives a description of the applied numerical techniques and discusses the accuracy of both the numerical and the physical approximations of the obtained results.
where L denotes the linear matrix operator and the state vector u is defined as u~= (v1 = rvr,
3.1. Numerical methods In the external region (approximated by a vacuum) the perturbed magnetic field, b~,satisfies the equations ~ X b,~=0, and ~7 b~=0, which can be solved analytically yielding b~in terms of modified Bessel functions [6]. This solution then has to be connected to the solution in the plasma column which is determined numerically. Since the equilibrium quantities do not depend on 0 and z in the considered equilibrium model for coronal loops, the perturbed quantities can be Fourier analyzed in these coordinates and each Fourierterm can be studied separately. The following separation ansatz is made for the perturbed quantities 1: f(r, 0, z; t) =f(r; t) ei(mof~hIc2). (6) The number k = TT/L defines a quantization factor to allow an integral number of half-wavelengths on the column (the ends of the loop are tied in the dense photosphere) with L the length of the loop in the z-direction (L = ~r, with the
(8)
v2 = it’8~ v3 = irv2, ~ = rp, b1 = irbr, b3 = rb2). The 0-component of the perturbed magnetic field has been eliminated with the divergence equation V b = 0 (provided m * 0). The system of differential equations (8) is solved in its weak form. The vector u(r) is a weak solution if for any function v(r) of the admissible Sobolev space satisfying the appropriate boundary conditions, the scalar product
S. Poedts et aL
/ Numerical simulation of the stationary state of coronal loops
spatial localization of the solutions at the near-singularity and in the boundary layer. The components of u are approximated by a finite linear combination of local expansion functions h~(r),
79
cients. The matrix B is symmetric and positivedefinite but A is non-Hermitian. A and B possess a tridiagonal block structure owing to the use of finite elements as expansion functions.
N
~ a~’h~(r), k=1,...,6,
(9)
3.2. Boundary conditions and drivingterm
1=1
where the 6N (with N the number of gridpoints) unknown coefficients a~ are to be determined. The development of the spectral codes for ideal MHD has indicated that every component of the state vector, u(r), has to be represented by appropriate finite elements. If the same elements are
The solutions of the 6th-order system of ODE’s (8) have to satisfy appropriate boundary conditions. For physical reasons the solutions have to be regular at the axis (r = 0) and it can be shown that the regularity conditions on v,. and br form a complete set. Hence, in terms of v1 and b1, we
used for all components, poor discretization is achieved because of the condition that the transverse divergence of the velocity vanish exactly in every interval [9], i.e.
have v~(0)= 0 and b1(0) = 0. The fact that the linear oscillations of the plasma column are excited by an external source is taken into account by imposing appropriate boundary conditions at the plasma—vacuum interface and by connecting
l’dv1 r dr
(
cannot be satisfied. If finite elements of order n + 1 are chosen for v1 and b1 and elements of order n for the remaining components, then the divergence condition (10) can be satisfied and the spectrum is well approximated numerically without pollution [101. Linear elements for v1 and b1 and piecewise-constant elements for v2, v3, ~ and b3 (as used in ideal spectral codes) yield a good numerical representation for zero resistivity. For finite q, however, the component b3 has to be differentiated, which also calls for the use of finite elements that are (at least) linear. We, therefore, introduced higher-order elements. Cubic Hermite spline functions are used for b1 and v1 and quadratic finite elements for v2, v3, ~ and b3. In each case two orthogonal functions define a cornplete set [8,111. Owing to the choice of two orthogonal finite elements per interval and per component of ii, the number of unknowns a~is raised to 12N. In the Galerkin method, which is applied here, the basis functions h~(r) are used in the weak form. This leads to a set of linear equations of the following form:
in this way the solutions in the plasma column to the solutions in the external region. At the (perturbed) plasma—vacuum interface the perturbed total pressure and the three components of the perturbed magnetic field have to be continuous. In the finite element method the regularity conditions on the axis are called essential boundary conditions because they have to be imposed on the basis functions themselves. In other words, the Sobolev-space of the basis functions has to be limited to those functions that satisfy the regularity conditions. The boundary conditions at the plasma—vacuum interface are called natural boundary conditions since they are satisfied automatically after a slight modification of eq. (11), the weak form of the system (8) [8]. This modification consists of substituting appropriate expressions (derived from the boundary conditions) in the surface terms that appear in the weak form (11) by partial integrations. This operation yields contributions to the matrix A and to a “driving” term, where the latter contributions are proportional to the driving current in the antenna coil [12]. Hence, in an externally driven system the Galerkin procedure leads to a system of linear equations of the form
Aa=ito1,Ba,
(Am—iwpB)a=f,
—
—
+
my2)
=
0,
(10)
(11)
with a the vector of the 12N expansion coeffi-
(12)
with a the vector of the 12N expansion coeffi-
S. Poedts et a!. / Numerical simulation of the stationary state of coronal loops
80
cients and / the “driving term”. Owing to the local nature of the finite elements the modification of the coefficient matrix is very small and involves only the last subblock. The dimension of the coefficient matrix (Am ico~B)is 12N X 12N. However, the use of finite elements as expansion fuctions gives it a tridiagonal block structure so that it can be stored in “band storage mode”. The system (12) is then solved with the LINPACK routines CGBFA and CGBSL. —
The left-hand side term of eq. (13) is related to ~aflt’ the power emitted by the antenna [6],
rant =
—
J
Jant
e~dv.
V,,, is the volume of the vacuum region and J~is the electric current density in the antenna. In the stationary state all physical quantities oscillate with the frequency of the driver and the resistive energy balance reads
3.3. Accuracy = OD. (14) This means that the fraction of the energy supplied by the external source that couples to the plasma, exactly balances the rate at which energy is converted into heat by Ohmic dissipation. The resistive energy equation (13) and the resis-
Re(Pant)
When the system (8) is solved the energy dissipation rate can be calculated in order to determine the efficiency of the heating process. However, let us first discuss numerical accuracy of the code and the validity of the physical assumptions made. An excellent diagnostic for this which, in addition, provides us with the complete picture of energy conservation in the system, has been built-in in the code. The linearized resistive MHD equations yield the following energy equation: —
J =
V
(e *
Jpv* +
j
{
X b)
dV
av ~ dV v~V/P
—
(~K) .
<
b
—~--
j dV
+f~J2dV.
static equilibrium, adiabatic law). The numerical accuracy and the physical approximations are good when the eqs. (13) and (14) are satisfied within reasonable bounds. The numerical accuracy is determined by the distribution of the meshpoints and can be raised to any level by enlarging the number of gridpoints and/or mesh accumulation in the nearly-singular layer and at the plasma boundary. The relative error in the energy balance because of the physical assumptions made, is typically 0.001% for 2501 gridpomts (2500 intervals) with mesh accumulation. These assumptions yield, therefore, a very good approximation.
0b* ~ +b
tive energy balance (14) provide excellent tests for the numerical accuracy of the code and for the validity of the physical approximations made (ideal
~ (mOD) (13)
e and j are the Eulenan perturbations of the electric field and the electric current density, respectively. The * denotes the complex conjugate and V~,is the volume of the plasma. According to eq. (13) an inflow of electromagnetic energy produces a rise of the kinetic energy of the plasma (K), a change of the potential energy of the plasma (Wa) and heat by Ohmic dissipation (OD).
4. Results In this section it is shown that the efficiency of the resonant absorption process and the localization of the intrinsic heating are affected by the characteristics of the external driving source and the coronal ioop equilibrium and that resonant absorption of Alfvén waves is very efficient for typical parameter values in solar coronal loops. The efficiency of the resonant absorption process
S. Poedts et al.
/ Numerical simulation of the stationarystate of coronal loops
is expressed in terms of the fractional absorption, ía, which is defined as Ohmic dissipation rate
—
7732
6
Clearly, 0
240.0
—
0.2
=
0.4
0.6
.~
0.8
1.0
r
Wa 0.2o~
b
=
0. 18~ 0.16~ 0 i4
4.1. Energy deposition profile
in r = t~, indicating that the ohmic heating is restricted to a narrow plasma layer around the ideally singular layer [13]. The position of this resonant in the depends frequency layer (con) and the plasma wavenumbers (m on and the n) of the wave that is incident on the coronal loop and on the variation of the physical quantities in the loop, as toa(r) is given by toa(r) = ((m/r)B 9(r) + nkB2(r))/ 1/~j.The width of the resonant layer, depends on theet physical quantities in the ~coronal loop. Poedts al. [13] showed that
a
10 8
OD
When a singie wave with frequency w~ inpitches on a solar loop, the plasma response has a nearly-singular behaviour in r = r~where tot, = toa(1~j, with Wa, the local Alfvén frequency. By consequence, ~jJ2(r) has a pronounced maximum
81
° iO~
0.08 0.0
.6
.8
r Fig. 1. The ~j2(r)-profiIes ,~=iO~,
for m = 1, n
=
1, q
0 = 1, d ‘~a(’~ = 0.2, and v=0.5 and 4(a), and the corresponding profiles (b).
2(r) for 50 driving frequencies in the range of ~J the ideal continuum (equidistant) and normalizing this sum such that its integral over the plasma volume equals one, hence ~~J2(to r) 2(r) nJ
=
f~~qj2(to
r) dV
“p
— 1/3 ‘~res
The “intrinsinc” heat deposition profile is typi-
ul”3Itoa(1s)
where the prime denotes the derivation with respect to r. However, this idealized case never occurs in nature. Solar loops are excited by many different waves with different frequencies each of which yields an ~J2(r)-profile that is peaked around the corresponding resonant layer. In linear theory the energy deposition profile in a loop that is excited with a uniform driving spectrum can be obtained by adding up the ~j2(r)-profiles for each of the driving frequencies. Poedts et al. [13] studied the localization of the intrinsic dissipation for a given equilibrium and for given modenumbers m and n by adding up the contributions
cally localized in the outer regions of the plasma column. The numerical survey carried out by Poedts et a!. [13] reveals that in most cases the i~J2(r)-profileis only slightly affected by the equilibrium parameters, at least as long as these parameters are varied within realistic bounds. On some special occasions, however, the heat is deposited more homogeneously as illustrated in fig. la, which shows the ~J2(r)-profiles for m = 1, n = 1, Jo = 0.1, d = 0.2, ~j = iO~, and v = 0.5 and 4. For the smooth current density profile (~ = 0.5) the Ohmic dissipation is almost uniformly spread over the coronal loop while for the peaked current density (v = 4) the intrinsic heating is localized in
82
S. Poedts et at
/
Numerical simulation of thestationary state of coronal ioops
the outer plasma layers. This is because for v = 0.5 the profile of the local Alfvén frequency is monotonic and the chosen driving frequencies (equidistantly spread over [Wa,,~,, toam,,j) yield resonant surfaces everywhere in the coronal loop, while for v = 4 the toa-profile is non-monotomc and the resonant surfaces are shifted to the outer regions of the plasma column (see fig. ib). 4.2. Efficiency Poedts et al. [12] noticed that not every driving frequency yields equally good coupling of the external source to the plasma. Extremely high efficiency is obtained when the driving frequency is close to the real part of the frequency of a so-called “collective mode” or “damped quasi-mode” which is the remnant of a discrete mode (in this case the stabilized external kink mode) that disappeared in the continuous spectrum [14,15]. For a bandwidth around the optimal driving frequency ía max(la(wp)}, and in the following the range of driving frequencies for which the ía 0.95 max{ la(COp)} is called the “efficient” bandwidth; clearly an important parameter for the efficiency of resonant absorption as a coronal loop heating mechanism. The efficiency of resonant absorption also depends strongly on the wavenumbers of the waves that are incident on the coronal loops as iilustrated in fig. 2, where the fractional absorption is shown versus the driving frequency for n = 5 and m = — 5—5, and where the other parameters are fixed (see the figure caption). Of course, as is clear from the definition of coa(r) (see above), the
f i .~
-‘7
0.81 0.6
1
~
°2L~~~~’ / ~ -
00
0.20.4
0.6
5 o~
Fig. 2. The fractional absorptionp versus the frequency of the driver for v=2, d=0.2, j,,=0.1, ,~=10_8,n=5, and m= —5, —3, —1,1,3 and 5.
range of the ideal Alfvén continuum — and so the range of frequencies that can be resonantly excited — depends on m. By consequence, many different driving frequencies contribute to the heating of the same plasma layers as the corresponding resonant layers range from r = 0 to r = 1 for each value of m. The frequencies in fig. 2 are dimensionless and related to physical frequencies a by the transformation formula i r —1 a(Hz)=0.1417co( ~“ \ 10 m / B (0) \ i p (0) - 1/2 x ~ doi T io’~ _3) m where we assumed a mean atomic weight of 0.6 (for fully-ionised H plasma ~i= 0.5; in the solar atmosphere the presence of extra elements makes ~ 0.6 [16]). Hence, for a loop with a radius of 5 x 106 m, a magnetic field of 0.005 T (50 Gauss) on the axis and a number density on the axis of 3, for example, the ideal Alfvén con1016 m tinuum reaches from 0.00896 Hz to 0.04009 Hz for m = —3, from 0.02689 Hz to 0.05345 Hz for m = 1, and from 0.04482 Hz to 0.06681 Hz for m = 5 (for the parameters used in fig. 2). So the range of frequencies at which this coronal loop can be resonantly excited agrees very well with the observed motions in the solar atmosphere which have periods of typically 100 s. The m = 1 mode indeed yields the strongest coupling and so the most efficient heating. For this mode ía = 1 for the optimal driving frequency and ía 90% for a large bandwidth. The m = —1 mode couples almost as efficient to the plasma as the m = 1 mode and this feature can be generalized: max( ía(top)} and the range of the efficient bandwidth are almost independent of the sign of the poloidal wavenumber. The fractional absorption gradually decreases as m increases and for m = 5, max( í (to )) amounts to 28%. For m = 3 the fracp tional aabsorption is 48% for the optimal driving frequency. When the wavenumber m is fixed and n is varied, it is found that the maximal value of ía
) (~
(for the optimal driving frequency) does not depend . on the value of n when m is positive. For negative values of m, however, the optimal efficiency is higher for the higher values of n [13].
S. Poedts et al.
/ Numerical simulation of the stationary state of coronal loops
OD 1.000
~
0 995
-~
o.goo
83
on the axis). For typical solar loop parameter values the energy supplied by the exciting waves is very efficiently dissipated and converted to heat.
/
5. Discussion and conclusions
0.98~ ,.—.~.—=.-=,= 6 7 B
9
10
11
12
The heating of solar coronal loops by resonant absorption of Alfvén waves has been investigated Fig. 3. The Ohmic dissipation rate versus the resistivity (in . . arbitrary units, chosen such that OD 1 for Rm 1012) for in the framework of hneanzed compressible restscop 0.195, m 2, a =1, q 0 Q5, p =1 and d =1. tive MHD. The resonant absorption of the waves that are incident on the coronal loops is numerically simulated in straight cylindrical loop models, Figure 3 shows how strikingly insensitive the externally excited by a periodic source. The s/aOhmic dissipation rate in the stationary state is to tionary state of this driven system and the Ohmic the electric conductivity of the plasma: enlarging dissipation rate in this state are determined by the magnetic Reynolds number Rm (m 1) by six means of a very accurate code based on the finite orders of magnitude changes the Ohmic dissipaelement technique which is extremely suitable for tion rate (in dimensionless units) by less than this purpose because of the strong localization 1.5%. Moreover, the energy deposition rate (due to the nearly-singular and boundary layer reaches an asymptotic value for Rm> 108. Since behaviour) of the solutions. We followed Grossthe resistivity of coronal loop plasmas is very mann and Smith [7] and computed the “intrinsic” small (i~ 10 8), the dissipation rate is almost dissipation; i.e. the power spectrum of the external independent of i~ for the relevant values of this source that excites the magnetic loops in the solar parameter. However, the heating rate due to resocorona is not specified in the numerical simulanant absorption does depend on the other equitions (it is assumed to be uniform). The actual libnum quantities. We did not succeed to confirm heating of a coronal loop can only be determined Grossmann’s and Smith’s [7] expression for the when the driving spectrum of the loop is known. ideal decay rate (their eq. (47)) which should corHowever, our numerical simulations give a very respond to the heating rate in a dissipative system. good picture of the efficiency of the resonant Moreover, we couhi not find an alternative forabsorption process in solar loops and yield valmula to express the Ohmic dissipation rate in uable information on the viability of this process those equilibrium parameters. It seems there are as a coronal heating mechanism. one or more parameters missing. One such paramThe energy supplied by incident waves on eter could surely be the proximity of the driving coronal loops, with frequencies within the range of frequency to the frequency of the “collective the Alfvén continuum, is very efficiently dismode”. In an attempt to elucidate the effect of the sipated in typical solar coronal loops. However, equilibrium parameters on the efficiency of resothe role of resonant absorption in coronal heating nant absorption in solar ioops, Poedts et al. [13] depends dramatically on the magnitude of the computed la(c~p)-profilesfor different choices of background magnetic field. Indeed, the energy the equilibrium parameters. They found that the dissipation rate due to resonant absorption is pro3 (as can be found by transforefficiency of the mechanism, by portional max{ la(t0p)} andheating by the range of measured the efficient ming backtotoB2(0) quantities with dimensions (see seebandwidth of driving frequencies, is substantially tion 2)) and so the magnitude of the background influenced by the parameter d (the plasma density field plays a crucial role in the amount of energy of the loop surface) as wall as by the parameters v that is converted into heat per unit of time. This (related to the shear of the equilibrium magnetic result is very promising as the observed magnetic field (see section 2.2)) and J 0 (the current density fields are rather strong in the solar corona. Other Log(R4
.
=
=
=
=
=
84
S. Poedts et al.
/
Numericalsimulation of the stationary state of coronal loops
indications in favour of resonant absorption are the localization of the heat deposition and the fact that the efficiency of the heating mechanism is independent of the length of the coronal loop. The intrinsic dissipation is localized in the outer plasma layers indicating that the absorption process is most efficient in the outer regions of the coronal loops, consistent with the temperature profiles ohserved in some loops [17]. Parker [18] remarked that the surface brightness of the active regions is almost independent of the dimensions. The only way where the dimensions of the coronal loop enter the numerics in our work is via the quantization factor k in the Fourier decomposition in the longitudinal direction (eq. 6). Given the fact that longer loops have larger radii and so the aspect ratio (and so k) is approximately the same for all coronal loops, the process of resonant absorption is equally efficient in all coronal loops to matter what their dimensions are. We conclude that resonant absorption is a viable candidate heating mechanism for coronal loops.
Acknowledgements S. Poedts is Research Assistant of the Research Council of the K.U. Leuven whose support is greatly acknowledged. We thank J.P. Goedbloed, J.V. Hollweg, B. Roberts and S. Antiochos for stimulating discussions and suggestions. The K.U. Leuven—MPIPP collaboration is supported by the European Economic Community Scientific Cooperation Contract No. ST2J-0275-C.
References [1] R. Rosner, L. Gollub and G.S. Vaiana, Ann. Rev. Astron. Astrophys. 23 (1985) 413. 12] J.L. Linsky, Sol. Phys. 100 (1985) 333. [3] C. Jordan and J.L. Linsky, in: Exploring the Universe with IUE Sateffite, ed. Y. Kondo (Reidel, Dordrecht, 1987)the p. 259. [4] J.A. lonson, Astrophys. J. 226 (1978) 650. [5]
161
J.V. Hollweg and G. Yang, J. Geophys. Res. 93 (1988) 5423. S. Poedts, W. Kerner and M. Goossens, J. Plasma Phys.
42 (1989) 27. [7] W. Grossmann and R.A. Smith, Astrophys. J. 332 (1988) 476. [8] G. Strang and G.J. Fix, An Analysis of the Finite Element Method (Prentice-Hall, Englewood Cliffs, NJ, 1973). [9] K. Appert, D. Berger, R. Gruber and J. Rappaz, J. Cornput. Phys. 18 (1975) 284). [10] J. Rappaz, Num. Math. 28 (1975) 15. [11] W. Kerner, K. Lerbinger, K. Gruber and T. Tsunematsu, Comput. Phys. Cominun. 36 (1985) 225. [12] S. Poedts, M. Goossens and W. Kerner, Sol. Phys. 123 (1989) 83. [13] S. Poedts, M. Goossens and W. Kerner, Astrophys. J., in press. [14] D. Van Eester, M. Goossens and S. Poedts J. Plasma Phys, submitted. [151 J.P. Goedbloed, Physica D 12 (1984) 107. [16] E.R. Priest, Solar Magnetohydrodynaniics (Reidel, Dord1984). Sol. Phys. 43 (1975) 327. [171 recht, P.V. Foukal, [18] E.M. Parker, Heating of the stellar corona, in: Coronal and Prominence Plasmas, ed. AL. Poland, NASA-CP-2442 (1986).
75
NUMERICAL SIMULATION OF THE STATIONARY STATE OF PERIODICALLY DRiVEN CORONAL LOOPS Stefaan POEDTS
~,
Marcel GOOSSENS
Astronomisch Instituut, Katholieke Universiteit Leuven, Celestynenlaan 200 B, 3030 Heverlee, Belgium
and Wolfgang KERNER2 Max-Planck-Institutfür Plasmaphysik, Euratom Association, Boltzmannstral3e 2~D-8046 Garching be! München, Fed. Rep. Germany
The heating of solar coronal loops by resonant absorption of Alfvén waves is studied in the framework of linearized, compressible, resistive MHD by means of numerical simulations in which the ioops are approximated by straight cylindrical, axisymmetric plasma columns with equilibrium quantities varying only in the radial direction. The incident waves that excite the loops are modelled by a periodic external source. The stationary state of this driven system is determined numerically with a finite element code. The finite element technique is extremely suitable to compute the nearly-singular solutions and yields very accurate results. The efficiency of the heating mechanism and the energy deposition profile in this stationary state strongly depend on the characteristics of both the external driver and the equilibrium. A numerical survey of the relevant parameter space shows that resonant absorption is very efficient for typical coronal parameter values and appears to be a viable candidate heating mechanism for solar loops.
1. Introduction The search for the non-radiative heating mechanism(s) active in the solar corona is going on for several decades by now and its importance and significance has still grown in the last decade with the observations of hot, solar-like, coronae in many classes of stars [1—3].The aim of the investigations presented in this paper is to clarify the role resonant absorption of MHD waves can play in the heating of the dense magnetic loops that are observed in the solar corona. The ideal MHD spectrum of oscillation frequencies of an inhomogeneous plasma contains two continuous parts: an
2
Present address: Max-Planck-Institut fur Plasmaphysik, Euratom Association, Boltzmannstral3e 2, D-8046 Garching bei Munchen Fed. Rep. Germany. Present address: JET Joint Undertaking, Theory i~, Abingdon, Oxfordshire 0X14 3EA, England.
0010-4655/90/$03.50 © 1990
-
Alfvén continuum and a slow magnetosonic continuum, with corresponding modes that possess non-square-integrable spatial singularities. Therefore, when such a plasma is driven periodically by an external source with a frequency within the range of the continuous spectrum, the energy would accumulate unbounded (in ideal MHD) in an ever diminishing plasma layer around the singular point, as time progresses. However, due to dissipative effects, this unbounded growth will be prevented and after a finite time the driven system will attain a stationary state in which all physical quantities oscillate in time with a constant amplitude and with the frequency of the external driver. In this stationary state the energy dissipation rate in the plasma exactly balances the power input from the external source. Resonant absorption of Alfvén waves was first studied in the framework of controlled thermonuclear fusion research as a supplementary heat-
Elsevier Science Publishers B.V. (North-Holland)
76
S. Poedts et a!.
/ Numerical simulation
ing mechanism to bring laboratory plasmas into the ignition regime. After lonson [4] proposed resonant absorption of MHD surface waves as a means of heating solar coronal loops, this heating mechanism has been studied by many authors in this context (see the introduction of Hollweg and Yang [5] for a summary review). In the present paper, the heating of solar coronal ioops by resonant absorption of MHD waves is studied in the framework of linearized, resistive MHD. The stationary state of a cylindrical coronal loop model, excited by a periodic external source, is determined numerically and the efficiency of the heating mechanism is expressed in terms of the fraction of the power supplied by the external source that is actually dissipated and converted into heat. We concentrate on Alfvén wave heating because the resonant absorption of slow magnetosonic waves is much less efficient: for the relevant parameter range the fractional absorption never exceeds 1% in the slow continuum. The efficiency of resonant absorption of Alfvén waves and the corresponding energy deposition profile strongly depend on the frequency and the mode numbers of the external driving source and on the characteristics of the plasma. For typical coronal parameter values the resonant absorption process is very efficient in the stationary state and the heat deposition is localized in the outer plasma layers of the ioops. The physical model that is used to study the phenomenon of resonant absorption of Alfvén waves is presented in the next section. The numerical techniques that are applied to solve the equations of linearized resistive MHD and the accuracy of the obtained results are discussed in section 3. Section 4 contains an overview of the results of our numerical simulations. A discussion of these results and the conclusions of our investigations are presented in section 5.
2. Physical model In this section we first describe and motivate the cylindrical configuration and the class of equilibria we consider to model a solar coronal loop that is excited by waves that pitch in on it. Next,
of the stationary state of coronal loops
the equations of linear, resistive MHD are presented and briefly discussed. 2.1. Configuration Poedts et a!. [6] developed a computer code to simulate resonant absorption of Alfvén waves in laboratorium devices (tokamaks) in the context of controlled thermonuclear fusion research. They considered a configuration that consists of a cylindrical, axisymmetric plasma with radius ,~,, surrounded by a vacuum and by a perfectly conducting wall at radius r~.The plasma column is assumed to be excited by a radiating external “antenna” (an external coil) situated in the vacuum region, with radius r,~(r~ r~ r~).The boundary conditions applied at the plasma—vacuum interface relate the fields at the interface to those of the external antenna such that the current in the antenna coil scales the whole solution. Clearly, solar coronal loops are not surrounded by a vacuum region and a perfectly conducting wall and they are not excited by a current in an external coil either. Therefore, for the coronal loop case, the problem has been reformulated to simulate the resonant absorption of a wave with particular wavenumbers and frequency and with unit amplitude that is incident on the ioop. The role of the vacuum region has been reduced by placing the driver at the plasma surface (ra i~,,)and the wall has been taken away from the plasma column (r~ oo). In addition, the solution is renormalized such that the wave velocity field perturbation is of unit amplitude at the plasma surface, i.e. v( r.,,) II = 1. The cylindrical geometry closely approximates that of the loops observed in the solar corona which have a large aspect ratio, i.e. e L/~u~,,>> 1, with L the length of the loop. For the calculations presented in this paper we took ~ = 20, a typical value for coronal loops. The renormalization of the solution reduces the calculation to what Grossmann and Smith [7] call the “intrinsic absorption”, which is a measure for the efficiency of the heating mechanism in the particular loop considered. The determination of the actual heating of a solar coronal loop by resonant absorption requires the knowledge of the power spectrum of the external source and con—
—~
S. Poedts et al.
/ Numerical simulation of the stationary state of coronal loops
sists of taking the convolution of the “intrinsic” energy dissipation spectrum of the coronal loop with the power spectrum of the external source [7]. Since the power spectrum of the waves that are incident on the coronal loops is not known, only the efficiency of the heating coronal loopsheat by resonant absorption and the of corresponding deposition profiles can be determined. The procedure described above has the advantage to yield the source impedance which is a very useful quantity in a system that is excited by an external source, since it provides information about the fractions of the supplied energy that are dissipated and circulating in the system. Following Grossmann and Smith [7], the efficiency of the resonant absorption process is then expressed in terms of the fractional absorption, i.e. the fraction of the power supplied to the system by the external source that is actually dissipated and converted into heat. 2.2. Equilibrium model The heating of solar coronal loops by the resonant absorption of Alfvén waves is studied within the framework of linear, resistive MHD. The resistive MHD equations are linearized around an ideal static equilibrium. We consider a cylindrical plasma with direction equilibrium varying onlycolumn in the radial (thequantities main direction of inhomogeneity in solar loops). In the usual cylindrical coordinates (r, 0, z) the equilibrium force balance then reads dr
z
dr
d r ~ (rB 9)
=
o,
(i)
with P(r) and B(r) respectively the equilibrium plasma pressure and magnetic field (magnetic induction). Two of the three profiles P(r), B9(r), and B2(r) may be chosen arbitrarily; the remaining profile is determined by eq. (1). At the plasma—vacuum interface the total pressure (plasma pressure + magnetic pressure) and B~have to be continuous (the latter condition is automatically satisfied since B,. 0). The equilibrium mass density p(r) does not appear in the equilibrium equation and can be chosen arbitrarily. We con-
77
sider a class of equilibria given by the profiles B2 = 1, (2a) J =j~(i r2)~’, (2b) 2, (2c) pin=dimensionless 1 (1 d ) r units. The distance is normalized to the radius of the loop (re), and B 2 and p are normalized to their respective values at the axis, —
—
—
i.e. B2(0) and p(O). B9(r) can be computed from the equilibrium current density given by eq. (2b). The plasma pressure P then follows from the substitution of B2 and B9 in the equilibrium force balance (1). We consider equilibria with vanishing plasma pressure at the surface (P(1) = 0), such that there are no surface currents in the equilibrium model. This is the most realistic choice for the equilibrium. The constants d, J0 and v are free parameters. d denotes the density at the plasma boundary, d = p(l), and J0 is the current density at the axis (J0 = J(0)), with J0 = 2k/q0 (where k is the inverse aspect ratio of the loop and q0 the safety factor on the axis, q(r) = rkB2/ B9). The parameter iJ determines the profile of the electric current density and hence the shear of the equilibrium magnetic field since q(1)/q(O) = v + 1, for the considered class of equilibria. Notice that these equilibria2,are force-free. Theisplasma is anot function of r and of the beta, 2 P/B order fiof=1% for the parameter range considered in this paper. 2.3. Resistive MHD equations
The displacements resistive MHDabout equations that static governequithe linear this ideal librium are 0v p-~ = =
ab -~-
=
—
v’p
V
+
V/P
(t” x B)
+
(v x b) x B,
yP’17v,
x (v x B)
t,~’
xb
—
ç” x (~Vx b),
(3) (4) (5)
with v, p and b, the Eulerian variation of, respectively, the velocity, the plasma pressure and the magnetic field. ~ is the electric resistivity and y
78
5. Poedts et al.
/ Numericalsimulation of the stationary state of coronal loops
the adiabatic index (we assume y = ~). In this paper q is assumed to be constant over the plasma column. Equations (3)—(5) are expressed in dimensionless units where the time is normalized to the Alfvén-transit-time ta = with J’~= B2(0)/
‘i,/’c~
i~ö~). Equation (3)
is the momentum equation for a non-viscous plasma. Equation (4) is the ideal MHD equation for the variation of the internal energy (adiabatic law), while eq. (5), the induction equation, includes the Ohmic term (due to the finite electric conductivity of the plasma). Notice that there is no restriction to incompressible plasmas. Equations (3)—(5) form a system of 7 linear partial differential equations which is to be completed with appropriate boundary conditions.
the aspect ratio of the loop). n is the longitudinal mode number and m is the mode number in the 9-direction. When the dissipative system is excited by an external periodic source, it reaches a stationary state i© which all physical quantities vary periodically in time with the frequency of the external driver. This stationary state and the energy dissipation rate in it can be computed by substituting the temporal dependence
f ( r; t)
=
f(r) e~pt
(7)
in eqs. (3)—(5), where to1, is the frequency of the external source. With the substitutions (6) and (7) the system (3)—(5) reduces to a set of ordinary differential equations of order 6 in the radial coordinate r, which can be written in the form
3. Numerical procedure
Lu =0,
This section gives a description of the applied numerical techniques and discusses the accuracy of both the numerical and the physical approximations of the obtained results.
where L denotes the linear matrix operator and the state vector u is defined as u~= (v1 = rvr,
3.1. Numerical methods In the external region (approximated by a vacuum) the perturbed magnetic field, b~,satisfies the equations ~ X b,~=0, and ~7 b~=0, which can be solved analytically yielding b~in terms of modified Bessel functions [6]. This solution then has to be connected to the solution in the plasma column which is determined numerically. Since the equilibrium quantities do not depend on 0 and z in the considered equilibrium model for coronal loops, the perturbed quantities can be Fourier analyzed in these coordinates and each Fourierterm can be studied separately. The following separation ansatz is made for the perturbed quantities 1: f(r, 0, z; t) =f(r; t) ei(mof~hIc2). (6) The number k = TT/L defines a quantization factor to allow an integral number of half-wavelengths on the column (the ends of the loop are tied in the dense photosphere) with L the length of the loop in the z-direction (L = ~r, with the
(8)
v2 = it’8~ v3 = irv2, ~ = rp, b1 = irbr, b3 = rb2). The 0-component of the perturbed magnetic field has been eliminated with the divergence equation V b = 0 (provided m * 0). The system of differential equations (8) is solved in its weak form. The vector u(r) is a weak solution if for any function v(r) of the admissible Sobolev space satisfying the appropriate boundary conditions, the scalar product
S. Poedts et aL
/ Numerical simulation of the stationary state of coronal loops
spatial localization of the solutions at the near-singularity and in the boundary layer. The components of u are approximated by a finite linear combination of local expansion functions h~(r),
79
cients. The matrix B is symmetric and positivedefinite but A is non-Hermitian. A and B possess a tridiagonal block structure owing to the use of finite elements as expansion functions.
N
~ a~’h~(r), k=1,...,6,
(9)
3.2. Boundary conditions and drivingterm
1=1
where the 6N (with N the number of gridpoints) unknown coefficients a~ are to be determined. The development of the spectral codes for ideal MHD has indicated that every component of the state vector, u(r), has to be represented by appropriate finite elements. If the same elements are
The solutions of the 6th-order system of ODE’s (8) have to satisfy appropriate boundary conditions. For physical reasons the solutions have to be regular at the axis (r = 0) and it can be shown that the regularity conditions on v,. and br form a complete set. Hence, in terms of v1 and b1, we
used for all components, poor discretization is achieved because of the condition that the transverse divergence of the velocity vanish exactly in every interval [9], i.e.
have v~(0)= 0 and b1(0) = 0. The fact that the linear oscillations of the plasma column are excited by an external source is taken into account by imposing appropriate boundary conditions at the plasma—vacuum interface and by connecting
l’dv1 r dr
(
cannot be satisfied. If finite elements of order n + 1 are chosen for v1 and b1 and elements of order n for the remaining components, then the divergence condition (10) can be satisfied and the spectrum is well approximated numerically without pollution [101. Linear elements for v1 and b1 and piecewise-constant elements for v2, v3, ~ and b3 (as used in ideal spectral codes) yield a good numerical representation for zero resistivity. For finite q, however, the component b3 has to be differentiated, which also calls for the use of finite elements that are (at least) linear. We, therefore, introduced higher-order elements. Cubic Hermite spline functions are used for b1 and v1 and quadratic finite elements for v2, v3, ~ and b3. In each case two orthogonal functions define a cornplete set [8,111. Owing to the choice of two orthogonal finite elements per interval and per component of ii, the number of unknowns a~is raised to 12N. In the Galerkin method, which is applied here, the basis functions h~(r) are used in the weak form. This leads to a set of linear equations of the following form:
in this way the solutions in the plasma column to the solutions in the external region. At the (perturbed) plasma—vacuum interface the perturbed total pressure and the three components of the perturbed magnetic field have to be continuous. In the finite element method the regularity conditions on the axis are called essential boundary conditions because they have to be imposed on the basis functions themselves. In other words, the Sobolev-space of the basis functions has to be limited to those functions that satisfy the regularity conditions. The boundary conditions at the plasma—vacuum interface are called natural boundary conditions since they are satisfied automatically after a slight modification of eq. (11), the weak form of the system (8) [8]. This modification consists of substituting appropriate expressions (derived from the boundary conditions) in the surface terms that appear in the weak form (11) by partial integrations. This operation yields contributions to the matrix A and to a “driving” term, where the latter contributions are proportional to the driving current in the antenna coil [12]. Hence, in an externally driven system the Galerkin procedure leads to a system of linear equations of the form
Aa=ito1,Ba,
(Am—iwpB)a=f,
—
—
+
my2)
=
0,
(10)
(11)
with a the vector of the 12N expansion coeffi-
(12)
with a the vector of the 12N expansion coeffi-
S. Poedts et a!. / Numerical simulation of the stationary state of coronal loops
80
cients and / the “driving term”. Owing to the local nature of the finite elements the modification of the coefficient matrix is very small and involves only the last subblock. The dimension of the coefficient matrix (Am ico~B)is 12N X 12N. However, the use of finite elements as expansion fuctions gives it a tridiagonal block structure so that it can be stored in “band storage mode”. The system (12) is then solved with the LINPACK routines CGBFA and CGBSL. —
The left-hand side term of eq. (13) is related to ~aflt’ the power emitted by the antenna [6],
rant =
—
J
Jant
e~dv.
V,,, is the volume of the vacuum region and J~is the electric current density in the antenna. In the stationary state all physical quantities oscillate with the frequency of the driver and the resistive energy balance reads
3.3. Accuracy = OD. (14) This means that the fraction of the energy supplied by the external source that couples to the plasma, exactly balances the rate at which energy is converted into heat by Ohmic dissipation. The resistive energy equation (13) and the resis-
Re(Pant)
When the system (8) is solved the energy dissipation rate can be calculated in order to determine the efficiency of the heating process. However, let us first discuss numerical accuracy of the code and the validity of the physical assumptions made. An excellent diagnostic for this which, in addition, provides us with the complete picture of energy conservation in the system, has been built-in in the code. The linearized resistive MHD equations yield the following energy equation: —
J =
V
(e *
Jpv* +
j
{
X b)
dV
av ~ dV v~V/P
—
(~K) .
<
b
—~--
j dV
+f~J2dV.
static equilibrium, adiabatic law). The numerical accuracy and the physical approximations are good when the eqs. (13) and (14) are satisfied within reasonable bounds. The numerical accuracy is determined by the distribution of the meshpoints and can be raised to any level by enlarging the number of gridpoints and/or mesh accumulation in the nearly-singular layer and at the plasma boundary. The relative error in the energy balance because of the physical assumptions made, is typically 0.001% for 2501 gridpomts (2500 intervals) with mesh accumulation. These assumptions yield, therefore, a very good approximation.
0b* ~ +b
tive energy balance (14) provide excellent tests for the numerical accuracy of the code and for the validity of the physical approximations made (ideal
~ (mOD) (13)
e and j are the Eulenan perturbations of the electric field and the electric current density, respectively. The * denotes the complex conjugate and V~,is the volume of the plasma. According to eq. (13) an inflow of electromagnetic energy produces a rise of the kinetic energy of the plasma (K), a change of the potential energy of the plasma (Wa) and heat by Ohmic dissipation (OD).
4. Results In this section it is shown that the efficiency of the resonant absorption process and the localization of the intrinsic heating are affected by the characteristics of the external driving source and the coronal ioop equilibrium and that resonant absorption of Alfvén waves is very efficient for typical parameter values in solar coronal loops. The efficiency of the resonant absorption process
S. Poedts et al.
/ Numerical simulation of the stationarystate of coronal loops
is expressed in terms of the fractional absorption, ía, which is defined as Ohmic dissipation rate
—
7732
6
Clearly, 0
240.0
—
0.2
=
0.4
0.6
.~
0.8
1.0
r
Wa 0.2o~
b
=
0. 18~ 0.16~ 0 i4
4.1. Energy deposition profile
in r = t~, indicating that the ohmic heating is restricted to a narrow plasma layer around the ideally singular layer [13]. The position of this resonant in the depends frequency layer (con) and the plasma wavenumbers (m on and the n) of the wave that is incident on the coronal loop and on the variation of the physical quantities in the loop, as toa(r) is given by toa(r) = ((m/r)B 9(r) + nkB2(r))/ 1/~j.The width of the resonant layer, depends on theet physical quantities in the ~coronal loop. Poedts al. [13] showed that
a
10 8
OD
When a singie wave with frequency w~ inpitches on a solar loop, the plasma response has a nearly-singular behaviour in r = r~where tot, = toa(1~j, with Wa, the local Alfvén frequency. By consequence, ~jJ2(r) has a pronounced maximum
81
° iO~
0.08 0.0
.6
.8
r Fig. 1. The ~j2(r)-profiIes ,~=iO~,
for m = 1, n
=
1, q
0 = 1, d ‘~a(’~ = 0.2, and v=0.5 and 4(a), and the corresponding profiles (b).
2(r) for 50 driving frequencies in the range of ~J the ideal continuum (equidistant) and normalizing this sum such that its integral over the plasma volume equals one, hence ~~J2(to r) 2(r) nJ
=
f~~qj2(to
r) dV
“p
— 1/3 ‘~res
The “intrinsinc” heat deposition profile is typi-
ul”3Itoa(1s)
where the prime denotes the derivation with respect to r. However, this idealized case never occurs in nature. Solar loops are excited by many different waves with different frequencies each of which yields an ~J2(r)-profile that is peaked around the corresponding resonant layer. In linear theory the energy deposition profile in a loop that is excited with a uniform driving spectrum can be obtained by adding up the ~j2(r)-profiles for each of the driving frequencies. Poedts et al. [13] studied the localization of the intrinsic dissipation for a given equilibrium and for given modenumbers m and n by adding up the contributions
cally localized in the outer regions of the plasma column. The numerical survey carried out by Poedts et a!. [13] reveals that in most cases the i~J2(r)-profileis only slightly affected by the equilibrium parameters, at least as long as these parameters are varied within realistic bounds. On some special occasions, however, the heat is deposited more homogeneously as illustrated in fig. la, which shows the ~J2(r)-profiles for m = 1, n = 1, Jo = 0.1, d = 0.2, ~j = iO~, and v = 0.5 and 4. For the smooth current density profile (~ = 0.5) the Ohmic dissipation is almost uniformly spread over the coronal loop while for the peaked current density (v = 4) the intrinsic heating is localized in
82
S. Poedts et at
/
Numerical simulation of thestationary state of coronal ioops
the outer plasma layers. This is because for v = 0.5 the profile of the local Alfvén frequency is monotonic and the chosen driving frequencies (equidistantly spread over [Wa,,~,, toam,,j) yield resonant surfaces everywhere in the coronal loop, while for v = 4 the toa-profile is non-monotomc and the resonant surfaces are shifted to the outer regions of the plasma column (see fig. ib). 4.2. Efficiency Poedts et al. [12] noticed that not every driving frequency yields equally good coupling of the external source to the plasma. Extremely high efficiency is obtained when the driving frequency is close to the real part of the frequency of a so-called “collective mode” or “damped quasi-mode” which is the remnant of a discrete mode (in this case the stabilized external kink mode) that disappeared in the continuous spectrum [14,15]. For a bandwidth around the optimal driving frequency ía max(la(wp)}, and in the following the range of driving frequencies for which the ía 0.95 max{ la(COp)} is called the “efficient” bandwidth; clearly an important parameter for the efficiency of resonant absorption as a coronal loop heating mechanism. The efficiency of resonant absorption also depends strongly on the wavenumbers of the waves that are incident on the coronal loops as iilustrated in fig. 2, where the fractional absorption is shown versus the driving frequency for n = 5 and m = — 5—5, and where the other parameters are fixed (see the figure caption). Of course, as is clear from the definition of coa(r) (see above), the
f i .~
-‘7
0.81 0.6
1
~
°2L~~~~’ / ~ -
00
0.20.4
0.6
5 o~
Fig. 2. The fractional absorptionp versus the frequency of the driver for v=2, d=0.2, j,,=0.1, ,~=10_8,n=5, and m= —5, —3, —1,1,3 and 5.
range of the ideal Alfvén continuum — and so the range of frequencies that can be resonantly excited — depends on m. By consequence, many different driving frequencies contribute to the heating of the same plasma layers as the corresponding resonant layers range from r = 0 to r = 1 for each value of m. The frequencies in fig. 2 are dimensionless and related to physical frequencies a by the transformation formula i r —1 a(Hz)=0.1417co( ~“ \ 10 m / B (0) \ i p (0) - 1/2 x ~ doi T io’~ _3) m where we assumed a mean atomic weight of 0.6 (for fully-ionised H plasma ~i= 0.5; in the solar atmosphere the presence of extra elements makes ~ 0.6 [16]). Hence, for a loop with a radius of 5 x 106 m, a magnetic field of 0.005 T (50 Gauss) on the axis and a number density on the axis of 3, for example, the ideal Alfvén con1016 m tinuum reaches from 0.00896 Hz to 0.04009 Hz for m = —3, from 0.02689 Hz to 0.05345 Hz for m = 1, and from 0.04482 Hz to 0.06681 Hz for m = 5 (for the parameters used in fig. 2). So the range of frequencies at which this coronal loop can be resonantly excited agrees very well with the observed motions in the solar atmosphere which have periods of typically 100 s. The m = 1 mode indeed yields the strongest coupling and so the most efficient heating. For this mode ía = 1 for the optimal driving frequency and ía 90% for a large bandwidth. The m = —1 mode couples almost as efficient to the plasma as the m = 1 mode and this feature can be generalized: max( ía(top)} and the range of the efficient bandwidth are almost independent of the sign of the poloidal wavenumber. The fractional absorption gradually decreases as m increases and for m = 5, max( í (to )) amounts to 28%. For m = 3 the fracp tional aabsorption is 48% for the optimal driving frequency. When the wavenumber m is fixed and n is varied, it is found that the maximal value of ía
) (~
(for the optimal driving frequency) does not depend . on the value of n when m is positive. For negative values of m, however, the optimal efficiency is higher for the higher values of n [13].
S. Poedts et al.
/ Numerical simulation of the stationary state of coronal loops
OD 1.000
~
0 995
-~
o.goo
83
on the axis). For typical solar loop parameter values the energy supplied by the exciting waves is very efficiently dissipated and converted to heat.
/
5. Discussion and conclusions
0.98~ ,.—.~.—=.-=,= 6 7 B
9
10
11
12
The heating of solar coronal loops by resonant absorption of Alfvén waves has been investigated Fig. 3. The Ohmic dissipation rate versus the resistivity (in . . arbitrary units, chosen such that OD 1 for Rm 1012) for in the framework of hneanzed compressible restscop 0.195, m 2, a =1, q 0 Q5, p =1 and d =1. tive MHD. The resonant absorption of the waves that are incident on the coronal loops is numerically simulated in straight cylindrical loop models, Figure 3 shows how strikingly insensitive the externally excited by a periodic source. The s/aOhmic dissipation rate in the stationary state is to tionary state of this driven system and the Ohmic the electric conductivity of the plasma: enlarging dissipation rate in this state are determined by the magnetic Reynolds number Rm (m 1) by six means of a very accurate code based on the finite orders of magnitude changes the Ohmic dissipaelement technique which is extremely suitable for tion rate (in dimensionless units) by less than this purpose because of the strong localization 1.5%. Moreover, the energy deposition rate (due to the nearly-singular and boundary layer reaches an asymptotic value for Rm> 108. Since behaviour) of the solutions. We followed Grossthe resistivity of coronal loop plasmas is very mann and Smith [7] and computed the “intrinsic” small (i~ 10 8), the dissipation rate is almost dissipation; i.e. the power spectrum of the external independent of i~ for the relevant values of this source that excites the magnetic loops in the solar parameter. However, the heating rate due to resocorona is not specified in the numerical simulanant absorption does depend on the other equitions (it is assumed to be uniform). The actual libnum quantities. We did not succeed to confirm heating of a coronal loop can only be determined Grossmann’s and Smith’s [7] expression for the when the driving spectrum of the loop is known. ideal decay rate (their eq. (47)) which should corHowever, our numerical simulations give a very respond to the heating rate in a dissipative system. good picture of the efficiency of the resonant Moreover, we couhi not find an alternative forabsorption process in solar loops and yield valmula to express the Ohmic dissipation rate in uable information on the viability of this process those equilibrium parameters. It seems there are as a coronal heating mechanism. one or more parameters missing. One such paramThe energy supplied by incident waves on eter could surely be the proximity of the driving coronal loops, with frequencies within the range of frequency to the frequency of the “collective the Alfvén continuum, is very efficiently dismode”. In an attempt to elucidate the effect of the sipated in typical solar coronal loops. However, equilibrium parameters on the efficiency of resothe role of resonant absorption in coronal heating nant absorption in solar ioops, Poedts et al. [13] depends dramatically on the magnitude of the computed la(c~p)-profilesfor different choices of background magnetic field. Indeed, the energy the equilibrium parameters. They found that the dissipation rate due to resonant absorption is pro3 (as can be found by transforefficiency of the mechanism, by portional max{ la(t0p)} andheating by the range of measured the efficient ming backtotoB2(0) quantities with dimensions (see seebandwidth of driving frequencies, is substantially tion 2)) and so the magnitude of the background influenced by the parameter d (the plasma density field plays a crucial role in the amount of energy of the loop surface) as wall as by the parameters v that is converted into heat per unit of time. This (related to the shear of the equilibrium magnetic result is very promising as the observed magnetic field (see section 2.2)) and J 0 (the current density fields are rather strong in the solar corona. Other Log(R4
.
=
=
=
=
=
84
S. Poedts et al.
/
Numericalsimulation of the stationary state of coronal loops
indications in favour of resonant absorption are the localization of the heat deposition and the fact that the efficiency of the heating mechanism is independent of the length of the coronal loop. The intrinsic dissipation is localized in the outer plasma layers indicating that the absorption process is most efficient in the outer regions of the coronal loops, consistent with the temperature profiles ohserved in some loops [17]. Parker [18] remarked that the surface brightness of the active regions is almost independent of the dimensions. The only way where the dimensions of the coronal loop enter the numerics in our work is via the quantization factor k in the Fourier decomposition in the longitudinal direction (eq. 6). Given the fact that longer loops have larger radii and so the aspect ratio (and so k) is approximately the same for all coronal loops, the process of resonant absorption is equally efficient in all coronal loops to matter what their dimensions are. We conclude that resonant absorption is a viable candidate heating mechanism for coronal loops.
Acknowledgements S. Poedts is Research Assistant of the Research Council of the K.U. Leuven whose support is greatly acknowledged. We thank J.P. Goedbloed, J.V. Hollweg, B. Roberts and S. Antiochos for stimulating discussions and suggestions. The K.U. Leuven—MPIPP collaboration is supported by the European Economic Community Scientific Cooperation Contract No. ST2J-0275-C.
References [1] R. Rosner, L. Gollub and G.S. Vaiana, Ann. Rev. Astron. Astrophys. 23 (1985) 413. 12] J.L. Linsky, Sol. Phys. 100 (1985) 333. [3] C. Jordan and J.L. Linsky, in: Exploring the Universe with IUE Sateffite, ed. Y. Kondo (Reidel, Dordrecht, 1987)the p. 259. [4] J.A. lonson, Astrophys. J. 226 (1978) 650. [5]
161
J.V. Hollweg and G. Yang, J. Geophys. Res. 93 (1988) 5423. S. Poedts, W. Kerner and M. Goossens, J. Plasma Phys.
42 (1989) 27. [7] W. Grossmann and R.A. Smith, Astrophys. J. 332 (1988) 476. [8] G. Strang and G.J. Fix, An Analysis of the Finite Element Method (Prentice-Hall, Englewood Cliffs, NJ, 1973). [9] K. Appert, D. Berger, R. Gruber and J. Rappaz, J. Cornput. Phys. 18 (1975) 284). [10] J. Rappaz, Num. Math. 28 (1975) 15. [11] W. Kerner, K. Lerbinger, K. Gruber and T. Tsunematsu, Comput. Phys. Cominun. 36 (1985) 225. [12] S. Poedts, M. Goossens and W. Kerner, Sol. Phys. 123 (1989) 83. [13] S. Poedts, M. Goossens and W. Kerner, Astrophys. J., in press. [14] D. Van Eester, M. Goossens and S. Poedts J. Plasma Phys, submitted. [151 J.P. Goedbloed, Physica D 12 (1984) 107. [16] E.R. Priest, Solar Magnetohydrodynaniics (Reidel, Dord1984). Sol. Phys. 43 (1975) 327. [171 recht, P.V. Foukal, [18] E.M. Parker, Heating of the stellar corona, in: Coronal and Prominence Plasmas, ed. AL. Poland, NASA-CP-2442 (1986).