Modelling the influence of photospheric turbulence on solar flare statistics
Solar flares stem from the reconnection of twisted magnetic field lines in the solar photosphere. The energy and waiting time distributions of these events follow complex patterns that have been carefully considered in the past and that bear some resemblance with earthquakes and stockmarkets. Here we explore in detail the tangling motion of interacting flux tubes anchored in the plasma and the energy ejections resulting when they recombine. The mechanism for energy accumulation and release in the flow is reminiscent of self-organized criticality. From this model we suggest the origin for two important and widely studied properties of solar flare statistics, including the time-energy correlations. We first propose that the scale-free energy distribution of solar flares is largely due to the twist exerted by the vorticity of the turbulent photosphere. Second, the long-range temporal and time-energy correlations appear to arise from the tube-tube interactions. The agreement with satellite measurements is encouraging.
Parker conjectured that solar flares are driven by the random continuous motion of the footprints of the magnetic field in the photospheric convection parker88; parker89. This conjecture and the experimental observations of power laws in the energy exp1; exp2; exp3; hudson91 and waiting time norandom1; intertime1; intertime2 distributions stimulated a new way of looking at violent bursts. In particular, the fact that the energy distribution is a power law, a property that flares share with diverse physical phenomenastocks, such as avalanches and earthquakes arcangelis1, led to the formulation of flare occurrence models inspired in self-organised criticality (SOC)avalanche; avalanche2; nicodemi; model3; georgoulis1998; vlahos2004; dimi2011; dimi2013; morales2009; morales2008self; kras2002. Although these models reproduce the power-law behaviour in the distribution of flare energies, they predict a Poissonian distribution waiting times, which implies that flares result from an uncorrelated process, contrary to experimental observations norandom1; intertime1; intertime2. This point was stressed by Boffetta et al.norandom1 who, by implementing shell models for turbulence, reproduced the observational power-law decay of the waiting time distribution. However, their exponents are not universal, depending instead on the model parameters.
Based on a different approach, the waiting time distribution can also be reproduced in terms of a piecewise Poissonian process wheatland1. More recently, the existence of correlations between flare energies and waiting times has also been investigated energyintertime1; intertime2; energyintertime2; lucilla2. In particular, Lippiello et al. lucilla2 found that the observed time-energy correlations are not simply attributable to obscuration effects. Here the term obscuration indicates an observational limitation that can be at the origin of the incompleteness of the catalogue, for example, the occurrence of a large flare can hide the detection of smaller flares occurring nearby.
An approach more closely inspired in magnetic reconnection was adopted by Hughes et al.model3, who proposed a dynamical model of solar flares as cascades of reconnecting magnetic loops, with multiple loops that are randomly driven at their footprints and interact with each other. Despite some discrepancy with experimental observations, they showed a relation of the distribution of magnetic loops with a scale-free network, which conceptually supports self-organised criticality.
The formulation of a theoretical model able to reproduce both the flare statistics and the behaviour of time-energy correlations remains unsolved. A complete model would require a fully developed realistic convection zone, a stratified atmosphere above it and the study of the interplay between magnetic fields and flows. Additionally, the model must be three-dimensional, including explicit resistivity in order to control the reconnection between the colliding magnetic fields. At present, a numerical study of such a model is far out of reach due to the excessive requirement in computing time. On the other end, simplified models considering magnetic reconnections based on a purely statistical approach and neglecting completely the photospheric flow, fail in reproducing a variety of observational data avalanche; avalanche2; nicodemi; model3; georgoulis1998; vlahos2004; dimi2011; dimi2013; morales2009; morales2008self; kras2002.
In our study, we present an approach that represents a compromise between these two different scenarios. Here we introduce and study numerically a theoretical model of solar flare occurrence in terms of reconnection of magnetic flux tubes twisted by the photospheric turbulent flow, which reproduces satisfactorily the flare statistics and the behaviour of time-energy correlations. The motion of the tubes in the solar corona is mainly rooted in the photosphere, which is about km thick, corresponding to an extremely thin layer as compared to the solar radius. Therefore, in the model proposed here we consider the photospheric flow as a two-dimensional turbulent system following Kolmogorov scaling.
The turbulent fluid dynamics of the photosphere is simulated through a lattice Boltzmann modelLBE1 on a square lattice of size , with a forcing term that specifically reproduces the Kolmogorov energy spectrum regime (see Methods).
Anchored in the photospheric flow, the footprints of the magnetic flux tubes follow the local velocity field, and are twisted by the vorticity. The magnetic lines are modelled as lines (see green and pink lines in Figure 1) wrapped around the semi-circular flux tubes (see semi-transparent grey tori in Figure 1), forming, when twisted, a spring-shaped bundle. This representation, which is conceptually consistent with previous realistic models for the kink instability hood2, leads to an explicit relation between the length of the spring and the magnetic energy stored in the corresponding flux tube (see Eq. (6) in Methods).
Observational evidence supports the kink instability as the triggering mechanism for magnetic reconnection, and consequently, for solar flare occurrence priest2007magnetic; expkink1; expkink2; kliem2004. The kink instability is an ideal magnetohydrodynamical (MHD) instability, where magnetic reconnection is not a necessary ingredient, as a flux tube can destabilise by converting twist into writhe. However, in an active region (which is our region of interest), the kink instability will trigger a flare (due to the accumulated free energy), and therefore, we assume that each time a kink instability occurs in a magnetic tube, a flare is released. Accordingly, the kink instability of a flux tube occurs as soon as the intensity of its cumulative twist reaches a given critical value. In our model, the tube releases its total energy (vanishing from the simulation), and a new magnetic flux tube is inserted with the initial condition and located at a random position inside the simulated zone (see Methods). The critical twist , at which the magnetic reconnection occurs, can be obtained from stability analysis. Several values have been proposed from numerical simulations, theoretical models hood2 or deduced from experimental observations expkink1; expkink2, ranging from to , depending on the particular plasma conditions in the corona.
Every time a magnetic flux tube reconnects and releases its energy, we implement two possible scenarios: The reconnection heats up the surrounded plasma increasing the local coronal pressure, and therefore, increasing the critical twist of the surrounding magnetic flux tubes hood2. This causes a delay in the reconnection process, which is taken into account here by multiplying the cumulative twist of neighbouring tubes by a positive factor , while keeping the critical twist constant. This procedure implements tubeÂ-tube interactions and in our simulations we consider either random values or a constant value for , obtaining the same results. Conversely, we also consider the case which induces an avalanching process in the occurrence of solar flares. In this case, every time a magnetic flux tube releases its energy it increases the twist of the surrounded tubes, triggering new flares and generating a cascade of events. In both cases, the parameter has the important role of tuning the level of tube interactions, whereas non-interacting tubes correspond to . We are therefore able to detect how interactions affect the flare waiting time distribution. In what follows, we show that tube-tube interaction is the most relevant ingredient to reproduce the correct temporal organisation and the experimentally observed time energy correlations. We will also show that changing the value of the critical twist within the range of the experimentally meaningful values, does not affect the statistics of solar flares. Our model takes into account the well-known SOC mechanism of slow energy storage by changing the flaring threshold of tubes that are neighbours to the reconnected tube. This mechanism is equivalent to tuning the critical tube twist, a procedure usually implemented in SOC models avalanche; avalanche2; nicodemi; model3; georgoulis1998; vlahos2004; dimi2011; dimi2013; morales2009; morales2008self; kras2002. More precisely, by increasing the critical twist () , the photospheric flow has more time to add energy to the respective magnetic flux tubes before they release their energy as a flare. For the case , the similarity with SOC models is more evident, since our model can generate a cascade of events by increasing the twist of surrounded magnetic flux tubes, each time a flare is released. Finally, we should mention that our approach does not account for a number of additional features of flare occurrence and flaring active regions forbes2006cme; aulanier2013, e.g. systematic peculiar flows on top of the heavily suppressed quiet-Sun Kolmogorov flow velocity field, the presence of intense magnetic polarity inversion lines (PILs) in the photosphere, and slip-running magnetic reconnection in flares. However, the excellent agreement between numerical and observational data suggests that these features are not relevant to produce the observed statistical properties, even if they may be necessary for the complete understanding of solar flares.
Peak and Total Energy Distributions
Following the described dynamics, the occurrence time and energy released by each flare quantify the statistical properties of our model and are compared with observational data. We perform extensive numerical simulations for different system sizes , , , , and . For each value of , we first let the fluid evolve without the magnetic flux tubes until it reaches the turbulent regime at Reynolds numbers between and . The Reynolds number is computed by the relation , where is the root mean square velocity, and the kinematic viscosity (values of the Reynolds number and kinematic viscosity for each simulation can be found in the Supplementary Table 1). The fully developed turbulence condition is established when the energy spectrum of the flow follows the classical power-law behaviour with Kolmogorov exponent (see Supplementary Figure 1). At this point, we add magnetic flux tubes at random positions and with random initial energies, and record the occurrence time and energy of flares. From previous studies of satellite data by the Solar and Heliospheric Observatory (SOHO), deviations from the Kolmogorov exponent have been observed depending on the activity of a solar regionkolmogorov_MHD. However, we show that the choice of the exponent of the spectrum of the turbulent flow does not change appreciably our results (see Supplementary Figure 2).
As shown in Figure 2, the distribution of peak flare energies follows a typical power-law behaviour, , with an exponent . As expected, this scaling regime extends up to a cutoff energy that gradually increases with system size . We compare our results with soft and hard X-ray data from the GOES GOES and the BATSE BATSE catalogues, respectively. From the GOES catalogue, only flare events of class C and above (peak flux Watt m) observed between and are considered, which leads to a total of events, covering nearly two solar cycles. In the case of the BATSE catalogue, events in the period between and were considered, with a peak flux larger than counts sec cm. The Kolmogorov-Smirnov test ensures that both samples, numerical and observational data, follow the same scaling behaviour with a -value of (confidence level of ). The excellent agreement between numerical and observational data shown in Figure 2 confirms the validity of the theoretical approach. The exponent of the power law for our numerical results is also in excellent agreement with previous experimental studies energyintertime1; arcangelis1. Notice that in our model for , flares are instantaneous events, therefore we cannot measure separately the peak energy and the total energy associated to each event, or else the energy of a flare is a peak flux energy. Different is the case , where avalanching is induced and events are therefore not instantaneous. In order to implement a unified procedure for models with different , numerical data are compared to peak flux energies from experimental observations. For the case of , one can also study the total energy distribution and the duration of each event. In Figure 3, we can see that the numerical results are in very good agreement with observations, showing the correct duration and total energy distributions. The exponents for the total energy distribution and the duration of flares, and , respectively, are also in good agreement with previous studies georgoulis1998.
Note that our model results in steeper distribution functions for the total energy than for the peak released energy of the modeled events. On the other hand, numerous observational works report the opposite, that is, clearly flatter distribution functions for the total energy. Theoretical works, at least those relying on SOC, seem also to predict analytically that total-energy distributions functions are flatter than those of the peak energy released.
Waiting Time Distributions
We also investigate the statistical patterns of the waiting time, defined as the distribution of time delays between the end of an event and the beginning of the next one. As shown in Figure 4, the numerical results from our theoretical approach are also in very good agreement with experiments. The distribution is not a simple exponential, suggesting that flare occurrence is not a purely uncorrelated Poisson process. In order to closely compare the different numerical and observational catalogues, we have rescaled the waiting times, , by the average event rate in each catalogue arcangelis2, i.e., by the inverse of the average waiting time, , where is the number of events in the respective catalogue. Here and are the times at which the last and first events in the catalogue occur, respectively. As shown in Figure 4, rescaled distributions for numerical and observational data collapse onto a universal curve well fitted by fitting-time, , where and are constants, and denotes the exponent of the power-law regime of the distribution for large waiting times. This result is reasonable when compared with previously reported observational values, fitting-time and norandom1. We have also performed the Kolmogorov-Smirnov test, finding that both samples, numerical and observational data, come from the same distribution with a -value of (confidence level of ).
It is interesting to investigate the role of different ingredients of the theoretical model on the statistical properties of energy released and waiting times, in order to identify the main triggering mechanisms for the occurrence of solar flares. We start by considering that, instead of being driven by the turbulent flow, the magnetic flux tubes might move along purely random trajectories and the cumulative twist is calculated by assigning a random vorticity at each footprint. Results in Figure 5 show that the suppression of the turbulent flow leads to an energy distribution that is exponential rather than a power law. Next, we consider the case where there is a single magnetic flux tube evolving in the turbulent flow, eliminating the possible role of interactions among tubes. We observe in Figure 5 that the power-law regime in the peak energy distribution is recovered. These two results strongly suggest that the ingredient responsible for the power-law in the energy distribution is the turbulent motion of the footprints anchored into the photosphere, and not tube-tube interactions. We finally consider the case of several tubes having different degree of interaction, i.e. either interacting ( and ) or non-interacting () tubes. Results shown in Figure 5 confirm that interactions do not modify the distribution of solar flare energies. Interestingly, models with and without avalanching exhibit the same scaling exponent for the peak flare energy distribution, suggesting that indeed small flares share similar statistical properties with major flares.
We now consider the waiting time distribution for the previous cases. Indeed, in Figure 6, we see that the case of random footprint motion and the case of a solitary tube are well described by a Poissonian distribution (dashed line). This implies that, although essential to the recovery of a scale-free energy distribution, the turbulent fluid flow alone is not able to provide the right temporal organization of solar flare occurrence. If more tubes are considered, the distribution starts to deviate from a Poissonian one. For coexisting but non-interacting tubes (, ), the turbulent flow in the photosphere is able to induce time correlations between them, although not sufficiently to reproduce the observational results. Indeed, the physical correlations are fully recovered only for interacting tubes (, and ). From our results, we can conclude that, whereas the turbulent photospheric flow is the main mechanism responsible for the energy distribution, the interaction between magnetic tubes is what introduces the right temporal correlations in the process.
We further investigate the statistical features of time-energy correlationslucilla2. For each catalogue, we analyse how the flare energies are organised in time, by evaluating the probability that a flare with energy is followed by a flare with energy larger than under the conditions that their temporal distance is smaller than a certain threshold , . For each catalogue, this probability fluctuates wildly due to statistical noise. Therefore to eliminate this noise, we evaluate the same probability also in a synthetic catalogue generated by reshuffling the flare energies with respect to their occurrence time, such that energy and time are uncorrelated by construction. We then consider the difference between the conditional probabilities, , evaluated in the two data sets. This difference is different from zero only if significant time-energy correlations are present in the original catalogue. In particular, if is larger than zero, it is more likely to find two consecutive flares satisfying both conditions ( and ) in the real rather than in the reshuffled catalogue (see Methods). By using the same technique we also compute the conditional probability difference to observe a flare energy larger than a given threshold after an waiting time smaller than . We consider the behaviour of both conditional probability differences for a range of parameters , and .
In Figure 7 we see that the probability differences are very well described by our model with . In particular, for both, numerical and observational results, is different from zero beyond error bars. This implies that it is very likely that for close-in-time flares the second one will have slightly larger energy than the previous one (the maximum is for ), as far as their separation in time is shorter than approximately hours. These energy correlations decrease as the temporal separation increases. Conversely, it is very unlikely to observe in real catalogues close-in-time flares where the second one has a smaller energy ( for )). Furthermore, in Figure 7b curves for are different from zero beyond error bars and decrease with increasing . This implies that the probability to find couples of successive flares with the second flare having energy higher than decreases, if one includes events separated by a larger in the analysis. For large numerical results deviate from the observational ones. This can be due to the finite size of our simulation, which imposes an upper limit to the flux tube sizes, and therefore limits the maximum energy of flares (see also Figure 2, where finite-size effects in the energy distribution are analysed). Note that the agreement between observational and numerical data is very good, suggesting that, in fact, the turbulent flow and magnetic flux tube interactions induce correlations between energy and time in the occurrence of solar flares.
It is important to notice that agreement with observational data is only obtained for interacting tubes (), not for non-interacting tubes, single tube, and random motion of tube footprints, where and are found. For the case , we do not obtain full quantitative agreement with time-energy correlations measured in experimental catalogues (see Figure 8). Data still predict that the next flare statistically has slightly larger energy than the previous one but in a narrower range and with a smaller probability. Moreover, the case is not able to reproduce the anti-correlations for . These results suggest that the relaxation mechanism corresponding to seems to be more appropriate to reproduce time-energy interaction, as compared to the avalanche process characterised by .
Finally, we have performed extensive simulations varying every parameter of the model, namely, the critical twist , the magnetic flux tube density via , and the size ratio of the magnetic flux tubes (where and are the inner and outer radii of the magnetic flux tubes). For all cases, results for flare statistics and time-energy correlations remain unchanged as shown in the Supplementary Figures 3, 4, and 5.
It is well accepted in the literature that solar flare occurrence is a process driven by magnetic reconnection. Since high quality satellite data became accessible, several studies have evidenced that the statistics of this phenomenon is complex: It exhibits scale-free energy distribution and a nontrivial waiting time distribution. A number of theoretical models attempted to reproduce such statistics with different approaches. This study implements magnetic reconnection in a model framework that enables us to test the role of the different physical ingredients on observed statistical patterns. In particular, we have shown that the energy distribution are ruled by the turbulent features of the flow in the photospheric plasma.
More precisely, that, if tube footprints simply diffuse in the corona in absence of turbulent flow, the observed distribution of flare energies would be exponential, namely it would exhibit a characteristic flare size. Moreover, the turbulent flow alone is not sufficient to fully reproduce the statistical patterns of real data. Indeed the evolution of a single tube or several non-interacting tubes in the corona, exhibiting the observed energy distribution, is not sufficient to account for temporal correlations.
The detailed analysis of the energy organisation in time indicates that turbulence and tube interactions are the essential physical ingredients controlling solar flare occurrence. Proving time-energy correlations is the first step towards any forecasting model. This could be formalised by implementing the phenomenological laws, as done for earthquakes Ogata (ETAS model), and would open a novel field of investigation.
Evaluation of the turbulent flow
For modelling the two-dimensional turbulent flow, we have used a two-dimensional lattice Boltzmann model of size with a cell configuration ( dimensions and discrete velocity vectors)LBE1. In order to induce turbulence, we have included the following forcing term in the 2D Navier-Stokes equation
where is a constant, , and varies in time such that each value of is used during an interval of time , and then is increased by one. As an initial value, we take . Here, is a tuning parameter to control the spectrum of the energy of the turbulent flow (in our simulations, ). The coefficient defines the minimum wave number (due to space discretisation limitations). This forcing term satisfies the incompressibility condition . The kinematic viscosity of the fluid is set to and the forcing coefficient to . As initial conditions for the fluid, we choose density , and velocity . Furthermore, we also impose a large scale dissipation mechanism to avoid vortex condensation2Dturbulence. All the values are in numerical units.
Once the fluid has reached the turbulent regime, we insert magnetic flux tubes in random positions. The position of the respective footprints for each tube , denoted by for the positive footprint and for the negative one, is a function of time and evolves as,
where is the velocity of the fluid at position .
Then, we can define and as the cumulative twist in the positive and negative footprint, respectively, evolving according to the equation,
Note that the component used to twist the magnetic flux tubes is the z-component of the vorticity, since the velocity lies on the two-dimensional space.
As initial conditions, the flux tubes have an outer radius of cells and zero twist, . If the positive footprint of a tube comes very close to its negative partner (less than two lattice nodes), we reset the tube to the initial condition (initial length and zero twist) and relocate it at a random position inside the simulated zone. The vanishing of these tubes can be seen as small reconnection processes with negligible released energy. Because of space discretisation in our numerical simulations, the position of each footprint is, in general, not located at a fluid grid point, therefore we use bilinear interpolation to calculate the velocity at the footprint position. Note that the motion of the footprints, see Eq. (2), as well as the twist, see Eq. (3), are additive (See Supplementary Figure 6), in the sense that they systematically inject electric currents and associated magnetic energy in the system.
The magnetic field lines are modelled as lines wrapped around semi-circular flux tubes, forming, when twisted, a spring-shaped bundle (see Figure 1). Therefore, we can assume that the total energy of a tube is given by the length of the magnetic line, which depends on the twisting and the size of the semi-circular tube. Thus, when a flux tube is not twisted, its energy equals (here and throughout the following calculation, we have omitted the proportionality constant to get the right units of energy). On the other hand, if a flux tube is twisted, the spring-shaped bundle can be parametrised by
where is the cross section radius of the semi-circular tube. The value depends on as follows: , where is a constant that controls the number of turns that the magnetic line makes around the semi-circular tube. In this coordinate system , the photosphere is located at the plane . The length of this parametric curve is given by the integral,
According to the kink instability theory hood2, the twist is defined as , where and are the tangential and perpendicular components of the magnetic field to the plane at the footprint (). Therefore, the ratio is equivalent to the ratio between the and components of the derivative of the parametric curve, with , and we can conclude that is the twist. In our model, denotes the cumulative total angle due to the vorticity of the fluid, .
The integral in Eq. (5) does not possess an analytical solution. However, we can assume that , obtaining a very simple expression:
Note that this expression is identical to the expression for an untwisted flux tube for any values of and .
Once a magnetic flux tube reaches the critical twist , the tube releases its entire energy and vanishes. In order to keep the tube density in a stationary state and produce a sufficient statistics, a new tube is placed at a random position inside the simulated zone with the initial condition ( and cells). When several magnetic flux tubes reach the critical twist within the same temporal interval , we sum the energies of the tubes into a single event. For the case , we have also evaluated the distributions keeping simultaneous events separate (see Supplementary Figure 5) and verified that the main properties of solar flare statistics remain unchanged. For the case of , since flaring occurs through an avalanching process, we perform the measurement of events as follows. Once a flare occurs, we stop the fluid and observe if it triggers other flares. We measure the peak flux energy as the largest flare that occurs within the avalanche process, and the total energy as the sum of all flares. The duration of flares is taken as the total number of released flares. Afterwards, the dynamics of the photospheric flow is restarted. Our model cannot accommodate helicity conservation during the magnetic reconnection process berger1999. Therefore, it is assumed that each helical kink instability ejects the unstable flux tube out of the simulation volume to infinity (physically, that would mean that each flare is eruptive). However, it seems that this effect is not relevant to reproduce the solar flare statistics.
Note that the energy stored by a tube scales with the tube length and therefore has an upper cutoff controlled by the system size. We have also run the simulations for different initial conditions, finding that our statistical results remain unchanged. In particular, for the cases where (but not necessarily ), we have solved numerically the integral in Eq. (5) considering terms up to order .
We have implemented our numerical code using CUDA C. The simulations for a fixed set of parameters run three weeks on a cluster of graphic cards, Nvidia Tesla M2075, each one containing GPU cores.
Conditional probability analysis
Each flare in the numerical and observational catalogues is characterized by its starting time and its peak-flux energy . From each catalogue we evaluate the conditional probability to find the energy of the next flare () being larger than times the energy of the previous flare (), if their temporal distance, , is smaller than a certain threshold, . For comparison, the same conditional probability is evaluated from a reshuffled sequence of the same energy time series. In such synthetic catalogues we expect that flare energies and occurrence times are totally uncorrelated. Keeping and fixed, we compute the quantity for independent realisations of the reshuffled catalogue, obtaining an ensemble of values which follows a Gaussian distribution with mean value and standard deviation . We then define . If the absolute value , a significant difference in the number of pairs of sequential energies satisfying both conditions exists between the real and the reshuffled catalogue. By using the same technique we also compute the conditional probability difference .
Acknowledgements.We acknowledge financial support from the European Research Council (ERC) Advanced Grant 319968-FlowCCS, and the Brazilian agencies CNPq, CAPES and FUNCAP. MM and AK would also like to acknowledge many fruitful discussions with Fabrizio Lombardi and Micha Wasem.
All authors conceived and designed the research, analyzed the data, worked out the theory, and wrote the manuscript.
Competing financial interests: The authors declare no competing financial interests.