A SIMPLE MODEL FOR THE EVOLUTION OF MULTI-STRANDED CORONAL LOOPS

EUV and soft X-ray observations of the solar corona reveal a high degree of structuring and intermittency, especially in and around active regions (ARs). Due to the frozen-in condition of the coronal plasma, the emitting material is ordered by the magnetic field in elongated features commonly known as coronal loops. The study of the structure and evolution of loops is therefore fundamental for understanding coronal dynamics, particularly in relation to the origin of coronal heating (Mandrini et al. 2000). Historically, the first attempts to physically describe coronal loops began with what is known as the Rosner et al. (RTV, 1978) approach, which considers loops as monolithic structures containing plasma in static equilibrium. Although soft X-ray loop observations seem to be consistent with the RTV assumptions (Porter & Klimchuk 1995), more recent observations from the Transition Region and Coronal Explorer (TRACE) suggest that EUV loops are actually too dense to be in static or quasi-static equilibrium (Aschwanden et al. 2001; Winebarger et al.2003). The properties of these loops are more consistent with the predictions of models based on impulsive heating (Winebarger et al. 2003; Klimchuk 2006).

TRACE observations also show a high degree of filamentation of the coronal structure, suggesting that what appear to be monolithic loops can actually be collections of unresolved individual strands with more or less independent evolutions. Although the question of whether coronal loops are isothermal or multithermal has been a matter of heated debate for several years (Aschwanden & Nightingale 2005; Schmelz & Martens 2006), recent observations seem to indicate that both kinds of loops actually occur (Schmelz et al. 2009). Combining all the available information, there are strong arguments in favor of many loops being made of unresolved strands that evolve more dynamically than quasi-static models predict (for a recent review on the subject see Klimchuk 2009).

There are two main families of mechanisms proposed to explain the origin of coronal heating: the so-called AC and DC heatings. Generally speaking, AC refers to the dissipation of MHD waves, while DC relates to the dissipation of highly stressed magnetic fields. The latter involves the release of free magnetic energy at locations of strong gradients in the magnetic field, i.e., current sheets. Ultimately, the source of energy is the distortion of the magnetic structures by convective motions. This is precisely the mechanism envisioned by Parker (1988). He proposed that, since the footpoints of magnetic strands are being displaced by photospheric granular motions, there must be a continuous increase of the magnetic stress between neighbor strands. These stresses accumulate until some threshold is reached and reconnection occurs. In the process, stored magnetic energy is released and the plasma is heated. Because of similarities with the usual flaring process and the energy scale involved, these reconnection events are usually called nanoflares (Cargill 1994; Cargill & Klimchuk 2004). One key point of the process described above is the existence of a threshold or critical value for the reconnection to occur. A theoretical support for the existence of this threshold is given in Dahlburg et al. (2005, 2009). They found that when a certain field misalignment is reached, the secondary instability provides the mechanism for a rapid energy release.

The intermittent and filamentary behavior, and the possible presence of power-law distributions (see e.g., Aschwanden & Parnell 2002), led some authors to propose that the solar corona might be a good example of a self-organized critical (SOC) system (Bak 1990). Beginning with the seminal work of Lu & Hamilton (1991), there have been a number of important studies involving the use of SOC models in relation to the coronal heating problem (for a review on SOC models see Charbonneau et al. 2001).

In our recent paper, López Fuentes et al. (2007), we studied the evolution of coronal loops observed with the Solar X-ray Imager (SXI) on board the Geosynchronous Operational Environmental Satellite 12 (GOES-12). Since this instrument observes the Sun almost continuously, it allowed us to follow the evolution of a set of coronal loops from birth to disappearance. We found that the loop evolution can be separated into three main phases, which we called rise, main, and decay. As we indicated there, and began to explore in Klimchuk et al. (2006), the observed evolution can be related to the development, maintenance, and destruction of a critical system. Here, we develop a simple model based on a cellular automaton (CA) that reproduces the observed evolution. We model the plasma response to the heating using the EBTEL (Enthalpy-Based Thermal Evolution of Loops) hydrodynamics code (Klimchuk et al. 2008). We then create synthetic light curves taking into account the temperature-dependent response of the SXI instrument. We compare these with the observed light curves and analyze how the different parameters of the model influence the properties of the obtained light curves, such as the characteristic intensities and timescales, which we studied in detail in López Fuentes et al. (2007) for the observed cases.

In Section 2, we summarize the results from López Fuentes et al. (2007) and show examples of observed light curves that will be later compared with the model. In Section 3 we describe the CA model, and in Section 4 we analyze our results. We discuss and conclude in Section 5.

In López Fuentes et al. (2007), we studied the evolution of a set of coronal loops observed with the SXI on board GOES-12 (Hill et al. 2005). This instrument has the advantage of observing the Sun almost continuously, allowing us to follow the whole evolution of loops. Previous similar studies based on Yohkoh/SXT or TRACE observations were limited by the occultation times of the satellites and times of interrupted observation due to high radiation levels (e.g., when passing through the South Atlantic anomaly), therefore reducing the continuous observation times to a small portion of the loop lifetime (see, e.g., Porter & Klimchuk 1995).

We analyzed a set of 457 images taken with the telescope in the filter position "open," which is most sensitive to emission from plasmas below 2 MK. The analyzed data span 3 days (2003 August 26–28) of observations with a cadence of approximately six images per hour. From these data we selected and tracked the evolution of seven coronal loops. Following a procedure that is thoroughly described in López Fuentes et al. (2007), we obtained light curves of the intrinsic (background-subtracted) intensity of each of these loops. We found that the loop evolution can be separated into three parts: a rise phase during which the loop intensity slowly increases, a main phase during which the intensity remains approximately constant (or its long-term variation is much less pronounced), and a decay phase during which the intensity decreases until the loop is almost indistinguishable from the background. In Figure 1 (upper panels) we reproduce adapted versions of two of the light curves studied in López Fuentes et al. (2007). In the original analysis, we applied a procedure to remove the background contribution that left a residual component of long-term variation (see discussion in Section 3.1 in López Fuentes et al. 2007). For an easier comparison with the model light curves, here we removed the residual background by subtracting in each case a linear function of time, so that the initial and final intensities of the loops are reduced nearly to zero.

Zoom In Zoom Out Reset image size

We found that the durations and characteristic timescales of the rise and decay phases are longer than the expected cooling times. We concluded that this is consistent with two alternative scenarios. In the first one, loops are monolithic and the heating source varies very slowly with time, so the plasma responds quasi-statically. In this way, the three evolution phases would reflect the phases of the long-term change of a slowly varying heating source. In the second scenario, loops are made of unresolved and independent magnetic strands which are heated by impulsive events such as nanoflares. Since the observed loop emission is the summed contribution of the constituent strands, the different phases would be directly related to the variation of the number and intensity of the energy dissipation events that heat the strands. We argued that under this particular scenario the observed loop evolution can be related to the development, maintenance, and destruction of a self-organized system. In the next sections, we present and analyze a model that reproduces precisely this kind of behavior.

3.1. General Description

In this section, we describe the construction of a CA model that reproduces the main features of the observed loop evolution. Motivated by the arguments given in Section 1, the model simulates the process described there regarding the scheme: motion of magnetic strand footpoints → increase of magnetic stress → energy release by reconnection. Figure 2 shows a cartoon that illustrates the main idea.

Zoom In Zoom Out Reset image size

Let us begin with a bunch of thin parallel magnetic flux tubes or strands that connect two different regions of the solar surface (panel (a)). For simplicity the strands are drawn straight, so the top and bottom planes correspond to two relatively distant photospheric regions of opposite magnetic polarity. As time goes on, strand footpoints are horizontally dragged by random photospheric motions, so the strands become progressively more misaligned. This produces growing magnetic stresses and intensifying current sheets in the region separating adjacent strands (panel (b)). For two given neighbor strands, the misalignment angle increases until a critical value is reached (i.e., the secondary instability develops) and the energy is released in the form of an impulsive event. Consequently, the field between the two strands relaxes. This sudden change in the field may occasionally cause new instability conditions between these strands and their other neighbors. Under certain conditions, instabilities propagate further, producing avalanches of events involving many strands. This is the typical situation considered in SOC models. For clarity, we adopt the following terminology. The basic interaction between two strands that occurs each time the critical value is reached is called a "reconnection event." The total release of energy into a single strand by its reconnection with one or more neighbors over a short period of time is called a "nanoflare." This follows the convention used by the loop modeling community.

The system continues to evolve in response to the footpoint driving until there is a rather abrupt appearance of a critical state in which the energy input by the photospheric motions is statistically balanced by the dissipation due to reconnection events. In the critical state, the total energy of the system fluctuates intermittently around a constant level. Local intermittency results from the occurrence of nanoflares of different sizes.

The idea is then to develop a numerical algorithm that quantitatively reproduces the above scheme. First of all, we need to translate footpoint displacements into magnetic field stresses. To accomplish this, we begin with the simple scenario shown in Figure 2(c). Let us concentrate on one particular strand and, for simplicity, suppose that only one footpoint (the lower one) is being displaced. If the particular strand has a vertical magnetic field strength B_v, any horizontal displacement of its footpoint will produce a new horizontal component of the field

$$\begin{equation} \delta B_{h}=B_{v}\tan \theta, \end{equation} \tag{ 1 }$$

where θ is the angle that the now inclined strand forms with the vertical direction. It is easy to see that if the strand has a length L and the footpoint has been displaced a distance d, then tan θ = d/L. Therefore,

$$\begin{equation} \delta B_{h} = B_{v}d/L. \end{equation} \tag{ 2 }$$

The distance d is a characteristic length scale of the photospheric motions, which can be associated with the sustained displacement of a strand footpoint during the turnover time of a convective granule (≈10³ s). Assuming a typical photospheric velocity v ≈ 1 km s⁻¹ results in d ≈ 1000 km. After a number of successive displacements, the strand would acquire a horizontal field

$$\begin{equation} B_{h} = \Sigma \delta B_{h}. \end{equation} \tag{ 3 }$$

At this point, we have assumed that B_h increases with each and every step. This is only true if the new step is in the same direction as the previous step or if the change in direction causes the strand to wrap around a neighbor. The latter will generally be true if the random walk step size is larger than the characteristic separation of neighbor footpoints. However, even then, some steps will cause the strand to backtrack and B_h will decrease. We treat this more realistic situation below.

As a consequence of the random walk, the horizontal field B_h will gradually increase until the critical condition referred to above is reached. If we call θ_c the critical angle of misalignment between two adjacent strands, it is easy to see that there is an associated critical value of the difference in their horizontal field components:

$$\begin{equation} B_{c}=B_{v}\tan \theta _{c}. \end{equation} \tag{ 4 }$$

Note that θ_c refers to the relative angle between the adjacent strands and not the inclination from vertical, θ, in Equation (1).

Assuming one-dimensional displacements (see Section 3.2), we relate the misalignment between two neighbor strands (with indices p and q) to the difference of their horizontal field components ΔB_h = B_h(p) − B_h(q). The fulfilment of the critical condition implies |ΔB_h|>B_c. If we assume that all the field is equally distributed after reconnection, the new (primed) horizontal field in each strand will respectively be

$$\begin{equation} B^{\prime }_{h}(p)=B_{h}(p)-\Delta B_{h}/2\; \hbox{and}\; B^{\prime }_{h}(q)=B_{h}(q)+\Delta B_{h}/2. \end{equation} \tag{ 5 }$$

The released energy corresponds then to the difference in the total magnetic energy before and after reconnection. It can easily be demonstrated that

$$\begin{equation} \Delta E=(E^{\prime }-E)=\left(\frac{1}{8\pi } \right) \frac{(\Delta B_{h})^{2}}{4} \end{equation} \tag{ 6 }$$

per unit volume, where we have assumed that the energy is divided equally between the plasma in the two strands.

3.2. Model Implementation

For the numerical implementation of the model we use a one-dimensional array of N_m points. Each point corresponds to a single magnetic strand and, initially, all strands have horizontal magnetic field B_h = 0. At each time step, a small random amount of field δB_h is added to each site according to a statistical distribution that simulates the photospheric displacements of strand footpoints. Although the field increase in the model is one dimensional, we apply a distribution based on two-dimensional footpoint displacements. The idea is illustrated in Figure 3(a). If one compares the position of a strand footpoint at time step (t − 1) with its position at time step t (gray arrow), the probability that during the following time step (t + 1) (black arrows) the footpoint moves farther away increasing the horizontal magnetic strength is approximately 3/4, while the probability that it moves back on its previous track, canceling (or almost canceling) the previous displacement, is roughly 1/4. Therefore, we add the δB_h values according to the distribution shown in Figure 3(b). Values from the Gaussian distribution centered on d₀ are three times more probable than values given by the Gaussian centered on −d₀. The distance d₀ = 1000 km is the average horizontal displacement of a strand footpoint during the turnover time of a granular cell (in our scheme, the duration of each time step). The width of the Gaussians is 20% of d₀. As explained in the previous section, the random displacements taken from this distribution translate into magnetic inputs through Equation (2).

Zoom In Zoom Out Reset image size

It is easy to see that under this simple one-dimensional scheme all strands will eventually acquire positive values of B_h. To keep the simplicity of the model and to be consistent with the fact that neighbor strand footpoints tend to move in different directions, we alternate the sign of the displacements d among neighbor strands. This means that, identifying each mesh site with an index i, we add the random values obtained from the distribution to strands with even i, and we subtract them from strands with odd i. Thus, in the long run even strands acquire positive B_h while odd strands acquire negative B_h. This alternation of δB_h simulates in one dimension the magnetic stress produced by footpoints moving in different random directions.

The efficiency of the process that injects "horizontal component" of the magnetic field through the motion of the footpoints depends on the average distance between neighboring strands. If the strands are closer than d₀, then the process can be considered efficient at tangling the strands and accumulating magnetic stress. However, if there is no injection of new flux, then the distance between strands will increase as the footpoints disperse. The driving process will become progressively less efficient. To simulate the decay of the driving mechanism due to footpoint dispersion, at a chosen point during the system evolution we make the distribution shown in Figure 3(b) slowly evolve to the symmetrical distribution shown in Figure 3(c). In the final state, the driver distribution is completely inefficient at creating new B_h. As we will see below, this evolution of the driver leads to the decay of the system. The start time and the duration of the evolution from the initial to the final distribution are set as parameters of the model.

In Section 3.1 we defined B_c = B_vtan θ_c. Therefore, tan θ_c is one of the relevant parameters of the model. At each time step, after adding the random δB_h values to all the strands, we compare B_h on each strand with its immediate neighbors. For a strand having index i we compute the difference ΔB_h(i) = B_h(i) − B_h(i + 1). We do this for all the strands in the array. We use periodic boundary conditions, so strand i = N_m is compared with strand i = 1, and in that case ΔB_h(N_m) = B_h(N_m) − B_h(1). Next, we identify all the indices where |ΔB_h(i)|>B_c. For all the sites where this happens, we redistribute the horizontal field according to Equation (5).

As we discussed in Section 3.1, each field redistribution is associated with a reconnection event that frees an energy ΔE, in erg cm⁻³, according to Equation (6). This energy is used to heat the plasma in both strands involved in the event. Before proceeding to the next time step, we repeat the previous test and corresponding field redistribution as many times as necessary until all |ΔB_h(i)| values are less than B_c. Once this state is reached, the system is allowed to evolve to the next time step, a new set of random δB_h values is added to the B_h of the strands, and the process described above is repeated.

During the evolution of the system we keep record of the energy dissipated in each strand at each time step. Here, we assume that the reconnection events (ΔE) occur on a timescale that is much shorter than the time step duration (≈10³ s). As we defined in Section 3.1, the sum of all the ΔE "energy packages" that heat the plasma in the strand during a given time step is called a nanoflare. The final output of the CA model is the array of nanoflare energies dissipated in each strand as a function of time.

In the above description we assumed that when the critical condition (|ΔB_h|>B_c) is fulfilled, all the available field ΔB_h is distributed (reconnected) between the two strands. This is not realistic. On the Sun, each strand of the tangled field is wrapped around many other strands. It is only the horizontal field component along the section of mutual contact that is free to reconnect and heat the plasma. We account for this in our simple model by introducing a parameter f. The initial response to the reconnection is to change the horizontal field locally by an amount ΔB_h/2. This change is subsequently spread out along the entire length of strand. The adjustment happens at the Alfvén speed and is very rapid. Assuming a minimum Alfvén speed of $1000\, \rm km\, s^{-1}$ for the solar corona gives a time of $100\, \rm s$ or less. The end result is that the horizontal field changes by fΔB_h/2 along the entire strand. It can be demonstrated that in this case the energy release on each of the two strands per unit volume is

$$\begin{equation} \Delta E = \left(\frac{1}{8\pi }\right)f\left(1- \frac{f}{2} \right) \frac{(\Delta B_{h})^{2}}{2}. \end{equation} \tag{ 7 }$$

Note that the parameter f affects the time required for a pair of adjacent strands to once again reach the critical condition through continued driving. If f is small, the strands will never be far from critical and nanoflares will recur frequently. If f is large, it will take many steps to once again reach the critical level of stress. Of course the nanoflares will be weaker in the former case than in the latter, so the temporally averaged heating will not be greatly different. There are small differences though, reflected in the f dependence of Equation (11), discussed later.

For the analysis, we vary the values of the model parameters B_v, L, tan θ_c, f, N_m, and τ. The parameter τ is the duration of the nanoflares as defined later in Section 3.3. Table 1 shows the different values used. For the vertical magnetic field we use values in the range 20–120 G, typical of AR loops (Mandrini et al. 2000). We vary the length of the loops also considering typical AR values (20–120 Mm). The tangent of the critical angle (tan θ_c) varies between 0.1 and 1.0, corresponding approximately to angles: 6°, 11°, 22°, 31°, 39°, and 45°. The reconnection factor f varies between 0.1 and 1.0, meaning that 10%–100% of the available magnetic stress is released by reconnection. The total number of strands, N_m, goes from 20 to 1000, which is roughly the number of elemental magnetic flux tubes believed to be present in a coronal loop. We vary one parameter at a time keeping the rest of them at the intermediate values shown in bold face in Table 1.

Table 1. Numerical Values Used for the Parameters of the Model

Bv (G)	20	40	60	80	100	120
L (Mm)	20	40	60	80	100	120
tan θc	0.1	0.2	0.4	0.6	0.8	1.0
f	0.1	0.2	0.4	0.6	0.8	1.0
Nm	20	50	100	200	500	800	1000
τ	100	500	1000	2000	3500	5000

Download table as:  ASCIITypeset image

As we stated above, the output of the CA model is an array of nanoflare energies dissipated in each strand as a function of time, with time steps of ≈10³ s. Figure 4 helps to illustrate the typical evolution of the system. For the example we use the parameter combination shown in bold face in Table 1 and a total duration of 72 time steps (${\approx} 20\, \rm hr$). The upper panel shows the evolution of the mean energy of the system, defined as the sum of the magnetic energies of all the strands (B²_i/8π) divided by the number of strands. The middle panel shows the evolution of the mean nanoflare energy, which is the sum of the energies of all nanoflares occurring at each time step divided by the number of strands. The evolution is such that the number of reconnection events increases from zero (rise phase) to a quasi-stationary level (main phase). In the example of Figure 4 this transition occurs around time step number 20. The system reaches this state abruptly, not asymptotically, and maintains it for as long as the driver follows the distribution of Figure 3(b). Note that once the system reaches the main phase, the magnetic energy (Figure 4, upper panel) maintains an approximately constant level, except for fluctuations due to the nanoflares. As soon as the driver begins to evolve toward the distribution of Figure 3(c), the number of events decays until they almost disappear at the end of the system's evolution (Figure 4, middle panel). Likewise, once the decay phase begins, the fluctuations of the system energy begin to decrease rapidly and stop. Summarizing, the evolution of the full system qualitatively follows the evolution of the observed coronal loops in terms of rise, main, and decay phases.

Zoom In Zoom Out Reset image size

3.3. Plasma Response and Light Curve Construction

Since our main motivation is to compare the output of the CA model with the loops studied in López Fuentes et al. (2007), the next step is to transform this output into an observed signal such as the light curve of a coronal loop. To do this, we simulate the plasma response using the EBTEL model developed by Klimchuk et al. (2008). The temperature and density of the coronal plasma tend to be reasonably uniform along the magnetic field, and EBTEL computes the evolution of the spatially averaged values. This approximate zero-dimensional solution is quite similar to the exact solution given by one-dimensional hydrodynamic codes, but it takes 3–4 orders of magnitude less time to compute. EBTEL is therefore ideally suited for our study in which we simulate the evolution of thousands of individual strands.

The main input of EBTEL is the heating rate profile, i.e., the variation of the heating rate with time in the strand. Here, we assume that nanoflares have a triangular profile of a given duration τ. The upper panel of Figure 5 shows a series of "triangular" nanoflares that heat one of the strands of the CA model during a portion of the main phase. The time integral of each triangular heat function corresponds to the total energy provided by the particular nanoflare. The duration of the nanoflares shown in Figure 5 is 1000 s. This means that, in this case, the available heat is distributed (with a triangular profile) along a full time step of the CA model. The other two panels of Figure 5 show the output of EBTEL: the plasma temperature and density. During the time interval covered in the panels there are both isolated events (like the nanoflare starting at t = 4.3 × 10⁴ s, gray arrow) and "trains" of events occurring in a rapid succession (such as the one starting at t = 4.7 × 10⁴ s, black arrow). Note that while the plasma relaxes almost fully after the isolated event, during the succession of multiple events the plasma is still in a state of high density and temperature due to the previous nanoflare when the next event begins. As we will see in the following sections, this has consequences in the way the plasma emits and the observed loop intensity evolves.

Zoom In Zoom Out Reset image size

The loops studied in López Fuentes et al. (2007) have measured densities in the range 0.9–1.9 × 10⁹ cm⁻³ and temperatures between 1.2 and 2.1 MK. These values are reproduced by parameter values in the ranges given in Table 1.

Once we have the temperature and density evolution of each strand, we compute the observed intensity using the emission measure (EM) and the known temperature-dependent sensitivity of the instrument. The EM per pixel of a single unresolved strand is EM = n²ΔV, where ΔV is the volume of the strand portion covered by a pixel. This is roughly ΔV = l_pix A_s, where l_pix is the pixel dimension and A_s is the cross-sectional area of the strand. The contribution to the intensity in a pixel is then I = EMS(T), where S(T) is the SXI response function. We sum the contributions of all the strands to determine the intensity of the loop and, finally, we multiply this intensity by a factor that accounts for the proportion of the loop covered by a single pixel. The final output is the observed intensity of the loop per pixel per second as a function of time. The data are sampled using the instrument temporal resolution to simulate observed light curves. In Figure 4, lower panel, we show a light curve obtained in this way for the model example of the upper and middle panels. Comparing the lower and middle panels it can be noticed that the light curve follows (though smoothed) the general fluctuations produced by the summed contribution of the nanoflares.

The lower panels of Figure 1 show two model light curves obtained in the way described in Section 3.3. The cases shown have the same properties as the observed light curves of the upper panels. Given the random nature of the nanoflares, we do not expect an exact fit of the model to the observed data. For comparison, we instead consider the light-curve properties studied in López Fuentes et al. (2007); see Table 2.

The values of the model parameters used to construct the synthetic light curves of Figure 1 are, for Loop 1: B_v = 37 G, tan θ_c = 0.15, and f = 0.2; and for Loop 4: B_v = 20 G, tan θ_c = 0.2, and f = 0.1. For both cases we use N_m = 100 and τ = 1000 s. For the loop length and diameter we use the measures obtained in López Fuentes et al. (2007): L = 110 Mm for loop 1, L = 70 Mm for loop 4, and a diameter of 1.8 × 10⁹ cm for both.

We compute for both observed and synthetic light curves the mean intensity of the main phase (I_m) and the duration and timescales of the rise and decay phases. To obtain these values we first identify, by visual inspection, the approximate start and end times of each phase. Using the same definitions as in López Fuentes et al. (2007), we have for the rise phase timescale

$$\begin{equation} T_{\rm rise} = \frac{\bar{I}}{(\Delta I/\Delta t_{\rm rise})}. \end{equation} \tag{ 8 }$$

Here, $\bar{I}$ is the mean intensity during the rise phase, ΔI is the intensity variation, and Δt_rise is the duration. A similar definition is used for the decay phase. As a proxy for the amplitude of the intensity fluctuations, we compute the rms of the intensity variation during the main phase relative to its 10-point running average. In Table 2, we compare these properties for observed and modeled light curves. It can be seen that the model reproduces the light-curve properties consistently.

We remind the reader that the observed light curves shown in Figure 1 have been further processed (residual background removed, see Section 2) from the versions studied in López Fuentes et al. (2007). Therefore, the properties listed in Table 2 do not necessarily coincide with the corresponding values presented in the previous paper.

Regarding the durations and timescales it is interesting to see that, since the light curves have initial intensities very close to zero, $\bar{I}$ in Equation (8) is approximately ΔI/2. Therefore,

$$\begin{equation} T_{\rm rise} \approx \frac{\Delta t_{\rm rise}}{2}. \end{equation} \tag{ 9 }$$

We then expect the timescales to be approximately one half of the durations (compare the corresponding values in Table 2 for both the rise and decay phases). The relation from Equation (9) is not valid, in general, for the light curves studied in López Fuentes et al. (2007) in which there is still a remnant background component of the intensity.

To investigate how the model parameters determine the properties of the loops, we vary each of them separately according to Table 1 while keeping the rest fixed at the values shown in bold face (see Section 3). In Figure 6, we show log–log plots of the mean energy per unit volume released in reconnection events, ΔE (see Equation (7)), versus different model parameters. The lines are least-squares fits of the data, and their slopes are indicated in the panels. According to these plots, in our CA model the energy of the reconnection events approximately scales with B^2.01_v and (tan θ_c)^1.89. This can easily be explained by noting that in Equation (7), when reconnection occurs, ΔB_h ≈ B_c = B_vtan θ_c. Therefore, it is expected that ΔE ∝ B²_v(tan θ_c)².

Zoom In Zoom Out Reset image size

The horizontal field decreases by an amount fΔB_h/2 during a reconnection event. Several time steps are typically required to rebuild the field to the critical level whereupon another event occurs. From Equation (2) we see that the number of recovery time steps is given by

$$\begin{equation} N_{r} \approx \frac{f}{2} \frac{L}{d} \tan \theta _c. \end{equation} \tag{ 10 }$$

From this we obtain the mean nanoflare heating rate (averaged both in time and along the strand)

$$\begin{equation} \left< Q \right> \approx \frac{\Delta E}{N_r t_{s}} = \left(1 - \frac{f}{2} \right) \frac{v B_{v}^{2} \tan \theta _{c}}{L}, \end{equation} \tag{ 11 }$$

where t_s is the time step duration and v = d/t_s is the photospheric driving velocity. We have assumed that there is one reconnection event per nanoflare, whereas often there are two, but this is not important for identifying scaling relationships. In Figure 7, we show log–log plots of the mean heating rate versus B_v, tan θ_c, and L. The dependence on the parameters is B^2.02_v, L^−1.1, and (tan θ_c)^0.97, which is clearly explained by Equation (11). It is easy to see that <Q> does not depend on the number of strands, N_m, since both the driver and the critical condition are applied to the strands individually.

Zoom In Zoom Out Reset image size

We can also analyze scaling relations for the properties of the light curves obtained with the model. For each combination of the parameters from Table 1, we obtain synthetic light curves following the procedure described in Section 3.3. In Figure 8, we present log–log plots of the mean intensity of the main phase, I_m, as a function of the same parameters as in Figure 7. The scalings found in this case are B^2.37_v, L^−0.8, and (tan θ_c)^1.15. Clearly, I_m is related to the nanoflare energy, but the dependence is complex and is affected by both the detailed response of the plasma and the nonuniform temperature sensitivity of the instrument. I_m is not linearly proportional to the total radiation emitted by the strand. To understand the origin of the dependence, let us now consider the case of quasi-static equilibrium. This is a reasonable representation if the time interval between nanoflares is short compared to a cooling time. When such conditions apply, the conductive and the radiation terms of the energy balance equation are comparable in magnitude to each other and to the mean heating rate, <Q> (Vesecky et al. 1979):

$$\begin{equation} \left< Q \right> \approx \frac{2}{7} \kappa _{0} \frac{T^{7/2}}{L^{2}} \approx n^{2} \Lambda _0 T^{b}. \end{equation} \tag{ 12 }$$

Here, κ₀ is the coefficient of thermal conductivity, n is the electron number density, and Λ₀ and b are constants. From the first part of Equation (12) we can obtain for the temperature:

$$\begin{equation} T \propto \left< Q \right>^{2/7}L^{4/7} \propto B_{v}^{4/7} L^{2/7} (\tan \theta _c)^{2/7}. \end{equation} \tag{ 13 }$$

For the temperature range of the SXI loops studied in López Fuentes et al. (2007), a reasonable single-power value for the radiation loss is b ≈ −0.5 (RTV, 1978). Using this, and Equations (12) and (13) we have

$$\begin{equation} n \propto \frac{T^{2}}{L} \propto B_{v}^{8/7} L^{-3/7} (\tan \theta _c)^{4/7}. \end{equation} \tag{ 14 }$$

Zoom In Zoom Out Reset image size

The observed loop intensity is I = n²S(T), where S(T) is the response function of the instrument, which we can express as S(T) = S₀ T^a. We found that a ≈ 0.7 for the temperature range of the loops studied in López Fuentes et al. (2007). Replacing Equations (13) and (14) in the expression for the loop intensity, we finally obtain

$$\begin{equation} I_m \propto B_{v}^{2.69}\, L^{-0.66} (\tan \theta _c)^{1.34}. \end{equation} \tag{ 15 }$$

This expression is consistent with the scalings in Figure 8. Differences are to be expected since the strands are sometimes far from quasi-static equilibrium (Figure 5).

It might be expected that the intensity of the main phase depends on N_m, but this is not the case. The reason is as follows. The intensity contribution of each strand is proportional to the volume of the strand subtended by the instrument pixel, which is proportional to the strand cross section. Since we compute the strand cross section as the observed cross section of the loop (1.8 × 10⁹ cm) divided by N_m, an increase of N_m means a proportional decrease of each strand intensity, so the summed intensity of all strands remains unchanged.

As we mentioned above, another important property of the light curves is the amplitude of the intensity fluctuations, which we quantify from the rms of the intensity variation during the main phase. These fluctuations are related to the frequency and relative sizes of the nanoflares. As we briefly discussed in Section 3 the factor f, which relates to the fraction of the strand length that is involved in the reconnection, determines the waiting time for two neighbor strands to recover the critical condition. If f is small, the waiting time will be short and the nanoflares will tend to be smaller and more frequent. On the other hand, as f tends to 1 the nanoflares will be larger and less frequent. Therefore, it is expected that the intensity fluctuations increase with f. The upper panel of Figure 9 shows precisely that.

Zoom In Zoom Out Reset image size

The intensity fluctuations are also affected by the number of strands, N_m. Since the intensity of the loop is the summed contribution of all the strands, given the random nature of the release events, the larger the number of strands, the smoother the signal. This is shown in Figure 9, middle panel. We see that the rms intensity fluctuation decreases as the number of strands increases. The effect levels off, once a sufficiently large number of strands are present.

We also explore how the nanoflare duration, τ, affects the intensity fluctuations. The different values used are presented in Table 1, last row. The results are shown in Figure 9, lower panel, where it can be seen that the rms fluctuation also decreases as τ increases. The reason is that individual strands evolve more slowly as the nanoflares get longer. For long enough nanoflares the strands behave in a quasi-static manner. As is the case with N_m, the effect levels off once the nanoflares are sufficiently long.

Finally, although the dependence is much weaker than for the previous parameters, the rms variation tends to increase with all the parameter changes that increase the intensity, i.e., the increase of B_v and tan θ_c, and the decrease of L. This can be understood in terms of the time it takes the plasma to cool fully after a nanoflare. The cooling time varies inversely with the nanoflare energy and directly with the loop length (Cargill & Klimchuk 2004). Shorter cooling times result in more bursty emission and greater fluctuation amplitudes. The increase of the fluctuation amplitude with the increasing intensity is consistent with the observational results obtained recently by Sakamoto et al. (2008, and references therein) and discussed also in Vekstein (2009). They also found that the amplitude of the intensity fluctuations of Yohkoh/SXT loops is approximately 8% of the value of the loop intensity. We find similar values in our own observations and the results obtained with the CA model (see the main phase rms variation in Table 2).

Another important property of the light curves is the rise duration, Δt_rise, which is the time it takes to reach the main phase after the intensity first starts to rise. Note that the intensity does not increase immediately when footpoint driving begins (see Figure 4, middle and lower panels), because it takes time for the strands to reach the critical condition. The first strands to reach it are those that by chance have larger-than-average step sizes. We define the waiting time T_w to be the elapsed time between the start of the evolution (t = 0) and the main phase beginning. It can easily be seen that T_w is simply the number of steps required to reach the critical condition, |B_h| = B_c/2, where the factor of 2 is because odd and even strands are tilted in opposite directions. The average |δB_h| per time step that is injected by the driving is only half the value given by Equation (2), because one in four steps causes the stresses to decrease (Figure 3). Noting that B_c = B_vtan θ_c, we find that the waiting time is given by

$$\begin{equation} T_{w} \approx \frac{B_{c}/2}{\delta B_{h}/2} \approx \frac{L \tan \theta _{c}}{d} \end{equation} \tag{ 16 }$$

in units of the time step duration (1000 s). In Figure 10 we show log–log plots of T_w versus the parameters. The scalings are L^0.88 and (tan θ_c)^0.85. This is generally consistent with Equation (16), and we suggest that small deviations are due to the subjective nature of identifying the start of the main phase. Note that T_w is expected to be independent of B_v, which we have confirmed to be true from our simulations. We emphasize again that the waiting time is different from the rise duration. The rise duration is more closely related to the spread of the driver distribution (e.g., the Gaussian width in Figure 3(b)) than it is to the time required for the average strand to reach critical.

Zoom In Zoom Out Reset image size

The duration of the decay phase depends on both the form of the final driver distribution (Figure 3(b)) and the timescale for transitioning from the initial to final distributions. We leave the details for a later study.

We develop a pseudo-one-dimensional model based on a CA that reproduces the main properties of the evolution of coronal loops observed in soft X-rays. As we describe in Section 3, the model is based on the idea that coronal heating is due to the sudden reconnection of tangled magnetic strands (Parker 1988). Using typical solar ranges for the model parameters we qualitatively reproduce the observed intensity amplitude, fluctuations, and evolution timescales of soft X-ray loops studied in López Fuentes et al. (2007). The synthetic loops obtained with the CA+EBTEL model combination also reproduce the measured physical properties of the observed loop plasma such as density and temperature.

Although our model is related to usual models of self-organized criticality (SOC; see e.g., Charbonneau et al.2001), it differs in some important aspects. In most SOC models, the rules for the critical condition and the redistribution of the field among neighboring strands are based on mean values of the field strength. The strength of each strand is compared with the mean value of its neighbors, and if the critical condition is fulfilled the field is redistributed equally among those neighbors. One may classify this type of model as isotropic. As we describe in Section 3, our CA model is anisotropic in this respect. The other important difference is that SOC models use null boundary conditions. The mesh points (strands) at the boundary are set to have zero magnetic field (what we call B_h in Section 3) at all times. Here, we use instead periodic boundary conditions. It may be partly due to these differences that the nanoflares produced with our model do not follow power-law distributions, which are characteristic of SOC models (for an anisotropic SOC model associated with power-law distributions, see e.g., Morales & Charbonneau 2008). We must stress that our main goal here is not to obtain such distributions, but instead to explain the evolution of loops in the frame of the nanoflare model. It is worth pointing out, though, that our CA model reaches a steady critical state in much the same way as SOC models do (see Figure 4).

As we mention in Section 1, in recent years there has been a debate over the issue of whether coronal loops are isothermal or multithermal. The argument goes that if loops are isothermal they are likely to be monolithic structures, while multithermality implies that loops are composed of different unresolved strands evolving independently. There seems to be evidence in favor of both kinds of observations (see, e.g., Schmelz et al.2009), so a consensus is forming that both scenarios actually occur. Let us analyze this possibility in terms of the model presented here. For instance, if the reconnection mechanism is not very efficient (e.g., a small f factor, see Section 3) and all nanoflares are approximately the same size, then the heating injection will be difficult to distinguish from a single source of constant heating. In that case, the strand will keep more or less the same temperature along the main phase of the evolution. If all the strands that comprise the loop have similar physical conditions (it is expected that they all have approximately the same field strength, for example), then the loop will be nearly isothermal. These loops can easily be taken as examples of quasi-static evolution. On the other hand, if nanoflares occur less frequently, so the waiting time between consecutive nanoflares is not negligible, the strands have time to cool to lower temperatures between one nanoflare and the next. Since nanoflares occur randomly, at any given time during the evolution of the loop there will be strands simultaneously having different temperatures. That would correspond to a multithermal loop observation.

There is yet a third possibility. If the nanoflares occur only once in each strand, without repeating, and if they all occur at approximately the same time, then the loop will appear roughly isothermal at any given moment. The temperature will decrease on a cooling timescale, however, and the loop will be short lived. We call this scenario a short-duration nanoflare "storm" (Klimchuk 2009). Many warm EUV loops can be explained in this way. Hot soft X-ray loops tend to live much longer than a cooling time, on the other hand (López Fuentes et al. 2007). If they are indeed isothermal, then they can only be explained by rapidly repeating nanoflares or by truly steady heating.

The three alternative scenarios described above relate to the question weather the same loop is observable by instruments that are sensitive to different temperatures. In the case of nanoflares that repeat rapidly in each strand of the loop bundle, the loop will be visible only over a narrow range of temperatures. It will be seen either in hot soft X-ray emission or in warm EUV emission, but not in both. In the case of nanoflares that repeat with a time delay longer than the cooling time, the loop will be simultaneously visible in both emissions. In the final case of a short-duration nanoflare storm, the loop will be visible first in soft X-rays and then in the EUV.

Observational evidence has been presented for all three possibilities. The question of which is most prevalent is not yet answered. Loops produced by short-duration storms have been clearly identified and can be found at many locations within ARs (e.g., Ugarte-Urra et al. 2006, 2009). On the other hand, soft X-ray loops are most pronounced in the central cores of ARs, where EUV loops are less common (Schmieder et al. 2004). There are also reports of soft X-ray and EUV loops that are almost but not exactly co-spatial (Nitta 2000; Nagata et al. 2003). We attempted to study the warm emission associated with the SXI loops studied in López Fuentes et al. (2007). Unfortunately, TRACE data are not available, so we had to rely on lower cadence observations from SOHO/EIT. Since only four images per day were taken in each channel, we could only compare the hot and warm emission at a few times during the loop evolution. We found that some of the SXI loops had an EIT counterpart but others did not. Clearly, more systematic observations of loops in different wavelengths are needed.

It is possible that there are two well-differentiated populations of AR loops. Hotter and longer-lived loops are preferentially located at the core of ARs (Antiochos et al. 2003, Ugarte-Urra et al. 2009), while the periphery is dominated by longer, shorter-lived, cooler loops (Del Zanna & Mason 2003). If this is indeed the case in general, then the different locations and evolutions of these loop classes may be related to different physical properties at their origin. The model we present in this paper applies specifically to hot loops. Whether it can explain warm loops at the periphery of ARs is unclear. We can nonetheless speculate on the difference in these two regions. The strength of the coronal magnetic field is generally greater in the core of an AR and this will tend to produce stronger heating. Our model also depends fundamentally on the properties of footpoint shuffling. It may be that these properties vary systematically across the AR. This is an important question that we hope will be answered with the latest magnetogram observations from SOT/Hinode. These observations will hopefully also confirm our proposal that the slow decay phase of hot loops is associated with the dispersal of magnetic footpoints to the point where magnetic field tangling is no longer efficient.

In Section 4, we found a series of scaling laws relating observed properties of the loop evolution to different physical parameters of the model. Some of these parameters, such as the loop length, can actually be obtained directly from the observations. The strength and tilt of the magnetic field at the coronal base can be inferred by combining photospheric magnetogram observations with models of the coronal structure (see, e.g., López Fuentes et al. 2006; Mandrini et al.2000). Therefore, these scaling laws are predictions of the model that can be tested observationally.

From the discussion in the previous paragraphs we conclude that among future important investigations are the systematic study of loop observations in different wavelengths and the analysis of possible scalings of the evolution parameters (timescales, mean intensities, fluctuations) with the loop physical properties such as length and magnetic field. From the theoretical side, the model presented here is still very simple, and more sophisticated versions should be developed, such as two-dimensional models that better account for the vector nature of the magnetic field and the fact that each magnetic strand is wrapped around multiple other strands. We have already started work in this direction.

The authors thank Paul Charbonneau for enriching discussions and the anonymous referee for useful comments and suggestions. This work was supported by the NASA Supporting Research and Technology and Living With a Star Programs.