A&A 420, 719-728 (2004)
DOI: 10.1051/0004-6361:20034570

On the distribution of magnetic energy storage in solar active regions

T. Fragos - E. Rantsiou - L. Vlahos

Department of Physics, Aristoteleion University of Thessaloniki, 541 24 Thessaloniki, Greece

Received 24 October 2003 / Accepted 4 February 2004

Abstract
A two-dimensional probabilistic Cellular Automaton is used to model the appearance of active regions at the solar surface. We assume that two main competing processes control the magnetic field evolution at the solar surface (1) the magnetic field is locally enhanced by the flux emergence and/or the coalescence of emerged magnetic flux and (2) it is diminished by flux cancellation or diffusion. The flux emergence follows a basic percolation rule; it is more probable at the points were magnetic flux already exists. The magnetic field is also enhanced when magnetic fields of the same polarity collide. The flux cancellation is due either to the gradual diffusion of the magnetic field, when it is isolated, or to the partial release of energy when opposite magnetic field lines collide. The percolation model proposed in this article is capable of reproducing the statistical properties of the evolving active regions. The evolving simulated magnetograms, derived from our model, are used to estimate the 3-D magnetic fields above the photosphere using constant $\alpha$ force-free extrapolation techniques. Based on the above analysis we are able to estimate a variety of observed statistical characteristics, e.g. the size and flux distribution of the magnetic fields at the solar surface, the fractal dimension of the magnetic structures formed at the photosphere, the energy release frequency distribution, the waiting time distribution of the sporadic energy releases and the statistical properties of the steep horizontal magnetic field gradients in the extrapolated coronal magnetic field. Our main conclusion is that the photospheric driver plays a crucial role in the observed flare statistics, and the solar magnetograms, when interpreted properly, carry important statistical information for the solar coronal activity (coronal heating, flares, CME etc.).

Key words: Sun: activity - Sun: magnetic fields - Sun: flares - Sun: photosphere - plasmas

1 Introduction

The physics of the solar corona contrasts with the physics of the convection zone (CZ), where the plasma pressure dominates the magnetic field. The magnetic field is accurately measured only at the photospheric level. What we observe are the remnants of the sub-photospheric activity and at the same time constitutes the time-dependent boundary for the chromosphere and coronal physics.

The most active phenomena above the solar surface are related to active regions. It is for these reasons that their study has attracted the attention of many observers and theorists in recent years. The main focus of current theoretical studies is on the question: How is the subphotospheric activity mapped onto the formation and subsequent evolution of active regions. Two interrelated questions are usually posed: (1) where are the magnetic fields formed; and (2) how are they manipulated by the convection zone?

The formation and evolution of large-scale magnetic flux tubes inside the convection zone is an important theoretical problem, which still remains open. It is believed today that flux tubes are formed at the base of the convection zone and rise to the surface through buoyant forces. During their buoyant rise, flux tubes are influenced by several physical effects, such as the Coriolis force, magnetic tension, drag and large-scale convective motion. The large-scale magnetic flux tubes undergo a dramatic evolution before reaching the solar surface and probably split, due to hydrodynamic forces, to form small-scale fibrils. The lack of understanding of the relative importance of these forces in the formation and evolution of active regions has lead many observers to develop large-scale statistical studies to undrestand the characteristics of active regions (Howard 1996).

Observations support the idea that active regions on the sun are formed by the emergence of a large number of separate small intense flux bundles of the order of 500 G and radii $\sim$200 km ($\sim$1018 Maxwells) when they first appear, but that soon are compressed to $4{-}2\times 10^3$ G over 100 km in the photosphere (Brants & Zwaan 1982; Brants 1985; Brants & Steenbeck 1985). The flux bundles show a strong tendency to cluster together to form pores and sunspots as long as fresh flux continues to emerge. This effect vanishes when emergence ceases (Zwaan 1985).

The separate magnetic flux tubes observed in the photosphere expand in the chromosphere to fill the available space. Hence the clustering observed at the photosphere compresses the close-packed field at higher levels and therefore is opposed by the magnetic stresses. It is obvious that the clustering does not appear spontaneous and it must be driven by hydrodynamic forces beneath the photosphere (Parker 1979, 1992).

Many models have been proposed for the formation and evolution of active regions, such as the rise of a kink-unstable magnetic flux tube (see Moreno-Insertis 1992) or the statistical description of the dynamical evolution of large-scale, two-dimensional, fibril magnetic fields (Bogdan & Lerche 1985). Schrijver et al. (1997a), using the methods originally invented by Bogdan & Lerche (1985) for the convection zone, constructed a set of rate equations, that are valid only at regions of weak gradients (quiet sun). Schrijver et al. (1997a) searched for a balance between the flux emergence and the competing processes of diffusion (partial cancellation, coalescence and fragmentation) Schrijver et al. (1997b) extended the magnetic carpet model to active regions. In this model, the collision frequency $\nu $ between the fibrils is assumed to vary quadratically with the number density of their concentration (Nt) i.e. $\nu=\ell N_t^2$. Their fragmentation rate K is assumed to be proportional to the local magnetic flux $K=k\Phi.$ The solutions deduced from the rate equations are very sensitive to the spectrum of the emerging flux concentrations that replace the cancelling flux. The solutions are also controlled by the ratio of the coefficients k and $\ell$which are both loosely related to the physical processes at the photosphere. Longcope & Kankelborg (1999) use the statistical results of Schrijver et al. (1997a) to study the interaction of randomly moving photospheric magnetic flux elements of opposite signs and the appearance of an X-ray bright point.

Independent models have also been developed using the anomalous diffusion of magnetic flux in the solar photosphere to explain the fractal geometry of the active regions (see Lawrence 1991; Lawrence & Schrijver 1993; Milovanov & Zeleny 1993) and a simple percolation model was also developed where clusters were formed by randomly placed fibrils at the photosphere (see Schrijver et al. 1992).

A new percolation model based on well-known observations was developed to simulate the formation and evolution of active regions by Wentzel & Seiden (1992) and Seiden & Wentzel (1996). In this model the flux emergence was based solely on the observational fact that "magnetic flux emerges where there is flux already''. Fibrils follow a random walk at the surface and collide with other fibrils, merging when they have the same polarity or cancelling when they meet fibrils of the opposite polarity (Wang et al. 1989; DeVore et al. 1985). The percolation model proposed by Seiden & Wentzel (1996) incorporated all the diffusion characteristics present in the magnetic carpet model of Schrijver et al. (1997a) but differs in one fundamental aspect: the process with which the magnetic flux emerges at the surface. The magnetic carpet introduces new flux when it is cancelled while the percolation model uses the rule described above. Both reach a steady state at a given activity level. The percolation model explains the observed size distribution of active regions and their fractal characteristics (Meunier 1999). Vlahos et al. (2002) analysed the sporadic energy release through flux cancellation (reconnection) when flux tubes of opposite polarities collide and analyzed the statistical properties of the energy released. They showed that the percolation model can easily interpret the statistical characteristics of the Ellerman bombs, coronal bright points, Ha bright points and transition region impulsive EUV emission.

In this article, we expand the percolation model proposed initially by Seiden & Wentzel (1996) and recently further developed by Vlahos et al. (2002) to estimate the dynamic evolution and the statistical properties of the magnetograms formed by our simulations. Using standard force-free extrapolation techniques, we estimate the evolving 3-D structures above the photospheres and we investigate in detail the statistical properties of the horizontal discontinuities.

2 Active region formation and evolution

2.1 Observational constraints

A variety of well-established observational constraints are important for any realistic model which will claim that it can explain the formation of active regions. We list a few of them here (see more details in Zirin 1988; Howard 1996, and the references cited therein):
1.
magnetic flux emerges and remains in the photosphere in the form of discrete structures called fibrils;
2.
magnetic flux emerges where there is flux already present;
3.
active (magnetized) and empty (non-magnetized) regions persist for relatively long times;
4.
active regions cover only a small fraction of the solar surface (around $\sim$$ 8\%$);
5.
fibrils undergo a random walk in the photoshere (Muller et al. 1994);
6.
the size and magnetic flux distribution exhibit two well-defined regimes, (1) a power law in the range $40{-}400~\rm Mm^2$ and an exponential rollover for the larger areas (Harvey 1993; Harvey & Zwaan 1993; Meunier 2003);
7.
the fractal dimension, $D_{\rm F}$, of the active region structures observed in the photosphere has been found, using high-resolution magnetograms, to be in the range $1.3<D_{\rm F}<1.8$ (Balke et al. 1993; Meunier 1999);
8.
properties of the magnetic field in active regions, such us the total magnetic flux of each active region and its maximum magnetic field, have frequency distributions that exhibit a well-defined power law at small values (Meunier 2003);
9.
the distribution function of the energies released in the Ellerman bombs exhibits a power-law shape with index $\sim$-2.1 (Georgoulis et al. 2002). Different attempts have been made to calculate the energy of the Ellerman bombs. The power-law index, however, is model dependent and it can change with the assumptions, as stressed by Georgoulis et al. (2002);
10.
the observed distribution of waiting times $\Delta t$ between X-ray solar flares of greater than C1 class listed in the Geostationary Operational Environmental Satellite (GOES) catalog exhibits a power-law tail $\sim$ $(\Delta t)^\gamma$ for large waiting times. The power-law index $\gamma$varies with the solar cycle; for the minimum phase of the cycle the index is $\gamma=1.4\pm 0.1$ but for the maximum phase of the cycle the index is $3.2\pm 0.2$ (see for example Wheatland & Litvinenko 2002).
The main question posed in this article is: Can we construct 3-D models that can simultaneously explain all the abovementioned observations?

2.2 The percolation model

The basic ingredient of the model presented here for the rise of a magnetic field to the surface and the formation of active regions is a well-established non-linear theory, the percolation theory. The percolation theory is well explored physically and mathematically (Stauffer & Aharony 1985; Deutscher et al. 1983). It has been used for the study of galactic spirals using stochastic self-propagating star formation (Seiden & Gerola 1982; Seiden 1983; Schulman & Seiden 1986) and the evolution of magnetic fields in the accretion disks around compact objects (Pavlidou et al. 2001).

We propose in this article that the main physical properties, as derived from the observations of the evolving active regions, can be summarized in simple Cellular Automata (CA) rules:

A 2-D quadratic grid with $256 \times 1280$ cells (grid sites) is constructed, in which each cell has four nearest neighbors. The grid is assumed to represents a large fraction of the solar surface, more precisely, as the the left and right boundaries are periodical, it represent a zone around the equator of the sun of 72 degrees width. Initially, a small, randomly chosen percentage ($1\%$) of the cells is magnetized (loaded with flux) in the form of positively (+1) and negatively (-1) magnetized pairs (dipoles); the rest of the grid points are set to zero. Positive and negative cells evolve independently after their formation, but their percentage remains statistically equal.

The dynamical evolution of the model is controlled by the following probabilities:

P: the probability that a magnetized cell is stimulating the appearance of new flux at one of its nearest neighbors. Each magnetized cell can stimulate (add one positive or negative flux unit) its neighbors only the first time step of its life. This procedure simulates the stimulated emergence of flux which occurs due to the observed tendency of magnetic flux to emerge in regions of the solar surface in which magnetic flux had previously emerged. The physical interpretation of this rule of the CA is the following: magnetic twists, or even knots, travel along larger-scale field lines, as shown by Parker (1979) for force-free fields. But usually twists and knots cannot travel freely because they caught up in the field structure due to line tying. Only occasionally, especially when the field is already somewhat simpler than normal, the twists may travel until they encounter other twists. Then both annihilate and the field is further simplified. The release of flux tubes creates sufficient newly vacated space that nearby fields can relax while expanding into this space. The relaxation provides new opportunities for the twists and knots on those field lines to travel. They have a renewed chance to reach sites where dissipation is possible. Clearly there is then some probability that these fields, newly relaxed and simplified, become sufficiently simple to constitute a flux tube that also rises to the surface (see Wentzel & Seiden 1992 for a qualitative discussion of the relation of the probability P with the physical processes expected below the photosphere.).

D$_{\bf m}$: the flux of each magnetized cell has a probability $D_{\rm m}$ to move to a random neighboring cell, simulating motions forced by the turbulent dynamics of the underlying convection zone (see Simon et al. 1995). If the moving flux meets oppositely polarized flux in a neighboring cell, the fluxes cancel (through reconnection), giving rise to a "sporadic energy release''. If equal polarities meet in a motion event, the fluxes simply add up.

D: the probability that a magnetized cell is turned into a non-magnetized on one time-step if it is next to a non-magnetized cell. Each magnetized cell has probability 1-(1-D)n to become non-magnetized, where n is the number of its non-magnetized neighbors. This rule simulates two effects, the direct submersion of magnetic flux and the disappearance of flux below observational limits due to dilution caused by diffusion into the empty neighborhood. It gives us also a control over the lifetime of the flux-tubes on the photosphere.

E: the probability that a non-magnetized cell is turned into a magnetized one spontaneously, independently of its neighbors, simulating the observed spontaneous emergence of new flux. Parker (1992) suggested that flux tubes may start to rise in response to thermal plumes in the convection zone. Plumes not only punch through the magnetic field but carry some field with them which then continues to rise to the surface. If plumes were the only case of rising flux tubes, then the resulting surface structures would not exhibit any of the correlations that are apparent in the active regions. Thus in our model the flux emergence caused by plumes corresponds to the spontaneous emergence controlled by the probability  $\mathcal{E}$. Magnetic flux appears always in the form of dipoles. Whenever a new flux tube appears on the grid by stimulation or spontaneous emergence, an opposite flux tube also appears, to form a dipole. The distance between the two poles is $l+{\rm d}l$, where l is a constant; in the results presented here we use l=10 and ${\rm d}l$takes a random value between 0-5. The emergence of the oppositely polarized flux takes place in such a way that all active regions have the same orientation. During each timestep, we apply the four rules successively scanning the grid four times and we store the changes of each of the four scans in a dummy array. At the end of the timestep, we record the sporadic energy release and update our main array.

Our model can incorporate differential rotation of the grid. We decided not to include it, since after several tests we came to the conclusion that it has no effect on our statistical results. Seiden & Wentzel (1996) also indicated that differential rotation is unimportant for the young active regions studied.

2.3 Results

2.3.1 Size distribution and fractal dimension
The parameters used for the results reported here are P=0.179, D=0.005, $D_{\rm m}=0.05$ and E=10-6. Our model does not critically depend on E, as long as this probability is very small but not zero. When E=0 the model gradually dies out. The parameters $(P,D,D_{\rm m})$ are chosen such that, when following the evolution of our model and recording the magnetized cells, the percentage of the active cells is stabilized to a value close to the observed one (around $ 8\%$) (see Fig. 1). The time was calibrated by comparing the evolution of the clusters of active points in our model with that obtained using the observed active regions (see Sect. 2.3.1).

The level of the active cells where the stabilization is reached depends critically on the parameters  $(P,D_{\rm m},D)$. The observed level of activity at the photosphere places a constraint on their values. In the rest of the article these values will be kept fixed and will be called "the standard values''. These values are in no way unique; slight changes in any of the parameters do not dramatically affect the behavior of the model as long as we readjust the other two free parameters. The level of activity will increase dramatically if the critical percolation coefficient is reached. In Fig. 1 we present the evolution of our model using three different values of P.

  \begin{figure} \par\includegraphics[width=6.7cm,clip]{pop2.eps} \end{figure} Figure 1: The percentage of magnetized cells as a function of time. After a short spike, lasting 1000 time steps, the system is stabilized. Dashed line: P = 0.181, solid line: P = 0.179, dotted line: P = 0.177. The method used to calibrate the time in our model is presented in Sect. 2.3.1.
Open with DEXTER

The active regions for four different time steps are presented in Fig. 2. (The method used to calibrate our model in space and time will be explained later in this section.) These pictures are the magnetograms produced from our model. All the pictures are taken were the system had reached stabilization (see Fig. 1). The emphasis in these "magnetograms'' is not placed on the detailed and accurate representation of all the observed characteristics, but on the overall statistical properties of the bipolar structures.

  \begin{figure} \par\includegraphics[width=7.6cm,clip]{PHOTOSPHERE.eps} \end{figure} Figure 2: A $256\times 256$ part of the modelled grid corresponding to a $320~{\rm Mm} \times 320~\rm Mm$ part of the solar photosphere, for four different time steps separated by time intervals equal to 25 h. The parameters used are the "standard values'' mentioned in the text. The method used to calibrate our model is presented at the end of Sect. 2.3.1.
Open with DEXTER

We can estimate several observed statistical characteristics of the active regions using a time sequence of magnetograms produced from our model. We start with the size distribution, the distribution of rise times to maximum development (Fig. 3) and the fractal dimension.

  \begin{figure} \par\includegraphics[width=7cm,clip]{size_dist.eps}\par\vspace*{3mm} \includegraphics[width=7cm,clip]{rise_time.eps} \end{figure} Figure 3: a) Size distribution from the simulated magnetograms. b) Distribution of rise times to maximum development.
Open with DEXTER

We define an active region as a cluster of contiguous active cells. Neighboring but not adjacent clusters are considered as different active regions. Wherever in this work we needed to use a cluster counting method to derive our statistics, we used the algorithm presented by Stauffer & Aharony (1985), modified suitably for each case. The size distribution function of the clusters of the active cells using the standard values for the free parameters was approximated by a least square power law fit of the form

\begin{displaymath}N(s) \sim s^{-b}, \end{displaymath} (1)

where s is the total area of the active regions in pixels, i.e. the number of pixels making up an active region, and $b=1.89\pm 0.03$ up to a typical area of 150 pixels (Fig. 3a). Above this value, we observe an exponential roll-over. This plot resembles remarkably well the observed distributions reported by Harvey (1993) and Meunier (2003). The index of the slope is almost identical to that found by Harvey (1993). We make sure that the size of each active region is measured only once in its lifetime by using magnetograms that are separated by very long time intervals. Every active region is thus measured once at a random time during its lifetime. The size distribution shows a sharp change to an exponential distribution for values around 150 pixels in our model.
  \begin{figure} \par\includegraphics[width=5.7cm,clip]{flux.eps}\vspace*{2mm} \in... ...s}\vspace*{2mm} \includegraphics[width=5.7cm,clip]{max_mean_B.eps} \end{figure} Figure 4: a) Frequency distribution of $\vert\Phi \vert$; the power law index is $1.5\pm 0.05$. b) Frequency distribution of flux density; the power law index is $1.77\pm 0.06$. c) Frequency distribution of $\vert B_{\rm max}\vert$; the power law index is $1.16\pm 0.04$. d) Frequency distribution of $\vert B_{\rm max}/B_{\rm mean}\vert$; the power law index is $1.8\pm 0.06$.
Open with DEXTER

Comparing this graph with the data, we can calibrate the x-axis, since 150 pixels in our scale are equivalent to $400~\rm Mm^2$ in the data, or 1 pixel corresponds to $2.67~\rm Mm^2$. Seiden & Wentzel (1996) calibrated the time scales involved in their model by using the observed time to maximum development of the active regions and compared this result to the percolation model. We have used the same method to calibrate our model. The time of raise to maximum development is defined as "rise time''. The distribution of rise times (Fig. 3b) to maximum development begins with a power law up to 21 timesteps and then becomes exponential. The observations analyzed by Harvey (1993) show the same form with a power law extending up to 1.1 days. So the time scale corresponds to 1.26 h/timestep. Therefore, our phenomenological model is calibrated in space and time through the existing data. We have used this result in Fig. 2.

The fractal dimension of the set of active cells in the model has been measured using the method of box counting (Falconer 1990). The image is covered with coarse-grained boxes of uniform scale L. If N(L) is the number of boxes containing at least one active cell, then the limit

\begin{displaymath}D_{\rm F}=\lim_{L \to 0}{\frac{\ln N(L) }{\ln(1/ L)}} \end{displaymath} (2)

is called the box counting dimension or simply the fractal dimension. We divide the grid into boxes of linear size L, and varying L we plot the number of boxes in which at least one active cell exists, N(L), versus L in log-log scale. Then the absolute value of the slope of a straight line fit is the practical definition of fractal dimension. A least square fit of the box-counting curve of the form $N(L) \sim L^{-D_{\rm F}}$ yields $D_{\rm F}=1.5\pm 0.1$. This result is in agreement with the observations (Meunier 1999) and serves as a constraint for our model.

A detailed statistical analysis of four more parameters related to the observed magnetic field at the photosphere, namely the absolute value of the total flux ($\vert\Phi \vert$) of each active region, the flux density ($\vert\Phi\vert/A$), the maximum magnetic field ( $\vert B_{\rm max}\vert$), and the ratio of the maximum magnetic field to the mean magnetic field ( $\vert B_{\rm max}/B_{\rm mean}\vert$) have also been performed. The frequency distributions of the four parameters exhibit power laws for small values of the parameters with indices $1.5\pm 0.05$, $1.77\pm 0.06$, $1.16\pm 0.04$ and $1.8\pm 0.06$ respectively, and well defined exponential turnovers for large values at 100, 40, 40 and 50 units respectively (see Fig. 4). These results resemble remarkably well the observations (Meunier 2003). Comparing the results from the model with the observations we can calibrate the unit of the magnetic field for the percolation model. The distribution of the maximum magnetic field exhibits a power law up to 40 units of magnetic field and then becomes an exponential. In the data presented by Meunier (2003), this turnover occurs at 1000 Gauss. So 1 unit of magnetic field of the percolation model corresponds to 25 Gauss.

2.3.2 Sporadic energy release
The next important ingredient of our model is the cancellation of magnetic flux due to collisions of oppositely polarized magnetic flux tubes in motion. This leads to the release of energy, whose amount we assume to be proportional to the difference in the square of the magnetic flux before and after the event. The area used in our model is assumed to be approximately constant before and after the collision. The realistic modelling of the sharp discontinuities developed during the collision of the fibrils is beyond the scope of this article (see Priest et al. 2002; Longcope & Kankelborg 1999). In Fig. 5a, we plot the released energy $E(t)=\frac{B(t)^2}{8\pi}\langle V \rangle$ as a function of time. $B^2(t)/8\pi$is the dissipated magnetic energy and $\langle V \rangle \sim \mathcal{L}^3$ the reconnection volume which is constant, with characteristic length $\mathcal{L}\approx \rm Mm$ = the linear cell dimension (the smallest dimension available in our model). Figure 5b shows the energy distribution of the recorded "sporadic energy release": It follows a power law, $f(E)\sim E^{-a}$, with $a=2.25 \pm 0.09$, up to 1029 erg. In Fig. 5c, we plot the distribution of the waiting times. To determine them, we used an artificial threshold for the small explosions that form a continuous timeseries and we record only the "important'' events, i.e. the ones above the "observational'' threshold. The explosive events have no duration, they last one timestep, thus we measured the peak-to-peak waiting times between subsequent "important'' events. The waiting time distribution exhibits a clear power law shape above a certain time scale. The index of the slope for waiting times above 50 h is $2.71 \pm 0.09,$ using a threshold of $8\times 10^{28}~\rm erg$ arbitrary units in energy. It is interesting to note that if we change the artificial threshold level in our data, the shape of the distribution remains unchanged, but the index of the power law tail changes (see Fig. 5d).

A variety of observational results have been confirmed with our model using a unique set of the free parameters reported initially (the standard values). We now discuss the 3-D extrapolation of the simulated magnetograms using the constant $\alpha$ force-free approximation (see Alissandrakis 1981; Démoulin et al. 1997).

  \begin{figure} \par\includegraphics[width=6.2cm,clip]{energy_time_serie.eps}\vsp... ...raphics[width=5.8cm,clip]{wait_peak_thr_slope.eps}\hspace*{2mm}} \end{figure} Figure 5: a) A time series of the released energy, b) the energy distribution, c) the peak to peak waiting time distribution, d) the characteristic exponent of the waiting time distribution as a function of the threshold used.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=7.5cm,clip]{FIELD_LINES.eps}\par\includegraphics[width=7.5cm,clip]{carpet_sheets.eps} \end{figure} Figure 6: a) The 3D topology of the magnetic field and b) regions in the 3D space with steep gradient in the horizontal magnetic field component ( $\Delta \vec{B}_\perp/\vec{B}_\perp>0.2$). The boxes used have dimensions $320{\rm (Mm)}\times 320({\rm Mm})\times 80\rm (Mm).$
Open with DEXTER

3 Force-free magnetic fields

The magnetograms derived from our simulations (see Fig. 2) are used as boundary conditions. We assume that the magnetic field is force-free,

 \begin{displaymath} \nabla \times \vec{B}= \alpha \vec{B}, \end{displaymath} (3)

and that the field is decreasing with height and vanishes at infinity. We also take $\alpha$ to be constant everywhere on the 2-D grid.

Equation (3) is linear and allows us to obtain solutions with the use of Fourier Transforms (FT). We take the z-axis to be perpendicular to the surface of the sun. The magnetogram gives us the vertical component of the magnetic field on the surface (z=0plane). We assume that the Fourier transform of our solutions is decreasing exponentially with height:

\begin{displaymath}\hat{ \vec{B} }(u,{v},z)=\exp(-kz) \hat{ \vec{B}}(u, {v},0) \end{displaymath} (4)

where u and v are the variables in the Fourier transform domain and the "hatted" quantities are the FT of the solutions. We finally end up with the equation which is solved numerically (see Alissandrakis 1981; Démoulin et al. 1997):

\begin{displaymath}\hat{\vec{B}}(u,{v},z)=\hat{\vec{G}}(u,{v},z) \hat{B}_z (u, {v},0)\end{displaymath} (5)

where


\begin{displaymath}\hat{G}_x=-{\rm i}(uk-{v} \alpha) \exp(-kz)/2 \pi q^2 \end{displaymath}


\begin{displaymath}\hat{G}_y=-{\rm i}({v}k+u \alpha) \exp(-kz)/2 \pi q^2\end{displaymath}


\begin{displaymath}\hat{G}_z=\exp(-kz) \end{displaymath} (6)

and q2=u2 + v2 , $k=\pm (4 \pi q^2- \alpha^2)^{1/2}$.

We use the percolation model to estimate the vertical component of the magnetic field at the boundary and obtain the three components of the magnetic field in the 3-D space above the modelled photosphere. We took care that the simulated magnetograms we used for extrapolation were flux-balanced.

For the visualization of our results, we use the magnetogram as the lower boundary of our image and we choose the points on the magnetogram with the most intense magnetic field to be the starting points for the magnetic field lines to be drawn (see Fig. 6a).

The free parameter $\alpha$ is responsible for the curling of the field lines. As $\alpha$ increases the field lines get more and more distorted and curled. We use a canonical value for $\alpha=0.005$ in this article (since $\alpha$ has the dimension of inverse length, we measured the length in terms of the radius of a typical sunspot on the boundary, and we have chosen $\alpha$to be a fraction 0.3 of this inverse length unit (see Sakurai 1981)). Problems related to the distortion at the edges due to aliasing and their solution have been discussed in the literature (see Alissandrakis 1981; Démoulin et al. 1997). A typical extrapolation of the magnetogram produced by our percolation model is given in Fig. 6a. We should emphasize once again that it is beyond the scope of our study to give an accurate and detailed representation of the coronal active region magnetic field topology. The emphasis is placed on the statistical properties of the structures formed.

  \begin{figure} \par\includegraphics[width=8.8cm,clip]{vector_grad.eps} \end{figure} Figure 7: The magnetic field topology at four different heights: the "green carpet'' represents the z-component of the magnetic field while the vectorgrams represent the horizontal components. The yellow contours are the regions were steep horizontal magnetic gradients appear.
Open with DEXTER

3.1 Horizontal magnetic field discontinuities

We next identify the regions were the horizontal magnetic field $B_\perp^2=B_x^2+B_y^2$ forms steep gradients. For each point of the 3-D space above the magnetogram, we compute the relative gradient of the horizontal component of the magnetic field:

\begin{displaymath}\frac{\Delta B _{\perp}}{B_{\perp}} \equiv \frac{\left\vert B... ...+B_{\perp i,j-1}}{4}\right\vert}{\vert B_{\perp i,j}\vert}\cdot\end{displaymath} (7)

In Fig. 6b we plot only the region where $ \frac{\Delta B _{\perp}}{B_{\perp}}$ is above a critical value. In this specific plot the critical value used was $\sim$0.2 (see Parker 1988). We believe that at these steep magnetic gradients the plasma may be unstable and the probability that a few of these regions initiate fast and sporadic energy release is much higher compared to the rest of the active region. Let us assume that the discontinuity will drive a local current over a critical value ( $j>j_{\rm c}=n_0ev_{\rm c})$ i.e.

\begin{displaymath}\nabla \times B\sim \frac{4 \pi}{c} n_0ev_{\rm c}\end{displaymath}

or

\begin{displaymath}\frac{\triangle B_\perp}{L} \sim \frac{4 \pi}{c} n_0 e v_{\rm c}.\end{displaymath}

Deviding both sides by B we have

\begin{displaymath}\frac{\triangle B_\perp}{B} \; \sim \frac{4 \pi}{c} \frac{n_0 e v_{\rm c}}{B}L.\end{displaymath}

Using standard values for the parameters $n_0 \sim 10^9 ~{\rm cm}^{-3},\; B \sim 5 \times 10^3~{\rm G},\; v_{\rm c}\sim 10^7~{\rm cm/s}^{-1} ~(\mbox{sound speed}), \; L\sim 10^6~\rm cm$ we find $\triangle B_\perp/B \sim 0.2{-}0.4$. The unstable current will enhance the local dissipation of the magnetic field and may ignite a flare.

A detailed representation of the magnetic field topology and the regions with steep gradients is shown in Fig. 7. In this picture we plot the intensity of the z-component of the magnetic field for four different heights above the 2-D grid, using the force-free extrapolation method. For each of these plots, we overlay the horizontal component of the magnetic field (represented by the vectors) and the regions where the gradient of the horizontal component is steep, represented by the yellow contour plots.

It is obvious that in the proximity of the photosphere the steep gradients are localized and concentrated in small volumes with a strong magnetic field. As we move away from the photosphere the gradients are spread in larger, cigar-shaped volumes. The statistical properties, in each time step, of these structures will be examined in the next section.

  \begin{figure} \par {\hspace*{2mm}\includegraphics[width=6.9cm,clip]{volume_dist... ...pace*{3mm} \includegraphics[width=7.4cm,clip]{av_energy_dist.eps} \end{figure} Figure 8: a) Volume distribution of the structures of steep magnetic field gradients and b) distribution of the available energy of the same structures.
Open with DEXTER

3.2 Statistical results

Using twenty uncorrelated frames from simulated magnetograms we estimated the 3-D force-free extrapolation of the magnetic field. We impose an arbitrary threshold on the horizontal magnetic field gradient ( $\mid\Delta B_\perp /B_\perp\mid~>0.2$) and record the structures developed in the 3-D simulation box (see Fig. 6b).

Region ( $D_{\rm F})$. The distribution of the volumes of the structures that are created is plotted in Fig. 8a. The frequency distribution of the volumes follows a power-law with index -1.65. Additionally we plot the probability distribution of the "available energy''

 \begin{displaymath} \frac{1}{8\pi}\int{(\Delta B_\perp)^2} {\rm d}V. \end{displaymath} (8)

The distribution follows a power law with index $-1.56\pm 0.05$ up to $3\times 10^{28}$ erg (see Fig. 8b).

Finally, we measure the fractal dimension of the regions with steep gradients in the 3D space. We used the same method of box counting, imposing the necessary changes to apply it in three dimensions. The fractal dimension was estimated to be $1.73\pm 0.05$ (see Fig. 9). We analyzed the fractal nature of the unstable volumes for two reasons: (1) it was already pointed out by McIntosh & Charbonneau (2001) and McIntosh et al. (2002) that the geometric characteristics of the energy-releasing volumes play an important role in the coronal heating estimates; and (2) we believe that there is a correlation between the fractal dimension of the unstable volumes and the fractal characteristics of the active region ( $D_{\rm F})$.

We analyzed also the behavior of the statistical properties of the steep-gradient regions by varying the parameter $\alpha$ from 0.01 to 0.09. The fractal dimension remains constant when we used different values for $\alpha$ (from 0.01 to 0.09) in the force free extrapolation method. We found that the volume distribution of the structures does not depend on $\alpha$ but rather on the height from the photosphere.

The distribution of the available energy inside the structures reported earlier remains also a power-law but the index changes from 1.53 to 1.66 when $\alpha$ changes between 0.01-0.09. The height of the box does not effect the distribution because the energy is mainly concentrated close to the photosphere.

It is not obvious which of the above structures, if any, will eventually release energy, but it seems that patterns formed at the photosphere load the 3-D active region with structures which follow specific statistical laws.

  \begin{figure} \par\includegraphics[width=7cm,clip]{fractal_3D.eps} \end{figure} Figure 9: The fractal dimension of the regions with steep gradients in the 3D space is estimated to be $1.73\pm 0.05$.
Open with DEXTER

4 Discussion and summary

In this article, we use a 2-D probabilistic percolation model for the active region formation. The main feature of the model is the competition between two very important processes (a) the percolation-driven flux emergence; (b) diffusion of magnetic flux due to cancellation or spread of magnetic energy. These processes are controlled by three arbitrary constants, which are calibrated by the fact that only a small fraction ($ 8\%$) of the photosphere is covered by magnetic flux. The choice of these parameters was made once and remained the same throughout this article. The results presented here do not depend sensitively on the choice of the free parameters, e.g. the shapes of the various distributions remained unchanged. Using the above model, we were able to reproduce a long list of observations:
1.
The size distribution of the active regions yielded by the model has the observed shape.
2.
The fractal dimension of the active regions produced by our model is inside the observed range.
3.
Several modelled statistical characteristics of the photospheric magnetic field, e.g the frequency distribution of the total magnetic flux ($\mid$$\Phi$$\mid$), the frequency distribution of the flux density, the distribution of the maximum magnetic field ($\mid$ $B_{\rm max}$$\mid$) and the frequency distribution of $\mid$ $B_{\rm max}/B_{\rm mean}$$\mid$ agree with the observations.
4.
Recording the energy release during the cancellation we were able to give an explanation of the sporadic energy release close to the solar surface and to suggest a possible interpretation for the Ellerman bombs.
5.
The waiting time distribution of the sporadic energy releases follows a well defined power law with an energy threshold-dependent index.
6.
We used a series of magnetograms produced by the percolation model and applying the standard constant $\alpha$ force-free extrapolation techniques we have estimated the 3-D magnetic structures above the photosphere. In these structures, we have located the horizontal magnetic gradients which are above a critical threshold. The analysis of these structures resulted in a power law distribution of the unstable volumes and of magnetic energies included in these volumes. We also estimated the fractal dimension of the unstable volumes, which plays a crucial role in the coronal heating resulting from the stored magnetic energy.
One important question related to the percolation model is the role of the free parameters. Are they crucial for the results reported here? We have explored this question in detail and found that the shapes of the distributions remain unchanged, but the power law indices may change.

The results presented in this article suggest that the observed statistical properties of the magnetic fields at the photosphere (size distribution, fractal dimension, etc.), which are reproduced by the percolation model, produce naturally a loading process which is similar to the one used by Georgoulis & Vlahos (1996). We propose that the energy release in solar active regions follows closely the Self Organized Criticality (SOC) model driven by a power law loading. The percolation model provides the action at the photosphere and the subsequent loading of the active region. The avalanches produced through this specific loading follow the rules of the SOC and yields the observed sporadic energy release in flares. The combined model (Percolation as the driver and SOC providing a model for the energy release) can explain the observed flare statistics.

Our main conclusion is that the activity of the corona is strongly coupled to the detailed balance between the magnetic flux emergence and diffusion in the stochastic photospheric flows. The loading of the coronal active regions and the topology of the magnetic structures formed determine the subsequent evolution of the coronal active regions. We propose that the combination of the photospheric active region formation (using our percolation model) as a dynamic boundary that imposes a power law loading with the energy release processes following the Self Organized Criticality rules can describe the global characteristics of the observed coronal activity (coronal heating, flares, CME).

Acknowledgements
We thank Drs. H. Isliker, A. Anastasiadis and M. Georgoulis for their comments on our article. We also thank the anonymous referee for constructive criticism. E.R. would especially like to thank Prof. F. Moreno-Insertis for his help in understanding the force-free extrapolation techniques during her visit (as part of the European mobility program (Erasmus)) to IAC. This work was in part supported by the Research Training Network (RTN) "Theory, Observation and Simulation of Turbulence in Space Plasmas'', funded by the European Commission (contract No. HPRN-eT-2001-00310).

References


Copyright ESO 2004