Multiscale statistical analysis of coronal solar activity

Multi-filter images from the solar corona are used to obtain temperature maps which are analyzed using techniques based on proper orthogonal decomposition (POD) in order to extract dynamical and structural information at various scales. Exploring active regions before and after a solar flare and comparing them with quiet regions we show that the multiscale behavior presents distinct statistical properties for each case that can be used to characterize the level of activity in a region. Information about the nature of heat transport is also be extracted from the analysis.


Introduction
The increasing number of space telescopes and space probes that provide information about phenomena occurring in the space is yielding enormous amount of data that need to be analyzed to get information about the physical processes taking place. In particular, images of the sun obtained from the Solar Dynamics Observatory (SDO) instruments and the new Interface Region Imaging Spectrograph (IRIS) mission have a remarkable high resolution that allow studies of the sun with great detail, not available before. The analysis of the images usually involves some processing if one is to extract information not readily obtained from the raw images. Some techniques used for image analysis and feature recognition in the sun have been described by [1]. A useful tool that gives information about the physical conditions in the emitting region is the emission measure (EM) which is obtained from images in several wavelengths; this in turn provides the local temperature of the emitting region (e.g. [2]). When the EM is computed for all the pixels in a given region of the sun, the temperature maps can be obtained. There are some computational tools developed in Solar-SoftWare (SSW) that take images from different instruments at the available 1 arXiv:1601.07681v1 [astro-ph.SR] 28 Jan 2016 filters to produce EM and temperature maps (e.g. for Yohkoh, SOHO, TRACE, SDO).
The images can be processed for pattern recognition using techniques such as wavelet analysis that sort the structures in terms of the scaling properties of their size distribution [8]. Another method, which is also frequently used in image compression, is the singular value decomposition (SVD) of the pixel matrix. This method is the basis of the more general proper orthogonal decomposition (POD) which has been applied in many contexts including: solar physics [17,18] image processing [16]; and turbulence models [12]. Application to plasma physics relevant to the present work include compression of Magnetohydrodynamics data [6]; detection of coherent structures [9], and multi-scale analysis of plasma turbulence [10,11]. In this work we apply POD techniques to a time sequence of images representing solar temperature maps. The specific method, referred to as Topos-Chronos [5,10], is used to perform spatio-temporal multi-scale analysis of different solar regions. Of particular interest is the analysis of solar flare activity. To this end we perform an analysis of temperature maps corresponding to a region of solar flare activity, and compare the results with the analysis done in the same region before the solar flare and in a "quiet" region where no flares were detected during the time of observation. Our POD multi-scale analysis is accompanied with a statistical analysis of the probability distribution function (PDF) of temperature fluctuations. In particular, in the search for statistical precursors of solar flare activity, we show that the small scale temperature fluctuations in pre-flare states exhibit broader variability (compared to the flare and quiet sun states) and the corresponding PDFs exhibit non-Gaussian stretched-exponential tails characteristic of intermittency. We also present a study heat transport during solar flare activity based on a simplified analysis that allows to infer the transport coefficients (advection velocity and diffusivity profiles) from the data.
The rest of the paper is organized as follows. In section 2 the fundamentals of POD methods are presented and the features of the Topos-Chronos method as the optimal representation of a truncated matrix are explained. Section 3 describes how the data for temperature maps are obtained for a solar flare event and for other solar regions. The results of applying the POD methods to the data are described in section 4 highlighting the statistical features of the multiscale analysis. This leads to the identification of intermittency in an active region previous to a flare. As a further application, section 4.2 presents the evaluation of heat transport coefficients based on the full temperature maps and on POD. Finally, in section 5 we give the conclusions of the analysis and an appraisal of the POD methods applied to solar physics.
2 Description of POD methods used in the work.
In this section we review the basic ideas of the POD method and its application in separating spatial and temporal information as used in the present work. 2

SVD
The Singular Value Decomposition (SVD) is, generally speaking, a mathematical method based on matrix algebra that allows to construct a basis in which the data is optimally represented. It is a powerful tool because it helps extract dominant features and coherent structures that might be hidden in the data by identifying and sorting the dimensions along which the data exhibits greater variation.
The SVD of the matrix A ∈ R m×n is a factorization of the form, where U ∈ R m×m and V ∈ R n×n are orthogonal matrices (i.e. U T U = V V T = I) and Σ ∈ R m×n is a rectangular diagonal matrix with positive diagonal entries. The N diagonal entries of Σ are usually denoted by σ i for i = 1, · · · , N , where N = min{m, n} and σ i are called singular values of A. The singular values are the square roots of the non-zero eigenvalues of both AA T and A T A, and they satisfy the property σ 1 ≥ σ 2 ≥ · · · ≥ σ N . [13,15] Another useful mathematical expression of SVD is through the tensor product. The SVD of a matrix can be seen as an ordered and weighted sum of rank-1 separable matrices. By this we mean that the matrix A can be written as the external tensor product of two vectors: u ⊗ v or by its components: where r = rank(A). Here u i and v i are the i-th columns of U and V respectively. Equation (2) is convenient when one wants to approximate A using a matrix of lower rank. In particular, according to the Eckart-Young Theorem [7], given eq. 1, the truncated matrix of rank k < r = rank(A) is the optimal approximation of A in the sense that A−A (k) 2 = min rank(B)=k ij is the Frobenius norm. In our analysis the matrix A(x, y) is the scalar field representing the temperature at each point x, y at a given time. In this sense we have temperature maps that are represented by the 2D spatial rectangular grid (x i , y j ) where i = 1, . . . , N x and j = 1, . . . , N y . The maximum number of elements in the SVD expansion (eq. (3)) is, therefore, k = min{N x , N y }. These maps change as a function of time.

Topos and Chronos
Topos and Chronos is the name given to the method based on the POD in which time variation is separated from the space variation of the data [5,9].
To illustrate this method, consider a spatio-temporal scalar field A( r, t) representing the temperature T, where r = (x, y). We introduce a spatial grid (x i , y j ) with i = 1, . . . , n x and j = 1, . . . , n y and discretize the time interval t ∈ (t in , t f ), as t k with k = 1, . . . , n t , t 1 = t in and t nt = t f . We represent the data as a 3D array A ijk = A(x i , y j , t k ). Since the SVD analysis applies to 2D matrices, the first step is to represent the 3D array as a 2D matrix. This can be achieved by "unfolding" the 2D spatial domain of the data into a 1D vector, i. e. (x i , y j ) → r i with i = 1, . . . , n x × n y , and representing the spatio-temporal data as the (n x n y ) × n t matrix A ij = A(r i , t j ), with i = 1, . . . , n x × n y and j = 1, . . . , n t .
The singular value decomposition of A(r i , t j ) is given by the tensor product expression, where N * = min[(n x n y , n t )]. In this sense, the POD represents the data as a superposition of separable space time modes.The vectors u k and v k satisfy the orthonormality condition given by, The n x n y -dimensional modes u k are the "topos" modes, because they contain the spatial information of the data set, and the n t -dimensional modes v k are the "chronos" since they contain the temporal information. From eq. (3), we define the rank-r optimal truncation of the data set A r as, where 1 ≤ r ≤ N * . The matrix A (r) is the best representation of the data due to the fact that, by construction, the POD expansion minimizes the approximation error, A − A r 2 .
3 Temperature data maps.
Using the combination of six filters for the coronal emission as measured by the Atmospheric Imaging Assembly (AIA) instrument [4,14] of the SDO, the dual maps for the emission measure (EM) and temperature are obtained, following the techniques developed by [3]. The first step is to select appropriate events that contain the phenomenology of interest. In our case, we focus on the propagation of a heat front associated with an impulsive release of energy, such as a solar flare. Therefore, we look for events satisfying the following criteria: (1) The AIA instrument consists of four detectors with a resolution of 4096×4096 pixels, where the pixel size corresponds to a space scale of ≈ 0.6 (≈ 420 km). It also contains ten different wavelengths channels, three of them in white light and UV, and seven in EUV, whereof six wavelengths (131, 171, 193, 211, 335, 94Å) are centered on strong iron lines (Fe VIII, IX, XII, XIV, XVI, XVIII), covering the coronal range from T ≈ 0.6 MK to 16MK. AIA records a full set of near-simultaneous images in each temperature filter with a fixed cadence of 12 seconds. More detailed instrument descriptions can be found in [14] and [4].
For each wavelength there is a corresponding temperature (T λ ) for the emitting region, so the images of the different intensities, F λ (x, y), can be used to obtain maps of the coronal temperature. These intensities represent the measured photon flux, which is determined by both the physical conditions in the emission region and the filter response. The former is usually represented by the so-called differential emission measure (DEM) which is used to obtain the temperature distribution of the plasma averaged along the line of sight.
The DEM distribution [dEM/dT ] is typically given by dEM (T ) = n 2 dh/dT [cm −5 K −1 ], where n(h(T )) is the electron density at height h and with temperature T . The DEM distribution DEM= d[EM (T, x, y, z)]/dT can be reconstructed from the six filter fluxes [F λ (x, y)], but [3] developed a method to compute the peak emission measure (EM) and the peak temperature T p for each pixel based on a simple representation of the DEM by a Gaussian function of temperature, for each filter. Then, model fluxes are computed from the expression where R λ (T ) are the filter response functions at a given wavelength, which are fitted to the observed flux F obs λ (x, y). This determines the best-fit values for EM and T p thus obtaining the temperature maps. Here we use the Automated Emission Measure and Temperature map routines developed by [3] that implement the method of multi-Gaussian functions based on six-filter data, following a forward-fitting technique. The result of applying this post-processing methodology on the data in Fig. 1 is shown in Fig. 2. The temperature map covers a temperature interval in the range log(T ) = 5.7 − 7.0 shown in the vertical bar. It highlights the hot regions while the EM map is more sensitive to the plasma density, since it is proportional to the square of the electron density (independent of the temperature). The temperature maps were computed for a time sequence covering 16.4 minutes with a cadence of 12 seconds, as provided by the AIA data. In total, 82 maps were generated.

Results of the analysis
The data generated in the previous section combining the emission at six wavelengths can be further processed to extract structures and dynamical features. The time evolution of the data space array, A(x, y, t), is analyzed with the POD 6 method described in Section 2. We first focus on the heat front that propagates from left to right in the region close to the flare site. From these maps we can extract information concerning: (1) the nature of the heat transport and (2) the physical conditions in the region. In particular, we seek to obtain information from the analysis about the energy distribution in the different spatial scales and how this compares with other regions with no flares. Fig. 3 presents the region where a heat front is moving at selected times in the full time interval; the vertical lines indicate the position of the leading-front at the initial and final times. We limit the attention to this region which is a rectangle with dimensions N x × N y with N x = N y = 32 pixels. The POD method can be applied in a systematic way to separate time and space features and study them with a multi-scale analysis, or to determine the underlying features of heat transport.
The same POD analysis performed on the temperature maps in Fig. 3 is also applied to two other temperature maps corresponding to (a) exactly the same physical region but about one hour before the appearance of the flare (to which we refer as pre-flare) and (b) a different region in the quiet sun for the same time interval as the first map (referred to as quiet-sun or QS). These other analyses are made in order to compare the properties of the three cases and be able to determine the features associated with flare activity in the solar corona.

Multiscale statistical analysis using Topos-Chronos decomposition
As mentioned in subsection 2.2, the POD represents the data as a superposition of separable space-time modes. However, if we plot the absolute value of the chronos mode as function of the indices j and k there is no evident regularity or periodicity found that may indicate some kind of cascade process. A similar behavior is found for the topos modes in which large scale structures are seen for low k while a granulated structure with small scales is found for high k. Thus, there is some correlation between small scales and high frequencies for high k modes and large scales and low frequencies for low k modes. This correlation can be further studied using Fourier analysis to extract the characteristic spatiotemporal scales of each POD mode. For the Topos modes, we first transform the one-dimensional vector r i back to space coordinates x l and y m and compute the two-dimensional Fourier transform,û k (κ x l , κ ym ). The characteristic length scale of the topos rank-k mode is then defined as the mean length scale of the corresponding Fourier spectrum, λ(k) = 1/ κ k , where κ k Odd columns: spatial modes u k (r i ); even columns: temporal modes v k (t j ), for k=1,2,3,4,20,30,40,46. One unit in the time axis equals 12 s. is defined as: Using the same procedure, we associate a characteristic time scale τ to each chronos v k (t m ) mode using τ (k) = 1/ f k with: wherev k (f m ) is the Fourier transform in time of v k (t m ). Figure 7 shows the characteristic length scale λ(k) versus the characteristic time scale τ (k) for all ranks. The tendency for increasing time scale with increasing space scales is apparent, although there is some dispersion for small space-time scales. This means that there is some sort of cascading where heat is transfered from large to small spatio-temporal scales. This is similar to the cascading process of turbulent systems (e.g. Futatani et al., 2009). The scaling can be fitted by a power law of the type λ ∼ τ γ and in this case γ = 0.11. For a diffusive-like process γ = 1/2, thus this results suggests that a sub-diffusive heat transfer may be taking place. The topos-chronos plots for the pre-flare and quiet-sun cases in Figs. 5 and 6 do not show a clear correlation of the time and space scales. One can notice that only for rank 1 there are dominant large spatial scales while for higher ranks there are only small scales present. In the chronos diagrams the dominant small frequencies seen at low k have higher frequency oscillations superimposed. In order to determine if there is any correlation between spatio-temporal scales we make the same Fourier analysis as for the case with flare. The resulting length and time scales for all ranks are shown in Figure 8 for the two regions, pre-flare and QS. It is clear that for these cases no correlation is found and thus, it is not possible to attribute any cascading-like process, as in the region affected by the flare.
Further information can be obtained from the reconstruction error in the POD representation, defined as the difference δT ti where T is the original temperature map at a given time t i , and T (r) ti is the reconstructed map up to rank-r at the same time t i . Since a measure of the total energy content in the data is given by E = ij A 2 ij = N * k=1 σ 2 k , the energy contained in the reconstruction to rank-r is k giving the partial energy contribution of the k th-mode. Thus, the Odd columns: spatial modes u k (r i ); even columns: temporal modes v k (t j ), for k=1,2,3,4,20,30,40,46. One unit in the time axis equals 12 s. reconstruction error |δT (r) ti | 2 contains the remaining energy in the leftover ranks E − E(r). In particular, for high r, δT (r) ti is associated with the energy in the small scales. Given that the singular values, satisfying σ k ≤ σ k+1 , decay very fast in the first few modes, as seen in the POD spectra of Figure 9(a) for the three cases analyzed, most of the energy is contained in the low-r reconstructions. To see this, Figure 9(b) shows the energy contribution percentage, E(r)/E, as a function of the reconstruction rank-r for the three cases under study. It is clear that the first mode already contains about two thirds of the energy and successive modes contribute less as the rank rises. The case with the flare has the energy distributed over more modes, reflecting the effect of the large scale flows produced by the flare.
To quantify the degree of the temperature intermittency, we construct the probability distribution function (PDF) by dividing the range of values of T − T ti and this can be done for different percentages of the energy contribution, i.e. taking the corresponding rank r from Fig. 9(b). We chose to use two energy fractions to compare them, namely a residual energy of 10% and 25% (i.e. 1 − E(r)/E). For each case, there is one PDF for each time t i having its own height and width, but we expect them to have similar statistical features. Thus, as it is customary, we rescale the PDF as P → σP and x → x/σ where σ is the standard deviation. Since σ changes with the different time frames, it is more convenient to use the average σ t over all times in order to have a single normalizing parameter. Thus, for every time the PDFs were rescaled with the average standard deviation as P → σ t P and δT r ti → δT r ti / σ t . It is found that the resulting rescaled PDFs are all quite similar having a small overall dispersion. Then we took the average over the time sequence of PDFs (47 frames) in order to have a single PDF that represents the statistical properties of the reconstruction error maps for each energy content. The results for the two energy percentages are shown in figure 10. It is observed that the PDFs of the pre-flare region are significantly broader than the other two, indicating the presence of very large temperature fluctuations.
In order to determine if the PDFs have non standard features like long tails, it is useful to fit some known function to the data. Normal statistics produce a Gaussian PDF, so we can fit a stretched exponential P (x) ∼ exp(−β|x| µ ) to the average PDFs in order to asses how far they are from a Gaussian (µ = 2). In Figure 11 we show the best fits to each of the six cases analyzed. For each PDF we present a fit of the central part and another fit of the tails (with smaller µ). Clear differences are observed among the Flare, Pre-flare and QS cases. However, for the different energy contents, the results do not show great variability. The interesting result is that pre-flare PDFs, in addition to being the broadest of the three (see Fig. 10), have the lowest values of µ and so they are the ones that depart the most from a Gaussian; the long tails indicate there is intermittency, i.e. a relatively large number of intense temperature fluctuations. On the other hand, the PDFs for the Flare have larger µ and thus are closer to a Gaussian; large fluctuation events are then less frequent. Our results agree with the study in [9] that observed lower levels of intermittency in the presence of sources. The work in [9] followed the transport of impurities in a plasma with decaying turbulence (without sources to sustain it) finding evidence of intermittency, but it disappears when there are sources that maintain the turbulence. In our case, the flare plays the role of a source that suppresses intermittency. As expected of a region with no solar activity, in the quiet sun case, the PDF does not exhibit long tails, i.e., the probability of large temperature fluctuations is small.

Heat Flux Estimate
In this section we focus on the region affected by the flare and calculate the thermal flux entering there, driven by the released flare energy. This is made using the time evolution of the original temperature maps and comparing the results with those obtained with the Topos-Chronos method. Next, the heat flux is used to explore the properties of heat transport in the corona.
Since the thermal front produced by the flare is a coherent structure it is expected that the main mechanism bringing energy into the region of interest is advective. Using this premise, we use the computed heat flux to estimate the velocity of the incoming plasma flow, as a function of the coordinate x along the flow (the velocity profile V (x)). However, there might be a diffusive component playing part in the process. To investigate this, we compare the heat flux profile with the temperature and temperature gradient profiles to look for convection or diffusion relations. Assuming that both processes coexist, we estimate velocity profiles taking different levels of diffusion. We also explore the scenario of diffusion being the main transport process by computing the profiles of the diffusivity D = D(x), but we find that it is not viable since D becomes negative.
For the analysis of the heat flux we make the following assumptions: (1) there are no sources of energy inside the region of the temperature maps since it was chosen to be slightly away from the flare site; (2) the main direction of motion is along the x-axis (y displacements are negligible and no information of the dynamics in the z-axis is known: we have an integrated view along z); (3) the internal energy of the plasma is u = 3 2 nk B T and no mechanical work is done by the plasma; (4) we use a constant electron density with the value for the lower solar Corona: n = 10 15 m −3 in the region of interest. Total Heat Flux using Original T-maps: The starting point is the heat transport equation in the absence of sources: Since the thermal front moves almost exclusively in the x direction, we assume that there is no significant variation in the y direction and therefore we take the average of eq. (9) over y. Moreover, we are interested in the energy transport over the whole time interval, τ , and therefore we take the time average of eq.(9) where f y = 1 L´L 0 f dy is the average value in the y-direction and f = 1 τ´τ 0 f dt is the time average. Integrating eq.(10) along the x-coordinate and using that ∂q y /∂y y = 0 we get: Assuming that the first term on the right hand side is zero, meaning there is no heat leaving the region at x = L, the average heat flux profile is where ∆ τ ξ ≡ (ξ(t = τ ) − ξ(t = 0))/τ . In Figure 13 we show the heat flux profile computed from the original temperature maps using (11). The heat flux can have an advective and a diffusive component Q = W v−D∇W (we do not consider here radiative transport) with W = nk B T the thermal energy density and D heat diffusivity. The relative importance of these contributions can be assessed by computing the profiles for the temperature and the temperature gradient and comparing them with the heat flux profile. The temperature profile is taken by averaging the original temperature maps in t and y and similarly for the temperature gradient. These are also shown in Fig. 13 which indicates that the T (x) profile is similar to the Q(x) profile suggesting that the advective component should be dominant one. In contrast, the negative temperature gradient has no clear correlation with Q(x) which seems to indicate that diffusion is not the main drive for the heat flux. From these observations we can compute the advection velocity as In Figure 14 the velocity profile computed in this way is shown for the case with no diffusion (D = 0) and when a small constant diffusivity is assumed to be present (D 0 = (1 and 8) × 10 11 m 2 s −1 ). Diffusivities larger than 8 × 10 11 m 2 s −1 will give sign changes in v which is unlikely. We could also explore the possibility of diffusive transport by computing the diffusivity, D = (v −Q/nk B T )T /∇T . In Figure 15 we show the diffusivity for v = 0 (no advection) and for constant values of the advection velocity. It is seen that D is negative in some regions which is not acceptable. This would disprove the dominance of diffusive transport. There can only be a small contribution superimposed on the advective transport.
Heat Flux using Topos-Chronos: In addition to the decomposition of the temperature fluctuations in optimal modes, the Topos-Chronos method can be used to do a multi-scale analysis of transport. The separation of space and time allows the extraction of prominent spatial structures at different times, ordered by rank. Similarly, temporal representation gives information about the time evolution of the structures at the corresponding rank. This separation is amenable to compute the dominant spatial temperature gradients, that drive the heat fluxes, and the time variation of the thermal energy, independently.
Substituting the representation for the temperature maps in eq. 4 into eq.(11) we obtain the Topos-Chronos decomposition of the heat flux profile: The time derivative as well as the time average affect only the temporal modes; the integral and the average in y affect only the spatial modes. When the space components u k (r i ) are transformed back to a 2D array using the coordinate transformation r i → (x l , y m ), we obtain the matrixû k (x l , y m ). Each mode of the heat flux is thus given by The heat flux profiles for the first four modes are shown in Figure 16. As expected, the dominant mode k = 1 has a profile similar to that of the raw temperature maps in Fig. 13. For the other ranks the flux is smaller but, interestingly, for k = 2, Q is negative, indicating there is a second order return heat flow. A similar analysis can be made for the k = 1 mode. The result is presented in Figure 17 where, as before, a correlation between Q(x) and T (x) is observed, Figure 14: Velocity profiles assuming different levels of constant diffusivity.
but not between Q(x) and the temperature gradient. Thus, advection is again the dominant mechanism and the advection velocity can be computed in the same way. The resulting velocity profile is shown in Figure 18 together with the velocity obtained from the original maps for comparison. The Topos-Chronos velocity is larger in most of the region because it lacks the contributions from higher modes, especially the negative k = 2 rank. This shows that keeping just the lowest mode can give incomplete information. Higher modes could also provide information about the small scale diffusive transport. However, an analysis like the one presented in Fig. 15 does not allow computing the diffusivity since, as before, it turns out to be negative at some points. This is because taking the average over the y direction is not a good procedure since the contribution of small scales is washed out and this is important for higher modes.

Discussion and Conclusions
Temperature maps obtained from SDO/AIA EUV images corresponding to a solar event near an AR in the Solar Corona were analyzed using POD methods and were compared with similar analyses of other maps, one of the same region but before the flare (pre-flare) and another of a quiet region during the same time period (quiet sun). The POD method separates time and space information (Topos-Chronos) thus allowing to determine the dominant space-time scales. The high-rank modes in the decomposition correspond to smaller spatial scales and in some degree with small time-scales. An interesting finding is that there is a rough correlation between Fourier time and spatial scales when the flare has occurred but not for pre-flare or QS states. This may indicate that flare-driven large scale heat flows tend to decay into smaller scales, in a cascade-like process. The Topos-Chronos method was also used to study the statistical properties of the temperature fluctuations. In particular, we reconstructed the temperature maps up to a certain rank, associated with a given energy content. The reconstruction error was thus associated with the small-scale temperature fluctuations. The probability distribution functions (PDFs) of these small-scale fluctuations were obtained for all times and then averaged in time for each of the three regions analyzed. One of the main results of the present paper is the observation that the PDF for the pre-flare state is significantly broader and has longer tails than the PDFs of the flare and quiet sun cases. The pre-flare activity seems to produce larger low-amplitude temperature fluctuations, characteristic of intermittency, which might herald the occurrence of the flare.
A multi-scale analysis of the heat flux was also performed for the region associated with the flare. The thermal flux profiles along the main (x) direction of the flow were computed using the original temperature maps and compared with the temperature variation along x, allowing to obtain the advection velocity profile. Diffusive transport is found to be not relevant. The same computations were applied for the Topos-Chronos decomposition finding also a prevalence of advection with a v(x) profile for the lowest mode that agrees with the one for the original maps. Higher modes are not so relevant for the advective flux but can provide information about small scale diffusion.