© ESO 2021

1. Introduction

The solar corona hosts many different eruptive and explosive events, which are manifests in the various flow patterns found in complex magnetic loop systems. In the past few decades, high temporal and spatial resolution imaging telescopes (Solar Dynamics Observatory (SDO), Solar Terrestrial Relations Observatory (STEREO), Solar and Heliospheric Observatory (SOHO), etc.) have made it possible to observe the intensity changes along the magnetic loops with characteristic spatial scales of at least tens of resolution pixels (for instance, SDO has a resolution of $1\stackrel{″}{.}2$ and SOHO of 5″). This has led to significant efforts to study the physical properties of the intermediate- and large-scale plasma flows at the scale of the visible size of the magnetic structure of interest (e.g., active regions, filaments, coronal bright points) induced during solar flares, including micro- or nano-flares. The main purpose of these efforts is to understand both the direct and indirect effects of energy releases into the ambient solar atmosphere (for an example of such investigation and a comprehensive outline of related literature, see Kamio et al. 2011). These flow phenomena help to establish the cause-effect relationship between properties of the magnetic structures and the plasma flows, which in turn determine the related physical conditions of the solar atmosphere (see Srivastava et al. 2019 and Reale 2014).

However, the observations of smaller scale flows are rather puzzling as most of these flows are of the order of a few times the spatial resolution limit, provided by the instruments on the abovementioned space missions. Recently, Williams et al. (2020) reported the discovery of a large population of so-called coronal strands in active regions made with the High-Resolution Coronal Imager (Hi-C 2.1, with $0\stackrel{″}{.}129$ pixel⁻¹ resolution, see Rachmeler et al. 2019), which is the advanced version of the original Hi-C NASA sounding rocket experiment performed in 2012 (Kobayashi et al. 2014). The existence of these strands reveals that large-scale magnetic objects like active regions have fine small-scale structures. In order to understand the energetics of the flows in the context of the overall energy budget of an active region, it is important to study the nature of the motion of the intensity blobs and to classify them in terms of their characteristic dynamic parameters. For instance, it is important to know which of them are engaged in oscillatory or aperiodic flow-type motions. An example case study of such oscillations is discussed in Philishvili et al. (2017), where oscillatory flows in the active region AR11429 are interpreted as standing longitudinal (acoustic) modes with periods within the range of 116−147 s. The nature of these oscillations drastically differs from those with long periods (> 1 h) reported by Dumbadze et al. (2017) and references therein. As to the aperiodic motions of the observed intensity blobs in the same active region (AR11429), these are the main objective of the present study. In particular, we aim to formulate a mathematically rigorous identification and statistical classification method using SDO/AIA 171 Å for flows in the small-scale strands of similar spatial scales as those observed with the higher resolution instruments (Williams et al. 2020). Our proposed method provides a set of identified aperiodic flow trajectories together with their estimated velocity bounds (taking into account corrections in pixel distances due to the projection deformation of the solar spherical surface in flat CCD images) and classification as either subsonic, sonic, or supersonic. We will do this by estimating the local background sound speeds by means of the differential emission measure (DEM) inversion algorithm of Plowman et al. (2013).

The investigation of flows in the solar corona with a variety of spatiotemporal scales is an important task as they carry information about their acceleration sources, which are related to the thermal and mechanical energy release during flows and sporadic energy release events as mentioned above. Dynamic mass flows in coronal magnetic loops are a typical attribute of a rapid energy deposition in the atmosphere. These flows affect the coronal heating rate by evaporating (Acton et al. 1982; Hirayama 1974; Neupert 1968) and condensing chromospheric plasma (Canfield et al. 1990; Fisher 1987) into the corona. Li et al. (2015) observed the chromospheric evaporation process in flaring loops and Nitta et al. (2012) found similar upward flow motions with 100−500 km s⁻¹ velocities associated with small transient brightening. These events can be separated into two different types (Li et al. 2015): subsonic flows (gentle evaporation), which are characterized by tens of kilometers per second velocities, and supersonic flows (explosive evaporation) which are characterized by velocities of hundreds of kilometers per second. The rapid acceleration of the plasma can be energetically supported by thermal (Shergelashvili et al. 2007) and magnetohydrodynamic (Shergelashvili et al. 2006) instability processes. In coronal hole bright points (Bagashvili et al. 2018) (representing small-scale, active-region-like structures) and in active regions (Joshi et al. 2020), quasi-oscillatory variability has been discovered to precede jet outflows, which can be interpreted as indications of the presence of instability processes. Furthermore, upflows are observed at the edges of active regions (Boutry et al. 2012). All these types of ‘intermediate scale’ flow patterns are proposed to contribute to the acceleration of the large-scale solar wind in the inner heliosphere (see Brooks & Warren 2010). It has been shown that the existence of magnetic reconnection sources in the corona can lead to extremely rapid acceleration, even forming discontinuous patterns of the solar wind (Shergelashvili et al. 2020). Shergelashvili & Fichtner (2012) and Vainio et al. (2003) discuss a causal chain in which the energy of intermediate-scale flows can be dissipated into heat and wave turbulence energy at the source region where the fast and slow solar wind patterns are formed. While the latter determine the space weather conditions in the interplanetary space and shape the entire heliosphere, it is plausible that the small-scale flows are also part of this chain. Small-scale flows are abundant in active regions and other magnetic structures, and represent carriers of energy for the formation of the multi-scale structures of magnetic objects, their distortions, and instabilities in the solar corona. Therefore, a thorough understanding of the dynamical properties of those small-scale flows, with the future perspective to upgrade towards full-scale statistical, machine, or deep learning models, is a crucially important problem.

In this paper, we propose a method to locate and track small-scale flow patterns in order to estimate their dynamic properties. The flow tracking algorithm is outlined in Sect. 2.1, the estimation of the local background temperature and sound speed is discussed in Sect. 2.2, and in Sect. 2.3 we perform simultaneous statistical inference on the flow velocities in order to classify them as subsonic, sonic, or supersonic. The results and corresponding discussion are presented in Sect. 3, and Sect. 4 consists of conclusive remarks.

2. Methodology

2.1. Plasma flow identification

In our analysis we consider the active region (AR) 11429, which was visible on the solar surface from 4 March 2012 to 15 March 2012. We used the high-resolution extreme ultraviolet (EUV) monitor images taken by AIA onboard the SDO. In order to locate and track the flow patterns, we considered the image data from the 171 Å spectral line. The AIA images, which have a 12 second time cadence, were analyzed during two intervals of time: (i) on 6 March 2012 from 12:22:00UT to 12:35:00UT when the active region was situated at −40 degrees from the central meridian of the visible solar disk (the limb data); and (ii) on 9 March 2012 from 03:00:00UT to 03:24:00UT when the active region was in the vicinity of the central meridian (the central data). On all AIA data we applied standard image preparation and cleaning techniques using the Solarsoftware (SSW) package.

It is an arduous task to locate and track plasma flows in real time by observing an 171 Å spectral line consisting of thousands of pixel intensities that change over time. We propose a method that automatically locates and tracks any amount of candidate plasma flows simultaneously over a spectral line, and yields the coordinates and times of each flow. This methodology is based on two properties of the plasma flows: (i) A local signature of a plasma blob must appear (as an enhanced intensity pattern) and (ii) it must move. We define I_k(x_i, y_i) as the intensity at the coordinates (x_i, y_i) at time t_k. Since movement is indicated by a change in intensity at the specific coordinates (x_i, y_i), we define the change in intensity as ΔI_k(x, y) = I_k(x, y) − I_k − 1(x, y). The two properties of plasma flow mentioned above indicate that a combination of movement ΔI_k and intensity I_k is considered when locating and tracking a flow. Therefore, we consider the linear combination γI_k(x, y)+(1 − γ)ΔI_k(x, y) for some 0 ≤ γ ≤ 1 as the quantity used to locate and track plasma flows. For time steps k = 2, 3, 4, … we perform the following iterative procedures:

Step 1. Introducing new intensity enhancements and tracking candidate flows.

We consider the M points (pixels) that have the largest values for γI_k(x, y)+(1 − γ)ΔI_k(x, y) and denote these as (x_1, k, y_1, k),(x_2, k, y_2, k)…, (x_M, k, y_M, k). In order to introduce the candidate intensity enhancements, we execute the following iterative procedure for i = 1, 2, …, M:

Let ${\mathbb{S}}_{i,k}=\{(x,\mathit{y}):\sqrt{{(x-x_{i,k})}^2+{(\mathit{y}-{\mathit{y}}_{i,k})}^2}\le 3\}$ and calculate

$$\begin{array}{c}\hfill x_{i,k}^{\ast }=\frac{\sum _{(x,y)\in {\mathbb{S}}_{i,k}}x{e}^{\lambda (I_k(x,y)-\gamma I_{k-1}(x,y))}}{\sum _{(x,y)\in {\mathbb{S}}_k}{e}^{\lambda (I_k(x,y)-\gamma I_{k-1}(x,y))}}\end{array}$$(1)

$$\begin{array}{c}\hfill y_{i,k}^{\ast }=\frac{\sum _{(x,y)\in {\mathbb{S}}_k}y{e}^{\lambda (I_k(x,y)-\gamma I_{k-1}(x,y))}}{\sum _{(x,y)\in {\mathbb{S}}_k}{e}^{\lambda (I_k(x,y)-\gamma I_{k-1}(x,y))}}\end{array}$$(2)

If $\sqrt{{(x_{i,k}^{\ast }-x_{i,k})}^2+{({\mathit{y}}_{i,k}^{\ast }-{\mathit{y}}_{i,k})}^2}>0.01$ then set $(x_{i,k},y_{i,k})=(x_{i,k}^{*},y_{i,k}^{*})$ and recalculate ${\mathbb{S}}_{i,k},x_{i,k}^{*}$ and $y_{i,k}^{*}$.

Otherwise if $\sqrt{{(x_{i,k}^{\ast }-x_{i,k})}^2+{({\mathit{y}}_{i,k}^{\ast }-{\mathit{y}}_{i,k})}^2}\le 0.01,$ then the centroid has converged.

Further, to enable the tracking of the advance of these introduced enhancements, we perform the same iterative procedure as above for i = M + 1, M + 2, …, M + N_k − 1, where it is implied that (x_{1, k − 1}, y_{1, k − 1}),(x_{2, k − 1}, y_{2, k − 1}),…(x_{Nk − 1, k − 1}, y_{Nk − 1, k − 1}) denote the previous time’s candidate flow points, where N₁ = 0 and (x_{i + M, k}, y_{i + M, k}) = (x_{i, k − 1}, y_{i, k − 1}) for i = 1, 2, …, N_k are set, where N_k is the number of candidate flows.

The hyper-parameters γ and λ directly affect how the method locates and tracks flows. The parameter γ specifies what feature is followed by the method, and the parameter λ specifies how the method follows it. For values of γ close to 0, the method will only follow the largest changes in intensity, and for values close to 1, the method will only follow the largest intensity enhancements. A central compromise of $\gamma =\frac{1}{2}$ is recommended since this balances the requirement for the method to follow large intensities and also follow changes in these intensities. For values of λ close to 0, the formula for (x^*, y^*) will yield a point that is at the local mean of I_k(x, y) − γ I_k − 1(x, y), and for a large value of λ, the formula for (x^*, y^*) will give a point that is at the local maximum of I_k(x, y)−γ I_k − 1(x, y). In general, if it is required that the maximum of I_k(x, y)−γ I_k − 1(x, y) be given a weight ω, then it follows that the parameter λ be chosen as λ_ω, which is the solution to the equation

$$\begin{array}{c}\hfill \frac{{\lambda }_{\omega }\underset{(x,y)\in {\mathbb{S}}_{i,k}}{max}\{I_k(x,y)-\gamma I_{k-1}(x,y)\}}{log\sum _{(x,y)\in {\mathbb{S}}_{i,k}}{e}^{{\lambda }_{\omega }(I_k(x,y)-\gamma I_{k-1}(x,y))}}=log\omega .\end{array}$$(3)

It is recommended that $\lambda ={\lambda }_{\frac{1}{2}}$ is chosen so that half of the weight is given to the local maximum of I_k(x, y)−γ I_k − 1(x, y). This, therefore, provides a compromise between a mean and maximum method of following a flow. The result of this iterative procedure, given by Eq. (2), is illustrated in Fig. 1. While the resulting movement of a candidate flow over time steps k = 2, 3, 4, … is illustrated in Fig. 2.

Fig. 1. Illustration of a centroid with a radius of three pixels converging over iterations 1 (left), 6 (center), and 12 (right). The blue point indicates the start point of the centroid, and the green point indicates the current centroid estimate.

Fig. 2. Illustration of a centroid moving from time step 15 (left) to time step 20 (center) to time step 25 (right). The blue point indicates the start point of the centroid, and the green point indicates the current centroid estimate.

Step 2. The candidate flow assessment.

The change in direction and change in distance of each of the candidate flows are assessed in order to determine whether they constitute a flow or not. We denote the direction (θ_i, k) and velocity (v_i, k) of each candidate flow i = 1, 2, …M + N_k − 1 as

$$\begin{array}{c}\hfill {\theta }_{i,k}=arctan\frac{y_{i,k}-y_{i,k-1}}{x_{i,k}-x_{i,k-1}}\end{array}$$(4)

$$\begin{array}{c}\hfill v_{i,t}=\frac{\sqrt{{(x_{i,k}-x_{i,k-1})}^2+{(y_{i,k}-y_{i,k-1})}^2}}{t_{i,k}-t_{i,k-1}}.\end{array}$$(5)

Consequently, a candidate flow is disregarded and removed if |Δθ_i, k|> Θ or v_i, k ∉ [v_min, v_max].

Altogether this method has four dummy parameters (M, v_min, v_max, and Θ) and two hyper-parameters (γ and λ) that are set by the user and they output a set of flow position coordinates {(x_i, k,y_i, k)}_{k = ki, ki + 1, …, ki + Ki} and occurrence times {t_i, k}_{k = ki, ki + 1, …, ki + Ki} for i = 1, 2, …, N, where N is the total number of identified flows, K_i is the number of time steps of the ith flow and k_i is the start time step of the ith flow. Ideally, M would be chosen to be as large as possible so that the greatest amount of candidate flows can be identified. The velocity parameters (v_min and v_max) control the minimum and maximum velocity of the candidate flows, and should be set to identify flows with reasonable velocities. Alternatively, if the user is interested in flows of a specific velocity, then these parameters can be initialized accordingly. The direction parameter Θ controls the maximum change in direction that a candidate flow is allowed to have at each observed time. The ranges of the minimum and maximum velocity and the angle Θ approve the candidate flows and filter out all others that do not satisfy the required conditions.

The proposed methodology is applied to the SDO image data from the 171 Å spectral line for the limb data and again for the central data. The dummy parameters are chosen as M = 100, v_min = 40 km s⁻¹, v_max = 200 km s⁻¹, and Θ = π/6. The hyper-parameters are set as γ = 0.5 and $\lambda ={\lambda }_{\frac{1}{2}}$ as we proposed above. The application of the methodology yielded a total of 5 flows for the limb data and a further 13 flows for the central data. These results are illustrated in Fig. 3 and are summarized in Table 1.

Fig. 3. Illustration of the flows identified by applying the proposed methodology to the SDO image data from the 171 Å spectral line for the date 6 March 2012 from 12:22:00 to 12:35:00 referred to as the limb data (top) and again for the date 9 March 2012 from 03:00:00 to 03:24:00 referred to as the central data (bottom).

The method can easily be extended to locate and track flows over multiple spectral lines simultaneously by substituting λ(I_k(x,y)−γ I_k − 1(x,y)) with ${\sum }_f{\lambda }_f\left(I_k^f(x,y)-\gamma I_{k-1}^f(x,y)\right)$ in the iterative procedure, where f indicates the spectral lines of interest. This would require the specification of an additional parameter value λ_f for each additional spectral line f.

2.2. DEM analysis and background sound speed estimation

We estimate the mean background temperature for each flow by applying the DEM inversion method of Hannah & Kontar (2012) to the six AIA channels, fitting Gaussian functions to the resulting DEM curves and taking the mean of the peaks of the fitted Gaussian functions. Following the algorithm developed by Plowman et al. (2013), for each flow we take the vector of pixels (the set of pixels are selected by the iterative procedure described above) along its trajectory at the time preceding the appearance of the intensity enhancement (blob) and extract the corresponding six AIA channel values along it. These six AIA channels are 94 Å, 171 Å, 193 Å, 211 Å, 304 Å and 335 Å. We then apply the DEM inversion method to obtain reconstructed DEM distributions for each flow, respectively. Gaussian functions are fitted to the DEM curves to estimate the peak temperature for each pixel at each time step. We then take the mean of these peak values of the temperature for the entire time span from the start of the observation until the occurrence of the flow to estimate the local background temperature. The local background sound speed of each flow can be estimated by using the formula

$$\begin{array}{c}\hfill Cs=\sqrt{\frac{\gamma k_{\mathrm{B}}T}{\stackrel{\sim }{\mu }{m}_p}},\end{array}$$(6)

where μ = 0.6 is the mean atomic mass factor of hydrogen plasma, γ = 5/3, k_B is the Bolzman constant, and T is the estimated background temperature. The error of the estimated background temperature T can be estimated by the shape parameter σ_T of the Gaussian functions fitted to the DEM curves (Plowman et al. 2013). Equation (6) is again used to estimate the error of the resulting background sound speed, denoted by ${\widehat{\sigma }}_{C_{\mathrm{s}}}$. The estimated background sound speeds and their corresponding errors are summarized in Table 1.

2.3. Flow velocity inference

In this section, we describe the statistical inference of the flow velocity detection and their classification. However, before doing this, it is reasonable to note that the 2D AIA images comply with a deformation of the spherical surface image of the Sun due to the projection effect. In order to take into account this deformation and recalculate the distances covered by AIA pixels, we use the approach that each pixel projection on the solar surface along the line of sight represents a locally planar non-rectangular section of the solar surface. If one then performs straightforward transformations of coordinates, the solution of this geometrical ansatz (derivation of it is rather cumbersome but straightforward, therefore for the sake of brevity we just demonstrate here the final result), it can be shown that the mentioned non-rectangular area would be seen as a parallelogram if it were situated right in the center of the visible solar disk image with the diagonals

$$\begin{array}{c}\hfill M_1=\frac{L}{cos(\theta )}\sqrt{2+{(tan(\theta )+\frac{tan(\varphi )}{cos(\theta )})}^2}\end{array}$$(7)

$$\begin{array}{c}\hfill M_2=\frac{L}{cos(\theta )}\sqrt{2+{(tan(\theta )-\frac{tan(\varphi )}{cos(\theta )})}^2}.\end{array}$$(8)

The angle between them is defined as

$$\begin{array}{c}\hfill cos(\alpha )=\frac{M_1^2+M_2^2-4L^2}{2M_1{M}_2},\end{array}$$(9)

where, L = R₀/R is the standard length of the pixel side, which is determined by the ratio of solar radius R₀ in kilometers and the actual solar radius R measured in pixels seen by the observer at the current moment in time. The angles ϕ and θ represent, respectively, the angular longitudinal and latitudinal distance of the pixel from the center of the solar disk, also with respect to the observer’s position, which are defined as

$$\begin{array}{c}\hfill sin(\varphi )=\frac{x-CRPIX1}{R},sin(\theta )=\frac{y-CRPIX2}{R}.\end{array}$$(10)

Here, x and y are pixel coordinates as defined above and CRPIX1 and CRPIX2 are the corresponding pixel coordinates of the visible solar center on the AIA_PREP-ed SDO image. From these definitions, the sides of the obtained parallelogram can be defined as

$$\begin{array}{c}\hfill L_1^{\prime }=\sqrt{\frac{M_1^2+M_2^2}{4}-\frac{M_1{M}_2}{2}cos(\alpha )}\end{array}$$(11)

$$\begin{array}{c}\hfill L_2^{\prime }=\sqrt{\frac{M_1^2+M_2^2}{4}+\frac{M_1{M}_2}{2}cos(\alpha )}.\end{array}$$(12)

We take the smallest side of the parallelogram $L_{\text{min}}=L_1^{\prime }$ and the largest diagonal L_max = max{M₁, M₂} as the minimum and maximum length scales, which are used for the estimation of the lower and upper bounds of the measured velocities. In this study, we do not perform any 3D reconstruction of the magnetic field. Therefore, another projection effect occurs related to the angle between the local orientation of the loop sustaining the observed flow and the line of sight. Although, the additional effect caused by the orientation of the magnetic field only can lead to the conclusion that the actual velocities can never be smaller (but they are equal or somewhat larger) than either the lower and upper bounds of them. The latter are estimated throughout the current work.

The above specifications imply that the observed velocities v₁, v₂, …, v_k are bounded by minima $v_1^{\text{min}},v_2^{\text{min}},\dots ,v_k^{\text{min}}$ and maxima $v_1^{\text{max}},v_2^{\text{max}},\dots ,v_k^{\text{max}}$ where

$$\begin{array}{c}\hfill v_j^{\min }=\frac{L_{\min }\sqrt{{(x_j-x_{j-1})}^2+{(y_j-y_{j-1})}^2}}{t_j-t_{j-1}}\end{array}$$(13)

$$\begin{array}{c}\hfill v_j^{\max }=\frac{L_{\max }\sqrt{{(x_j-x_{j-1})}^2+{(y_j-y_{j-1})}^2}}{t_j-t_{j-1}},\end{array}$$(14)

where j = 1, 2, ⋯, k.

In order to test whether a the apparent velocity of a flow is sonic, subsonic, or supersonic at time t, we perform the following hypothesis tests simultaneously:

$$\begin{array}{c}\hfill H_{01}:v^{\max }(t)=C_{\mathrm{s}}{H}_{11}:v^{\max }(t)<C_{\mathrm{s}}{H}_{21}:v^{\max }(t)>C_{\mathrm{s}}\end{array}$$(15)

$$\begin{array}{c}\hfill H_{02}:v^{\min }(t)=C_{\mathrm{s}}{H}_{12}:v^{\min }(t)>C_{\mathrm{s}}{H}_{22}:v^{\min }(t)<C_{\mathrm{s}},\end{array}$$(16)

where v^max(t) and v^min(t) are the upper and lower flow velocity bounds at time t and C_s is its background sound speed. The flow velocity upper and lower bounds v^max and v^min are approximated by a quadratic regression model function of time,

$$\begin{array}{c}\hfill v^{\max }(t)={\beta }_0^{\max }+{\beta }_1^{\max }t+{\beta }_2^{\max }{t}^2+{ϵ}^{\max },\end{array}$$(17)

$$\begin{array}{c}\hfill v^{\min }(t)={\beta }_0^{\min }+{\beta }_1^{\min }t+{\beta }_2^{\min }{t}^2+{ϵ}^{\min },\end{array}$$(18)

where ${\beta }_0^{\text{max}}$, ${\beta }_1^{\text{max}}$, ${\beta }_2^{\text{max}}$, ${\beta }_0^{\text{min}}$, ${\beta }_1^{\text{min}}$, and ${\beta }_2^{\text{min}}$ are real-valued coefficients and ϵ^max and ϵ^min are the error terms that are assumed to be normally distributed with mean 0 and standard deviation σ_ϵ^max and σ_ϵ^min respectively. The least squares estimators of v^max(t) and v^min(t) are given by (Hastie et al. 2001)

$$\begin{array}{c}\hfill {\widehat{v}}^{\max }(t)={\left[\begin{array}{c}1\\ t\\ t^2\end{array}\right]}^T{\left({\left[\begin{array}{c}\underline{1}\\ \underline{t}\\ {\underline{t}}^2\end{array}\right]}^T\left[\begin{array}{c}\underline{1}\\ \underline{t}\\ {\underline{t}}^2\end{array}\right]\right)}^{-1}{\left[\begin{array}{c}\underline{1}\\ \underline{t}\\ {\underline{t}}^2\end{array}\right]}^T{\underline{v}}^{\max },\end{array}$$(19)

$$\begin{array}{c}\hfill {\widehat{v}}^{\min }(t)={\left[\begin{array}{c}1\\ t\\ t^2\end{array}\right]}^T{\left({\left[\begin{array}{c}\underline{1}\\ \underline{t}\\ {\underline{t}}^2\end{array}\right]}^T\left[\begin{array}{c}\underline{1}\\ \underline{t}\\ {\underline{t}}^2\end{array}\right]\right)}^{-1}{\left[\begin{array}{c}\underline{1}\\ \underline{t}\\ {\underline{t}}^2\end{array}\right]}^T{\underline{v}}^{\min },\end{array}$$(20)

where $\underline{t}={[t_1,t_2,\cdots ,t_K]}^T$, ${\underline{t}}^2={[t_1^2,t_2^2,\cdots ,t_K^2]}^T$, ${\underline{\mathit{v}}}^{\max }={[{\mathit{v}}_1^{\max },{\mathit{v}}_2^{\max },\cdots {\mathit{v}}_K^{\max }]}^T$ and ${\underline{\mathit{v}}}^{\min }={[{\mathit{v}}_1^{\min },{\mathit{v}}_2^{\min },\cdots {\mathit{v}}_K^{\min }]}^T$ are the vectors of observational moments of time, the squares of them, and the corresponding observed values of the velocity upper and lower bounds respectively.

Under H₀₁ we have that ${\widehat{\mathit{v}}}^{\max }(t)-C_{\mathrm{s}}\sim N(0,{\sigma }_{C_{\mathrm{s}}}^2+{\sigma }_{{ϵ}^{\max }}^2\xi (t))$ and under H₀₂ we have that ${\widehat{\mathit{v}}}^{\min }(t)-C_{\mathrm{s}}\sim N(0,{\sigma }_{C_{\mathrm{s}}}^2+{\sigma }_{{ϵ}^{\min }}^2\xi (t))$, where in both cases

$$\begin{array}{c}\hfill \xi (t)=\left[\begin{array}{c}1\\ t\\ t^2\end{array}\right]{\left({\left[\begin{array}{c}\underline{1}\\ \underline{t}\\ {\underline{t}}^2\end{array}\right]}^T\left[\begin{array}{c}\underline{1}\\ \underline{t}\\ {\underline{t}}^2\end{array}\right]\right)}^{-1}{\left[\begin{array}{c}1\\ t\\ t^2\end{array}\right]}^T.\end{array}$$(21)

It follows that the test statistics under H₀₁ is given by

$$\begin{array}{c}\hfill T^{\max }(t)=\frac{{\widehat{v}}^{\max }(t)-C_{\mathrm{s}}}{\sqrt{{\widehat{\sigma }}_{C_{\mathrm{s}}}^2+{\widehat{\sigma }}_{{ϵ}^{\max }}^2\xi (t)}}\sim t_{K-1},\end{array}$$(22)

and the test statistics under H₀₂ is given by

$$\begin{array}{c}\hfill T^{\min }(t)=\frac{{\widehat{v}}^{\min }(t)-C_{\mathrm{s}}}{\sqrt{{\widehat{\sigma }}_{C_{\mathrm{s}}}^2+{\widehat{\sigma }}_{{ϵ}^{\min }}^2\xi (t)}}\sim t_{K-1},\end{array}$$(23)

where t_K − 1 denotes the Student-t distribution with K − 1 degrees of freedom and ${\widehat{\sigma }}_{{ϵ}^{\max }}$ and ${\widehat{\sigma }}_{{ϵ}^{\min }}$ denote the standard errors of ${\widehat{\mathit{v}}}^{\max }(t)$ and ${\widehat{\mathit{v}}}^{\min }(t)$ respectively, which are given by

$$\begin{array}{c}\hfill {\widehat{\sigma }}_{{ϵ}^{\max }}=\sqrt{\frac{1}{K}{\displaystyle \sum _{k=1}^K}{({\widehat{v}}^{\max }(t_k)-v_k^{\max })}^2}\end{array}$$(24)

$$\begin{array}{c}\hfill {\widehat{\sigma }}_{{ϵ}^{\min }}=\sqrt{\frac{1}{K}{\displaystyle \sum _{k=1}^K}{({\widehat{v}}^{\min }(t_k)-v_k^{\min })}^2}.\end{array}$$(25)

At time t with a 100(1 − α)% level of significance, H₀₂ is rejected and H₁₂ is accepted when T^min(t) > t_K − 1(1 − α)), H₀₁ is rejected and H₁₁ is accepted when T^max(t) < t_K − 1(α), H₀₁ is rejected and H₂₁ is accepted when T^min(t) < t_K − 1(α), and H₀₂ is rejected and H₂₂ is accepted when T^max(t) > t_K − 1(1 − α)). The amount t_K − 1(α) denotes the αth quantile of the Student-t distribution with K − 1 degrees of freedom. These hypothesis tests form the basis on which a flow is classified as supersonic, subsonic, or sonic.

2.4. Classification of flows

At time t with a 100(1 − α)% level of significance, a flow is classified as supersonic if H₀₂ is rejected and H₁₂ is accepted, otherwise as subsonic if H₀₁ is rejected and H₁₁ is accepted, otherwise as sonic if H₀₁ is rejected and H₂₁ and H₀₂ is rejected and H₂₂ is accepted. At a 100(1 − α)% level of significance, a flow can have different classifications at different times, and at other times no classification at all. In order to determine the predominant classification of the flow at a 100(1 − α)% level of significance, it is necessary to determine the cumulative durations that the flow was supersonic, sonic, subsonic, or none at this level of significance. Let Π^supersonic(α), Π^sonic(α), and Π^subsonic(α) be the cumulative durations that the flow is classified as supersonic, sonic, and subsonic respectively at a 100(1 − α)% level of significance:

$$\begin{array}{c}\hfill {\mathrm{\Pi }}^{\mathrm{supersonic}}(\alpha )={\displaystyle \int }I\{T^{\min }(t)>t_{K-1}(1-\alpha )\}\mathrm{d}t\end{array}$$(26)

$$\begin{array}{c}\hfill {\mathrm{\Pi }}^{\mathrm{sonic}}(\alpha )={\displaystyle \int }I\{T^{\min }(t)<t_{K-1}(\alpha )\}I\{T^{\max }(t)>t_{K-1}(1-\alpha )\}\mathrm{d}t\end{array}$$(27)

$$\begin{array}{c}\hfill {\mathrm{\Pi }}^{\mathrm{subsonic}}(\alpha )={\displaystyle \int }I\{T^{\max }(t)<t_{K-1}(\alpha )\}\mathrm{d}t,\end{array}$$(28)

where operator I denotes indicator function.

A flow is said to be predominantly supersonic, sonic, or subsonic at a 100(1 − α)% level of significance if, respectively, Π^supersonic(α), Π^sonic(α), or Π^subsonic(α) is the largest. However, if both Π^supersonic(α) and Π^subsonic(α) are of a similar size then it is not clear which classification is predominant. In such a case, the flow is said to be transonic since it transitions from one predominant classification to another over time.

The predominant classification of a flow depends crucially on two properties of the estimated velocity bounds: their closeness to the background sound speed C_s and the size of their respective standard errors ${\widehat{\sigma }}_{{ϵ}^{\max }}$ and ${\widehat{\sigma }}_{{ϵ}^{\min }}$. If ${\widehat{\mathit{v}}}^{\min }(t)-C_{\mathrm{s}}$ is larger then T^min(t) is larger, which means that the flow is classified as supersonic at a higher level of significance at time t. If $C_{\mathrm{s}}-{\widehat{\mathit{v}}}^{\max }(t)$ is larger then −T^max(t) is larger, which means that the flow is classified as subsonic at a higher level of significance at time t. If both $C_{\mathrm{s}}-{\widehat{\mathit{v}}}^{\min }(t)$ and ${\widehat{\mathit{v}}}^{\max }(t)-C_{\mathrm{s}}$ are larger then both −T^min(t) and T^max(t) are larger, which means that the flow is classified as sonic at a higher level of significance at time t. The size of the standard errors ${\widehat{\sigma }}_{{ϵ}^{\max }}$ and ${\widehat{\sigma }}_{{ϵ}^{\min }}$ inversely increase the level of significance of the flow’s classification. If the standard errors are larger, then the flow will be classified at a lower level of significance, and vice versa.

3. Results and discussion

In this paper we propose a three-stage process by which plasma flows are identified, their background temperatures and sound speeds are estimated, and their velocities are classified as either supersonic, subsonic, or sonic. We develop two novel methodologies. The first identifies plasma flows by means of a converging centroid algorithm and the second classifies their velocities as supersonic, subsonic, or sonic by means of statistical inference. This yields a total of eighteen identified plasma flows, 5 of which occurred in the limb data and 13 of which occurred in the central data.

The parameters values v_min = 40 km s⁻¹, v_max = 200 km s⁻¹, and Θ = π/6 ensure that the identified flows move in a consistent direction with velocities between 40 km s⁻¹ and 200 km s⁻¹. These parameter values can be adjusted based on the user’s requirements. The smaller v_min and the larger v_max and Θ are chosen, the larger the number of identified flows, and vice versa. The identified flows have a duration of around one minute and they cover a distance in the range between 4000 and 10 000 kilometers. Out of the 18 identified plasma flows, 11 are classified as predominantly supersonic, 3 as subsonic and 3 as sonic at a 70% level of confidence. Two of these flows have large supersonic and subsonic durations, which implies that these flows form a subclass of transonic flows. A lower number of flows are classified at higher levels of significance, since the confidence in the classifications of these flows is higher. At a 80% level of confidence, 12 flows are classified, 8 as supersonic and 4 as subsonic, and at a 90% level of confidence only 6 flows are classified, 3 as supersonic and 3 as subsonic. It is important to note that 5 of the 6 flows that are classified at a 90% level of significance have low standard errors (see Table 1) with the exception of Flow 11, which has an estimated velocity lower bound that is much larger than the background sound speed (see Fig. 4). This indicates that a low standard error contributes to the confidence by which a flow is classified, and vice versa.

Fig. 4. Velocity lower bound ${\widehat{\mathit{v}}}^{\min }$ and upper bound ${\widehat{\mathit{v}}}^{\max }$ together with the classification level of significance 100(1 − α)%. The red, blue, and green (solid line) curves indicate a supersonic, sonic, and subsonic classification (left ordinate axes), respectively at a 50% level of significance. The dotted line curves show the level of significance (rights ordinate axes). The horizontal dashed line represents the background sound speed. These flows all have a majority supersonic classification.

Fig. 5. Same as in Fig. 4 for subsonic flow patterns.

We utilize the upper and lower bounds of the flow velocities (upper and lower thick solid lines in each panel of Figs. 4–7; the velocity values are shown on the left ordinate axis) together with their corresponding significance parameter (dotted lines in the same panels, the values shown on right ordinate axis) in order to classify a flow as supersonic, subsonic, or sonic during the flow time span at a certain level of significance. If the flow is classified as supersonic for the majority of its duration at a certain level of significance, then the flow is said to be predominantly supersonic at that level of significance, similarly for sonic or subsonic. The flows are clearly identified as supersonic, subsonic, or sonic and are illustrated in Figs. 4–6, respectively. In two specific cases, which we show in Fig. 7, there is no clear predominant classification of the flow at a 70% level of significance as the supersonic and subsonic cumulative durations are both large, although one of them dominates (in Table 2 the corresponding classification indicators are marked with an asterisk ). Consequently these flows transition from subsonic to supersonic states or vice versa. In such a case, the flow is said to be transonic.

Fig. 6. Same as in Fig. 4 for sonic flow patterns.

Fig. 7. Same as in Fig. 4 for transonic flow patterns.

The method of flow identification and classification proposed in this work can be applied to SDO images of active regions in all available EUV spectral lines to reveal the characteristics of the small-scale flow patterns at different altitudes in the solar corona. For 171 Å images it enables us to differentiate flows with regard to their velocities and their relation to the background sound speed at different levels of statistical confidence. In many cases we observe that the flows accelerate or decelerate and in some cases transition from one classification to another. The velocities are estimated taking into account projection due to the sphericity of the solar surface and, consequently, we determine their upper and lower bounds v^max and v^min. We validate the applied distance correction method by matching the results of the velocity estimation close to the limb (flows 1–5) and at the visible central meridian (flows 6–18), which appeared to be comparable to each other. A further specification of the velocity curves can be achieved by evaluation of the vector magnetic field projection related uncertainties addressed in Falconer et al. (2016) and references therein and also reconstruction of the 3D structure of magnetic loops (Tiwari et al. 2017), and taking into account the orientation of different parts of the magnetic loops (presumably playing the role of the flow guides) with regard to the line of sight. This kind of modeling is beyond the scope of this study. However, building on the current results, one could make some intuitive indications. Plausibly, the effect of the magnetic loop curvature and its projection on the line of sight are mainly manifested in reduced apparent velocities compared to the real motion. This means that the method presented in this paper correctly classifies the supersonic flow patterns, while some of the sonic flows found might also be supersonic as their velocities can be underestimated because of the just-mentioned 3D field orientation. On the other hand, the proper analytical and numerical modeling of sources of information and the dynamics of intensity enhancement in magnetic loops, would shed light on the nature of the observed small-scale motions, perhaps attaching to them the role of tracers or ‘fingerprints’ of the small-scale magnetic structure of solar active regions.

The proposed method can be applied to any solar image data and is expected to perform better on higher resolution images such as Hi-C 2.1 data. Furthermore, the converging centroid algorithm of Sect. 2.1 can easily be extended so that flows can be located and tracked on multiple spectral lines simultaneously. This can be done by substituting λ(I_k(x,y)−γ I_k − 1(x,y)) with ${\sum }_f{\lambda }_f\left(I_k^f(x,y)-\gamma I_{k-1}^f(x,y)\right)$, where f indicates the spectral lines of interest. It is important to note that for each additional spectral line f it is necessary to specify an additional parameter value λ_f. This extension of the flow location and tracking methodology, together with further applications to high resolution solar images, such as Hi-C 2.1 data, is the matter of our forthcoming investigations.

4. Conclusions

The number and characteristics of plasma flows occurring on the surface of the Sun has an inherent effect on the physical conditions of the solar atmosphere. The results of the proposed methodology can be used to characterize the physical properties of the solar atmosphere, and in turn be used to analyze the cause-effect relationship between plasma flows and magnetic reconnection and energy release events. The spacial configuration of the flow presumably follows the complex structure of active regions, which recently have been confirmed and reported by high-resolution observations by Williams et al. (2020). It is known that the distribution of characteristic sizes of magnetic strands follows a certain pattern of scale invariance with noticeable slopes of power laws in the logarithmic scale. This manifests the fact that these distributions are the result of strongly non-equilibrium statistical processes (Maes et al. 2009, 2007) of magnetic strand formation and diffusion, which play the role of ‘seeds’ in coronal send-pile models that predict the presence of solar flares as avalanches of small-scale magnetic reconnection processes (see Baiesi et al. 2008, and references therein). This statistical method of flow velocity evaluation in EUV images can be of more wide use. For instance, coronal rains or other motions of the intensity enhancements, such as quiet Sun loops, filament channels (Panesar et al. 2020), other unstable coronal structures (Bagashvili et al. 2018; Joshi et al. 2020), or coronal rains (Vashalomidze et al. 2015; Antolin et al. 2015) can also be traced in a similar manner; only the centroid convergence and the following parameters need to be adjusted accordingly. In the future, such an advanced statistical analysis would be useful for identifying the characteristics of plasma flows as inputs to a machine learning or neural network algorithms, where the result would be a characterisation of an active region and the measurement of its evolution.

Acknowledgments

The work was supported by Shota Rustaveli National Science Foundation grants FR17_609, DI-2016-52. Work of E.P was supported under Shota Rustaveli National Science Foundation grant for young scientists for scientific research internships abroad IG/46/1/16. Work of B.M.S. was supported by the Austrian Fonds zur Forderung der wissenschaftlichen Forschung (FWF) under projects P25640-N27, S11606-N16. S.P. was supported by the projects C14/19/089 (KU Leuven), G.0D07.19N (FWO-Vlaanderen), C 90347 (ESA Prodex). M.L.K. additionally acknowledges support from the projects I2939-N27 and S11606-N16 of the Austrian Science Fund (FWF) and grant number 075-15-2019-1875 from the government of the Russian Federation under the project called “Study of stars with exoplanets”. We are thankful to the anonymous reviewer for the valuable remarks that led to the significant improvement of the manuscript content.

References

All Tables

All Figures

	Fig. 1. Illustration of a centroid with a radius of three pixels converging over iterations 1 (left), 6 (center), and 12 (right). The blue point indicates the start point of the centroid, and the green point indicates the current centroid estimate.
In the text

	Fig. 2. Illustration of a centroid moving from time step 15 (left) to time step 20 (center) to time step 25 (right). The blue point indicates the start point of the centroid, and the green point indicates the current centroid estimate.
In the text

	Fig. 3. Illustration of the flows identified by applying the proposed methodology to the SDO image data from the 171 Å spectral line for the date 6 March 2012 from 12:22:00 to 12:35:00 referred to as the limb data (top) and again for the date 9 March 2012 from 03:00:00 to 03:24:00 referred to as the central data (bottom).
In the text

Fig. 4. Velocity lower bound ${\widehat{\mathit{v}}}^{\min }$ and upper bound ${\widehat{\mathit{v}}}^{\max }$ together with the classification level of significance 100(1 − α)%. The red, blue, and green (solid line) curves indicate a supersonic, sonic, and subsonic classification (left ordinate axes), respectively at a 50% level of significance. The dotted line curves show the level of significance (rights ordinate axes). The horizontal dashed line represents the background sound speed. These flows all have a majority supersonic classification.

In the text

	Fig. 5. Same as in Fig. 4 for subsonic flow patterns.
In the text

	Fig. 6. Same as in Fig. 4 for sonic flow patterns.
In the text

	Fig. 7. Same as in Fig. 4 for transonic flow patterns.
In the text