Magnetohydrodynamic Turbulence of Coronal Active Regions and the Distribution of Nanoflares

We present results from numerical simulations of an externally driven two-dimensional magnetohydrodynamic system over extended periods of time, used to model the dynamics of a transverse section of a solar coronal loop. A stationary forcing was imposed to model the photospheric motions at the loop footpoints. After several photospheric turnover times, a turbulent stationary regime is reached that has an energy dissipation rate consistent with the heating requirements of coronal loops. The turbulent velocities obtained in our simulations are consistent with those derived from the nonthermal broadening of coronal spectral lines. We also show the development of small scales in the spatial distribution of electric currents, which are responsible for most of the energy dissipation. The energy dissipation rate as a function of time displays an intermittent behavior, in the form of impulsive events, that is a direct consequence of the strong nonlinearity of the system. We associate these impulsive events of magnetic energy dissipation with the so-called nanoflares. A statistical analysis of these events yields a power-law distribution as a function of their energies with a negative slope of 1.5, consistent with those obtained for flare energy distributions reported from X-ray observations. A simple model of dissipative structures, based on Kraichnan's theory for MHD turbulence, is also presented.


INTRODUCTION
Various scenarios for the heating of coronal loops in active regions have been proposed, most of which have in common the formation and dissipation of highly structured electric currents. The spontaneous formation of tangential discontinuities (Parker the development of an 1972, 1983), energy cascade driven by random footpoint motions on a force-free conÐguration Ballegooijen or the (van 1986), direct energy cascade associated with a fully turbulent magnetohydrodynamic (MHD) regime & Priest (Heyvaerts & Ferro Fonta n are a few examples. 1992 ; Go mez 1992) The main conclusion is that the formation of small scales in the spatial distribution of electric currents is necessary to enhance magnetic energy dissipation and therefore provide sufficient heating to the plasma conÐned in these loops.
The dynamics of a coronal loop driven by footpoint motions can be described by MHD equations. Since the kinetic (R) and magnetic (S) Reynolds numbers in coronal active regions are extremely large (R D S D 1010h12), we expect footpoint motions to drive the loop into a strongly turbulent MHD regime. Footpoint motions, whose length scales are usually much smaller than the loop length, will cause the coronal plasma to move in planes perpendicular to the axial magnetic Ðeld, generating small transverse components of the magnetic and velocity Ðelds. The nonlinear dynamics of these Ðelds can be adequately described using two-dimensional MHD equations. In this paper we present numerical integrations of these equations, considering the driving action of footpoint motions on a generic transverse section of a loop. In we model this coupling, and in°2°3 we describe the numerical technique used for the integration of the two-dimensional MHD equations. In we°4 report the energy dissipation rate and turbulent or excess velocity that we obtain, and describe the temporal series of magnetic and kinetic energies. The energy power spectra computed from our simulations are shown in The°5. spatial structures developed by this system, such as the formation of small scales in the distribution of electric currents, are presented in°6.
In we perform a statistical analysis of the dissipation°7 events and study the correlations between di †erent variables. In we discuss some of the implications of our°8 results. For instance, in we present a simple model for°8.1 the dissipative structures of a turbulent MHD system, and compare the predictions arising from this model with the results obtained from our simulations. In we discuss°8.2 the relevance of heating theories for the thermal balance of active regions, which is globally characterized by the scaling laws. Di †erent theories give rise to di †erent relationships between parameters of the active region loops, such as the pressure, temperature, magnetic intensity, and length. The possible connection between the dissipation events studied in this paper and the much larger ones involved in solar Ñares is discussed in Finally, in we summarize the°8.3.°9 main conclusions of this paper.

EQUATIONS FOR TWO-DIMENSIONAL MAGNETOHYDRODYNAMICS
where is the Alfve n speed, l is the kinev A \ B 0 /(4no)1@2 matic viscosity, g is the plasma resistivity, t is the stream function, and a is the vector potential. The Ñuid vorticity is w \ [+2t, j \ [+2a is the electric current density, and [u, v] \ z AE $u Â $v. For given horizontal photospheric motions applied at the footpoints (plates z \ 0 and z \ L ), transverse velocity and magnetic Ðeld components develop in the interior of the loop, given by and The RMHD equations can be regarded as describing a set of two-dimensional MHD systems stacked along the loop axis and interacting among themselves through the terms (see, e.g., & van Hoven 1996). plicity, hereafter we study the evolution of a generic twodimensional slice of a loop. We therefore model the v A L z terms given in equations as external forces. To this (1)È(2) end, we assume the vector potential to be independent of z and the stream function to interpolate linearly between t(z \ 0) \ 0 and t(z \ L ) \ (, where ((x, y, t) is the stream function for the photospheric velocity Ðeld. These assumptions yield (in These approximate expressions correspond to 0 eq. [2]). an idealized scenario in which the magnetic stress distributes uniformly throughout the loop. The two-dimensional equations for a generic transverse slice of a loop become where the external forcing f relates to the photospheric motions through This forcing can also be interf \ (v A /L )(. preted as an electric Ðeld generated by an electrostatic potential di †erence between the two photospheric plates (z \ 0 and z \ L ). We choose a forcing term that is narrowband in wavenumber space and stationary in time, as described in the next section. A comparable approach was used by et al. although with slight di †er-Einaudi (1996), ences in the forcing term (see also Velli, & Georgoulis, Einaudi 1998 For the nonlinear terms, we Ðrst compute the spatial derivatives in k-space. Next, we perform inverse FFTs and obtain the products required for the Poisson brackets in real space. Finally, we perform direct FFTs to obtain the Poisson brackets in k-space. To eliminate the aliasing in the transform operations, we make zero the Fourier components of modes with k [ N/3, where N is the number of grid points in each direction. The temporal integration scheme is a Ðfth-order predictor/corrector, in order to achieve an almost exact energy balance over our extended time simulations (about 250 turnover times). We model the external forcing in as equation (5) f k with to simulate the action of photospheric f 0 \ constant, granular motions on the loop Ðeld lines. We chose a narrowband and nonrandom forcing to ensure that the broadband energy spectra and the signatures of intermittency that we obtained (see below) are exclusively determined by the nonlinear nature of the MHD equations.
We turn equations into a dimensionless version, (3)È(4) choosing l and L as the units for transverse and longitudinal distances, respectively. We choose as the unit for v 0 \ f 0 1@2 velocities, since the Ðeld intensities are determined by the forcing strength. Since and The dimensionless dissipation coefficients are l 0 \ l/(lf 0 1@2) and g 0 \ g/(lf 0 1@2).

ENERGY DISSIPATION RATE
To restore the dimensions to our numerical results, we used typical values for the solar corona : L D 5 ] 109 cm, l D 108 cm, cm s~1, G, and n D 5 ] 109 v ph D 105 B 0 D 100 cm~3. In addition, the values for the dissipation coefficients allowed by our moderate-resolution code are l 0 \ g 0 \ 7 ] 10~3.
shows magnetic and kinetic energy versus time. Figure 1 After about Ðve photospheric turnover times (q ph \ l/v ph1 03 s), a statistically steady state is reached. The behavior of both time series is highly intermittent despite the fact that the forcing is constant and coherent. Intermittent signals display activity during only a fraction of the time, which decreases with the scale under consideration. To quantify the degree of intermittency in a time series S(t) (where S(t) can be any relevant physical quantity, such as the magnetic or kinetic energy or the dissipation rate), we adopt the method described by We compute a high-Frisch (1996). pass Ðltered version of S(t), where is the Fourier transform of S(t) and is the Ðlter SOE w w c frequency. The Ñatness of the Ðltered signal is deÐned as The characteristic feature of intermittent signals is that their associated Ñatness (as deÐned in grows eq. [9]) without bound with the Ðlter frequency. This is indeed the case for our time series of total energy and dissipation rates, up to the numerical resolution of the Ðlter frequency. The Ñatness is a measure of the departure from Gauss-F(w c ) ianity of the signal S(t) for a Gaussian [F(w c ) \ const. \ 3 signal, since the Ðltering process is a linear operation, which therefore does not alter its Gaussian character].
This kind of behavior is usually called internal intermittency, to emphasize the fact that the rapid Ñuctuations are not induced by an external random forcing. The forcing used by et al. also displays a slow time Georgoulis (1998) variation, but they introduced randomness in space by changing the places of the forcing vortices. Notwithstanding, many of their results are qualitatively consistent with the ones presented here, since spatial randomness is not an essential feature for the behavior of the average quantities.
also shows that the kinetic energy always Figure 1 remains at a small fraction of the magnetic energy in the stationary regime.
The energy dissipation rate is also a strongly intermittent quantity, as shown in For turbulent systems at Figure 2. large Reynolds numbers, the dissipation rate in the station-ary regime is expected to be independent of the Reynolds number Re For the rather moderate (Kolmogorov 1941 1997). addition, the dissipation rate is not completely independent of the spatial dependence of the forcing term. For instance, in our simulation only Fourier modes satisfying the ring condition 3 \ kl \ 4 are externally driven. If we also drive the modes k \ l~1(^3, 0), k \ l~1(^4, 0), k \ l~1(0,^3), and k \ l~1(0,^4), lying at the border of the ring, then the dissipation rate is larger and the behavior of the time series seems less intermittent. Since these modes correspond to purely one-dimensional spatial patterns, thus departing from the assumed isotropy for the photospheric velocity Ðeld, we exclude them from the forcing. The total dissipation rate is given by v^1.7 ] 1024 ergs s~1 where we have also indicated the scaling with the relevant parameters of the problem. Using an Alfve n time q A \ L /v A and a photospheric time we can rewrite the dis- The relevance of this functional dependence of the heating rate upon the physical and geometrical parameters of the system will be discussed in°8.2. We can transform the heating rate into an energy inÑux from the photosphere by simply dividing by twice the transverse area (because we have two boundaries), i.e., F \ v/ [2(2nl)2]. Note that the energy balance between the photospheric inÑux and Joule (and viscous) dissipation is a dynamic process. However, once the stationary turbulent regime is reached, the time averages of these energy rates (over timescales on the order of or larger than i.e., disre-q ph , garding the intermittency structure) are equal.
shows the quantitative value of the energy Equation (12) inÑux as well as the explicit dependence with the relevant parameters of the problem, This energy Ñux compares quite favorably with the heating requirements for active regions, which span the 1977).
From our simulations we can also obtain the "" excess ÏÏ velocity associated with the observed line broadening of a number of X-ray spectral lines in solar active regions (Seely et al.
see also & Strong We estimate this 1997 ; Saba 1991). velocity as the root mean square value of the turbulent velocity Ðeld, which is proportional to the square root of the kinetic energy : v excess The value from our simulation and the explicit dependence with the parameters of the problem is v excess . (14) This turbulent velocity is within the range of nonthermal line broadenings of ultraviolet spectral lines as measured by Doschek, & Feldman (10È25 km s~1) and Cheng, (1991) ties probably correspond to hotter and larger active regions, with a magnetic topology much richer than the one implied by the RMHD equations used for the present simulations. Also note that our heating rate (see is still a eq. [12]) factor of 5 smaller than the energy requirements for the hottest active regions. As discussed by et al. Georgoulis it is likely that in order to heat the largest active (1998), regions, full three-dimensional models are required.

ENERGY POWER SPECTRA
In spite of the narrow forcing and even though the velocity and magnetic Ðelds are initially zero, nonlinear terms quickly populate all the modes across the spectrum. The total energy (magnetic plus kinetic) power spectrum is plotted in The spectra correspond to di †erent corresponds to a Kraichnan spectrum In turb-(E k D k~3@2). ulent regimes, the e †ect of nonlinearities is to redistribute excitations in Fourier space in a virtually stochastic fashion.
Since dissipative e †ects are only nonnegligible at very large wavenumbers, a net Ñow of excitations toward large wavenumbers is established to compensate for the dissipation. Energy is therefore injected into the system at low wavenumbers, cascades toward large wavenumbers in the socalled inertial range, and is efficiently quenched in the dissipative range. The Kraichnan spectrum is expected to be FIG. 5.ÈZones of intense magnetic dissipation. White regions correspond to zones that concentrate 60% of the total dissipation, at time t \ 11 ] 103 s. satisÐed at the inertial range, where the role of external forcing and dissipation are negligible. For moderateresolution simulations like the one presented here, the inertial range is rather limited, gradually entering into the dissipative range at k [ 10, as shown in Figure 3.
In we plot the magnetic and kinetic spectra, Figure 3b showing that the spectra become close to equipartition when approaching the dissipative range (large k). (E k M D E k K) Nonetheless, note that the total kinetic energy (i.e., integrated over wavenumbers) remains much smaller than the total magnetic energy, as evidenced by Figure 1.

DISSIPATIVE STRUCTURES
shows the spatial distribution of electric Figure 4a current density for Intense positive currents are t \ 22q ph . shown by white regions, while intense negative current concentrations are indicated in black. The magnetic Ðeld is shown by arrows. This picture is a typical example of a strongly turbulent regime, with current concentrations of both O and X types. At the left center of the box there is a typical current sheet, of the type that result after the collapse of an X-type point. Moreover, the velocity Ðeld and vorticity displayed in conÐrm the standard Figure 4b picture of magnetic reconnection in current sheets, with Alfve nic jets emerging at the sides of the sheet and a quadrupolar distribution of vorticity. At the same time, Figures 4a  and show that in turbulent regimes there is also a rich 4b variety of dissipative structure of di †erent morphologies and degrees of complexity, in addition to the highly symmetric current sheet that we have just described. Figure 5 shows the spatial distribution and morphology of the most intense dissipative structures. The dissipation rate inside these structures is 60% of the total, although the Ðlling factor is only 6%. This highly inhomogeneous distribution of energy dissipation is a manifestation of the spatial intermittency that is characteristic of turbulent regimes.
The thickness of current sheets is determined by the magnetic Reynolds number. In numerical simulations of turbulent regimes, the Reynolds number is usually made as large as possible, so that the sheet thickness (which is the smallest feature expected in these simulations) is just marginally resolved. The widths of these current sheets, on the other hand, are determined by the dynamics, and therefore a broad distribution of widths can be observed.

DISTRIBUTION OF EVENTS
proposed that the energy dissipation of the Parker (1988) stressed magnetic structures takes place in a large number of small events, which he termed "" nanoÑares.ÏÏ The superposition of a large number of such events would give the global appearance of a spatially homogeneous and stationary heating process. From a turbulent scenario, it seems quite natural to relate this spiky (both in space and time) heating to the internal intermittency present in all turbulent regimes.
We therefore associate the peaks of energy dissipation displayed in with the so-called nanoÑares. We esti- Figure 2 mate the occurrence rate for these nanoevents, i.e., the number of events per unit energy and time P(E) \ dN/dE, so that is the total number of events per unit time and is the total heating rate (in ergs s~1) contributed by all events in the energy range A simple inspection [E min ; E max ]. of the v(t) time series shown in indicates that these Figure 2 events are in a highly concentrated or piled-up regime, i.e., that their event rate R multiplied by their typical duration is much larger than unity. In this regime, many dissipation events are going o † simultaneously at any given time. It is therefore impossible to perform a statistics of events, since we are unable to separate them.
As a way out to this difficulty, we deÐne an event in the following fashion : Ðrst we set a threshold heating rate on v 0 the time series displayed in on the order of its time Figure 2, average ; "" events ÏÏ are then excesses of dissipation that start when v(t) surpasses and Ðnish when v(t) returns below v 0 v 0 . Once a particular threshold is set, we perform a statistical analysis of the events, keeping track of their peak values, durations, and total energy content. The implicit assumption behind our working deÐnition is that the small fraction of events that emerge over the threshold are statistically representative of the whole set.
The occurrence rate as a function of energy (see Fig. 6a) displays a power-law behavior, in the energy range spanning ergs to E min^1 025 E max^2 ] 1026 ergs.
We also computed the distribution of events as a function of peak Ñuxes, which is a power law with slope a P \ 1.7 0.3, as shown in The slope we obtained is arising from di †erent choices for the threshold (which v 0 varies within a range such that the total number of events is maximum) is always smaller than the error sources mentioned above. A more technical procedure for choosing a threshold, which introduces a Ðtted s2 distribution to eliminate the noise, has been used in et al. Georgoulis (1998). This Ðtted threshold is slightly higher than the average value of the series.
Other important results from this statistical analysis are the correlations between the di †erent parameters of these events. In there is a scatter plot of energy released Figure 7a versus event duration. The data can be Ðtted by a power law where This result is consistent E D qcEq, c Eq \ 2.02^0.02. with the correlation reported by Petrosian, & McTier-Lee, nan from hard X-ray observations. The duration (1993) versus peak correlation is plotted in which can be Figure 7b, Ðtted by where These two corre-q D PcqP, c qP \ 1.12^0.03. lation indexes agree with the approximate relationship which should be exact if the dimensional relationship E D qP holds. Moreover, if these correlations are signiÐcant, then the di †erent power-law indexes and (a E , a P , a q ) should be mutually related through  (1998). We make the simplifying assumption that all the energy dissipation takes place inside current sheets of variable width j and Ðxed thickness which is the dissipation l d , length scale for this turbulent regime. The current sheets are formed between vortices of size j, and therefore their width is also j, which we assume as a free variable in order to perform a statistical analysis of dissipation events. We also assume that the energy spectrum follows a Kraichnan power spectrum, i.e., (see, e.g., & Welter Biskamp 1989) where v is the dissipation rate and is the total magnetic v A 2 energy per unit of mass, and where From this spectrum, we can obtain a R \ (v A L M )/g. typical Ðeld (in units of velocity) in a vortex of size j as b j so this is the typical Ðeld at the sides of a current sheet of width j. In the present model, we assume a hierarchy of dissipative structures labeled by their values of j, which are spatially distributed in such a way that their average separation is also on the order of j. The dissipation rate contributed by all the structures of width j is therefore given by since is the area Ðlling factor of current sheets of (l d j/j2) width j, and the current density is approximately j j D b j /l d . The lifetime of these current sheets is determined by the nonlinear timescale in which the vortices of size j break down and transfer their energy to smaller vortices as part of the turbulent cascade. This timescale is much shorter than the dissipation time of the current sheets, thus implying that two interacting vortices break down long before their energy is fully dissipated in the current sheet formed in between. The nonlinear timescale, following the Kraichnan model of MHD turbulence, is In this stationary turbulent regime, the number of ongoing dissipation events of width j per unit time is The energy dissipated in a single structure (current sheet) of size j is then Combining equations and we derive the dis-(27) (28), tribution of events as where which readily satisÐes a E \ 9/4 \ 2.25, where and . Therefore, according to this very simple model, the slope of the distribution of events is larger than 2, implying that small events dominate the heating process. There are several possible causes for this disagreement with the slope obtained from our simulation. (1) The theoretical model assumes that all the energy dissipation takes place in sheet structures formed by the coalescence of two vortices. Although our simulation shows that a fair fraction of the energy dissipation takes place in such structures, it is also clear that other structures, such as O-points, also contribute to the dissipation. (2) We assume that all current sheets have thickness which is likely to be an oversimpliÐcation. l d , Depending on the parameters and boundary conditions under which current sheets are formed, di †erent regimes of magnetic reconnection might occur, corresponding to different thicknesses. (3) To associate the lifetime of current sheets with the nonlinear timescale is also a simpliÐcation, since other processes might contribute to disrupt the current sheets once they are formed. (4) It might also be that our simulations with modest spatial resolution do not reÑect the statistical properties of dissipative structures in a large Reynolds number regime. In addition, because of the characteristics of the simulations presented here, the current sheets are resolved only marginally. Simulations with higher resolution might provide a better understanding of the dynamics of these current sheets, which in turn will help to generate a more complete theoretical model.
From equations and the correlation between (28) (26), duration and dissipated energy becomes where c Eq \ 4. The total event rate is then To evaluate whether we are in a pile-up regime of dissipation events, we can compare the total event rate R with the shortest duration events (from eq. [26]) : We obtain Rq ld \ R4@3 ? 1 .
Therefore, in a fully turbulent regime, the corresponding dissipation events always pile up.
In summary, from this rather simple model we can derive the following results : (1) a slope for the distribution of events, (see (2) a power-law relationship a E \ 2.25 eq. [29]) ; between total dissipated energy and duration, (see c Eq \ 4 and (3) a strong pile-up scenario for dissipation eq. [31]) ; events.

Heating Rate and Scaling L aws
It is interesting to compare the heating rate scaling shown in with the one arising from equation (11) (°4) ParkerÏs model in which the energy dissi-(Parker 1972), pated is assumed to match the energy injected by photospheric motions. When a Ðeld line is tilted by slow motions at its footpoints, a small transverse magnetic Ðeld component b is generated. The small tilt angle for a given Ðeld line after an elapsed time q is which is in turn of the b/B 0 , order of Therefore, Assuming that v ph q/L . b^B 0 v ph q/L . the energy accumulated over the time q completely dissipates, we obtain According to footpoint motions pump Parker (1983), energy into the coronal magnetic Ðeld until the angle between adjacent Ðeld lines reaches a critical value At a * . this point, the built-up energy is suddenly released. Therefore, since the heating rate per unit volume a * \ v ph q/L , scales as From our it emerges that q is an equation (11) (°4), average between the Alfve n and photospheric timescales It is interesting to note that the timescale [q^(q A q ph )1@2]. is also the growth time for nonlinear insta-q^(q A q ph )1@2 bilities in an RMHD system De Luca, & McCly-(Go mez, mont which are responsible for the energy transfer in 1993), k-space. Thus, the heating rate per unit volume, according to the present model, scales as The particular dependence of the heating rate of a loop on its magnetic Ðeld and its length L is important, since B 0 the heating rate in turn relates to the thermodynamic variables of the loop (such as pressure and temperature) through the so-called scaling laws Tucker, & (Rosner, Vaiana The scaling laws are derived assuming a static 1978). balance between a uniform heating rate, thermal conduction, and (optically thin) radiative losses. In addition, the gas pressure is assumed to be uniform along the loop, which is a good approximation for loops with heights smaller than 105 km. The relationships between the heating rate H of a loop of length L with the thermodynamic variables P (pressure) and T (temperature) can be quickly derived by assuming that heating, conductivity, and radiative losses are all comparable. Therefore, where is the Spitzer conductivity and is a i 0 ((T )^( 0 T~b power-law Ðt for the function of radiative losses for temperatures ranging from 106 to 5 ] 106 K.    (16) implies that the contribution to energy dissipation in a given energy range is dominated by the most [E min ; E max ] energetic events (i.e., According to this result, the E^E max ). relatively infrequent large-energy events contribute more to the heating rate than the much more numerous smallenergy events. Therefore, we speculate that the heating rate of a given active region is to some extent determined by the magnetic topology of that region. However, it is important to note that this is just one of the two possible scenarios for coronal heating.
pointed out that this regime Hudson (1991) (dominated by large events) would reverse if there were a turn-up of the slope to values larger than 2 toward the low-energy end. & Trottet recently reported Mercier (1997) indirect evidence for such a turn-up from type I starburst observations. et al. presented cellular Vlahos (1995) automata simulations with an anisotropic rule for the redistribution of the magnetic Ðeld, and reported an index of [3.5 for the peak Ñux distribution. The model presented in based on KraichnanÏs theory for stationary MHD°8.1, turbulence, also yields an index larger than 2. Therefore, it seems apparent that to elucidate whether large or small events are dominant, more observational and theoretical e †orts are required.

CONCLUSIONS
In the present paper we simulate the dynamics of a transverse section of a solar coronal loop through an externally driven two-dimensional MHD code. The relevant results of this study are as follows.
1. For an external forcing that is narrowband in wavenumber, we Ðnd that the system becomes strongly turbulent, and after about Ðve photospheric turnover times a stationary turbulent regime is reached.
2. The energy dissipation rate obtained for typical footpoint velocities is consistent with the power necessary to heat active region loops (F^2.1 ] 106 ergs cm~2 s~1).
3. The energy dissipation rate displays a highly intermittent behavior, which is a ubiquitous characteristic of turbulent systems. Temporal intermittency manifests in the form of discrete-like dissipation events in the time series of the total energy dissipation rate. Spatial intermittency is also apparent, since a large fraction of the dissipation at any given time takes place in a rather small fraction of the volume. 4. A statistical analysis of dissipation events performed on a long-term numerical simulation (about 250 turnover times, which is roughly 3 days in real time) shows a power-law event rate, proceeding as dn/dE D E~1.5, that is remarkably consistent with the statistics of Ñare occurrence derived from observations et al. (Hudson 1991 ;Crosby et al. 1993 ;Lee 1993 ;Shimizu 1995). 5. The peak energy release and duration of events also follow power-law distributions. Furthermore, all these quantities are mutually correlated by power-laws.
We would like to sincerely thank to our referee, Giorgio Einaudi, for his constructive comments, which contributed to improve an earlier version of this manuscript. We acknowledge Ðnancial support by the University of Buenos Aires (grant EX247), Fundacio n Antorchas, and NASA (grant NAGW-4644).