The flare event under study is a long-duration M9.2 solar flare located at the east limb, that occurred on March 7, 2015. We mostly focus on observations from the IRIS satellite, which, since its launch in 2013, has been providing unprecedented high-resolution images and spectra of the lower atmosphere and corona (De Pontieu et al., 2014; 2021). IRIS consists of: (1) a Slit-Jaw Imager (SJI) channel, acquiring images in four different filters showing plasma at photospheric, chromospheric and transition region temperatures with a 0.167″ resolution; and (2) a spectrograph channel, observing emission lines and continua formed over a broad range of temperatures (from the photosphere to the flaring corona), at very high spatial (0.33–0.4″), temporal (down to ≈ 1s) and spectral (2.7 km s⁻¹ pixels) resolution. The level 2 IRIS data used here are already calibrated and prepped, as described in the IRIS’ documentation and instrument papers (De Pontieu et al., 2014; Wülser et al., 2018).

The Hinode X-Ray Telescope (Golub et al., 2007; Kano et al., 2008) observed the flare with its flare response, using the Be-Thin and Be-Med filters and a resolution of 1.028″ per pixel. A two-filter observing program consisting of the Al-poly and Be-thin filters was taken before and after the flare response. The XRT data is calibrated and prepped using the xrt_prep routine, available in the SolarSoft suite of IDL programs (Freeland and Handy, 1998), which does the dark subtraction, and removes the pedestal and vignetting effects (Kobelski et al., 2014). The current version of xrt_prep also aligns the XRT images with AIA images by calibrating the time-dependent offsets between XRT and the Hinode Ultra Fine Sun Sensors (Yoshimura and McKenzie, 2015).

AIA (Lemen et al., 2012) provides full-Sun images with a resolution of ∼ 0.6″ per pixel and a cadence of 12 seconds for the EUV passbands. The AIA data are processed using the SolarSoft routine aia_prep, which de-rotates the images from the different AIA telescopes, aligns them, and gives them all the same plate scale.

Figure 1 shows an overview of the flare as observed by XRT in the Be-Thin filter (left columns), AIA in the 131Å filter (dominated by Fe XXI emission at ≈ 10MK, middle column) and in the IRIS SJI centered at 1330Å (right column). This SJI channel is dominated by emission from C II plasma at around 10–40kK, but it is also sensitive to hot plasma from Fe XXI during flares. The vertical lines in the XRT images indicate the location of the IRIS spectrograph slit during the observation. The white box in the middle images indicates the field of view of the IRIS SJI images on the right.

Figure 2 shows two images from the XRT flare response using the Al-Thick and Be-Thin filters. The right two panels of Figure 2 shows the temperature and emission measure, calculated from the filter ratio of the XRT two filters using the IDL routine xrt_teem_ch. We use coronal abundances to calculate the XRT temperature response functions because we are interested in the above-the-looptop plasma, which is likely directly heated, instead of resulting from chromospheric evaporation. Spectroscopic measurements in this region during another flare indicate that abundances are likely coronal (Warren et al., 2018).

In this work, we focus on the analysis of the Fe XXI observed by IRIS in high-temperature ≈10 MK plasma during flares. The IRIS observations during the flare event under study consist of an 8-step sparse raster (where the distance between consecutive slit positions is 1″), with a ≈30s exposure time, a raster cadence of ≈250s and a factor of two summing along the slit position. Figure 3 shows the intensity map of Fe XXI and spectra profiles at different times. For each slit position and time we averaged the Fe XXI spectra over 50 slit pixels to increase the signal-to-noise ratio. The centroid position (or Doppler velocity) of spectral lines and width (or FWHM) are obtained by using Gaussian fitting at each time. In fact, the line is very weak in the cusp/reconnection region, which is not surprising given the low density of the plasma there and the fact that the IRIS Fe XXI 1354.08Å line is a forbidden transition (e.g., Young et al., 2015). The non-thermal velocity was obtained by using the following formula: v nth = λ 0 c ⋅ 2 l n 2 ⋅ F W H M 2 − F W H M t h 2 − F W H M instr 2 ⋅ (1) where λ ₀ is 1354.08 Å, the rest wavelength of the line (e.g. Polito et al., 2015), c is 3 × 10¹⁰ cm s⁻¹, FWHM _instr is the IRIS instrumental width (26 mÅ) and FWHM _th is the thermal broadening of the line (≈0.43Å) assuming a formation temperature of logT [K] = 7.05.

Figure 4A shows that the non-thermal velocity gradually decreases over time from more than 100 km s⁻¹ to about 30 km s⁻¹ in around 2 hours. Figure 4B shows that the Doppler shift velocity appears to be close to a mean value(∼4 km⁻¹) for most of the time. The only exception is the first few minutes of the observation, where the Doppler velocity reaches larger values up to almost 40 km s⁻¹. We note that there is a lot of scattering between the values obtained in different IRIS slits. The IRIS level 2 data used in this work are already corrected for both the orbital and absolute calibration of the wavelength array, with an accuracy that is estimated to be ≈1 km s⁻¹ for most databases (Wülser et al., 2018). However, it is usually recommended to double-check the calibration manually using the centroids of strong neutral lines in case a small residual drift is still present. Such a sanity check cannot be performed in our data set because the raster is located off-limb and neutral lines are not observed.

Further, we note that the peak time of the M class flare occurred at about 22:22UT, but IRIS did not detect sufficient signals in the Fe XXI line before ∼22:55UT shown by the dashed line in Figure 4C. As shown by the light curves of GOES X-ray, RHESSI, SDO/AIA 131Å and Fe XXI flux in this figure, the eruption starts from about 22:00UT and continues for more than 4 h. On the other hand, high temperature plasma is already observed before this time in the AIA 131 Å channel, showing plasma at ∼10 MK. This lack of emission in the IRIS Fe XXI line may be caused by the fact that this line is particularly faint, as mentioned above. This also suggests that the non-thermal velocity might have been higher during the early impulsive phase of the flare, before we can observe any Fe XXI emission in IRIS.

We use the same 3D solar flare model as shown in Shen et al. (2022), which follows the classic CHSKP configuration of two-ribbon flares. In this model, we solve the initial and boundary value problem governed by resistive MHD equations using the public code Athena (Stone et al., 2008). We use the non-dimensional form of MHD equations and set the primary characteristic parameters as follows: length L ₀ = 1.5 × 10⁸m, magnetic field strength B ₀ = 0.001T, density ρ ₀ = 2.5 × 10¹⁴m⁻³, temperature T ₀ = 1.13 × 10⁸K, velocity V ₀ = 1.366 × 10⁶ m s⁻¹, and time t ₀ = 109.8s according to the general coronal environments. At the beginning of the MHD simulations, the magnetic field strength is uniform ( ∼ 10 Gauss) in the background region, and the plasma number density is around 8 × 10¹⁴m⁻³ near the bottom of the corona, above a denser layer representing the chromosphere. In the later time as shown in the following sections (e.g., Figure 5), the magnetic field strength slightly enhanced by ∼ 1.3 times in both flare foot-points and loop-top regions following the magnetic loops shrinking, while the plasma density can be ∼ 10 and ∼ 6 times higher in foot-points and loop-top regions, compared with the background corona. As the magnetic reconnection continually takes place, the plasma in the reconnection region and flare loop-top regions can be heated to higher than 10⁷ K at the most bursty phases and maintain the relatively high temperature ( > ∼ 8 MK) during the whole simulation period. It is worth mentioning that the plasma β, the ratio of the gas pressure to the magnetic pressure, is another dominant parameter in the MHD simulations and is set to about 0.1 in the background corona. The characteristics (such as magnetic strength and plasma density) are then scaleable variables by maintaining a constant β.

We perform combined 2D and 3D models in which the 2D simulation runs first to build the classic flare geometry and then the 3D model is started based on the resulting 2D configuration. In 2D cases, the system is initialized from a pre-existing vertical Harris-type CS along the y direction in mechanical and thermal equilibrium. Driven by the initial perturbation on magnetic fields (also see Shen et al., 2022), the CS becomes thinner due to the Lorentz-force attraction, and a pair of reconnection jets flow away from the reconnection X-point close to the initial perturbation center. The closed magnetic loops then appear as the reconnected magnetic flux accumulated at the bottom due to the line-tied boundary, in which the magnetic field lines are rooted at the boundary. Once the flare loops are well-formed after about t = 17.5t ₀, we start the 3D simulations by symmetrically extending all primary variables from the 2D domain to the third direction (z −). We then apply this time as the starting-point of evolution in the following analysis. The system then self-consistently evolves, and a set of reconnection-driven flare phenomena are seen, including Alfvénic bi-directional reconnection outflows, the termination shocks, and the turbulent interface region below the extended reconnection current sheet and above the flare loops (Figure 5A).

To compare our models with observations, we calculate synthetic Fe XXI 1354 Å line profiles that are observable by IRIS in high-temperature flare plasma (e.g., ∼11 MK). Once we compute the plasma properties (temperature, density, and velocity) on each cell from the 3D MHD simulation, the synthetic emission can be obtained along any chosen LOS (e.g., Lines A and B in Figure 5) by using the formula (e.g., Guo et al., 2017): I ν = h ν 4 π ∫ f ν n e n H g T e d l (2) where ν indicates the frequency of emission lines, and l is the integration path. n _e , n _H ∼ 0.83n _e are the electron and proton densities, which are obtained from MHD plasma density in the following calculations. g (T _e ) is the contribution function, which can be obtained from the CHIANTI database (Del Zanna et al., 2021) but needs account for proper element abundance (e.g., Feldman, 1992, in this paper). f _ν indicates the variation of velocity distribution function due to plasma flows along the LOS, given by f ν = 1 π 1 / 2 Δ ν exp − Δ ν + ν 0 v l / c Δ ν D 2 . (3) Here v _l is the plasma flow speed along the LOS, c is the speed of light. ν ₀ is the rest frequency, and Δ = ν − ν ₀ means the offset frequency accordingly. The thermal broadening can be calculated by Δ ν D = ν 0 c 2 k T m , (4) where k is the Boltzmann constant, T is the temperature, m is the atomic mass of the chosen ion (e.g., Fe in this work).

We first analyze the spectral properties along the sampling Lines A and B. As shown in Figure 5, Line A is located at a relatively higher altitude along the CS, which is close to the primary reconnection X-points indicated by the strongest current density (J _z , in Figure 5B) and flow stagnation regions (Figure 5C). On the other hand, Line B is chosen to go through highly turbulent above loop-top regions, where the SADs-like tenuous down flow features are well-developed as reported by Shen et al. (2022). At this time (5.2t ₀), Line A shows slightly higher temperature and lower density than that of Line B overall (Figure 6). Along the LOS, we note that the plasma is appreciably cooler than the expected Fe XXI 1354 Å line temperature (e.g., Log T ∼ 7.05 K). This effect occurs because the LOS path is outside the high-temperature CS at some positions. Due to the highly turbulent plasma flows, the temperature, density, and velocity along the sampling lines all show various perturbations. The domain perturbation velocity ranges from ∼ + 100 to ∼ − 100 km s⁻¹ as shown in Figure 6C.

Figure 6D shows the synthetic line profiles of the Fe XXI 1354 Å lines for the above two sampling lines. We assume that LOS is from z = −0.25L ₀ to z = +0.25L ₀ along each sampling line, and the Doppler shift velocity is then plotted based on this geometry accordingly. We note that the LOS direction may also be reversed in the observations, so the red and blue shifted Doppler velocity is a relative value in this plot. In this case, Line A shows a larger Doppler shift velocity (about 18 km s⁻¹) than that of Line B because the Alfvénic downward reconnection outflows around Line A dominate the overall flow behavior. We measure the line broadening width using the full width at half maximum (FWHM) of the synthetic emission profiles. Then the non-thermal velocity (v _nth ) can be obtained according to the following equation: FWHM = 4 ln 2 λ 0 c 2 k T e q m + v nth 2 , (5) where the λ ₀ is the central wavelength of Fe XXI 1354 Å line, logT _eq = 7.05 K is the typical formation temperature of this line. In this particular position, the non-thermal velocity of Lines A (∼95 km s⁻¹) and B (∼37 km s⁻¹) are consistent with the deduced values in recent observations (e.g., Warren et al., 2018), but the line A value is slightly larger than those found by Doschek et al. (2014). However, the synthetic line profile strongly depends on the space and time evolution of the solar flare system, which will be addressed in the next section. In the current analysis in this paper, we only focus on non-thermal broadening features but leave the Doppler shifts for future works because the integration depth along the LOS is limited to the simulation box scale, which is surely shorter than the realistic CME/flare scale. Thus, the synthetic Doppler shift of Fe XXI 1354 Å line can be significantly affected by several locally large momentum components (e.g., the well-developed flux ropes), and the model-deduced Doppler velocity can largely depart from the common observational results by IRIS (Figure 4B). Therefore, large-scale solar eruption models are required to make a detailed comparison with the observable Doppler velocity in the future.

In general, the actual temperature on each pixel (or cell) along the LOS cannot always be exactly at the T _eq , in either the observations and MHD models. Therefore, the deduced non-thermal velocity is unavoidably affected by the thermal broadening due to different temperatures, especially for high-temperature plasmas with logT > 7.05 K. In the following sections, we will follow this approach (where the equilibrium temperature is usually assumed when calculating the non-thermal broadening) to match the observational data analysis but give a more detailed discussion in Section 3.3.

We analyze the perturbation properties of primary variables along the LOS by using the Fourier transform method to investigate the nature of the broadening of spectral lines. In Figure 7, we perform the one-dimensional Fourier transform for density (ρ), velocity square (v ²), and kinetic energy (E _k ) along the Lines A and B. It is clear that these perturbations appear on all scales (or wave number (k) ranges), and the kinetic energy cascades from the large scale to the small scale as well. Furthermore, the spectrum matches the power-law tendency as the prediction in the classic turbulence theories, which could naturally result in the line broadening. Due to the limitation of the MHD grid sizes, we do not show here the spectrum distribution in high-k ranges where the possible inertial range is expected in turbulence theories. However, we can estimate the turbulence properties in the intermediate ranges quantitatively by fitting a power-law spectrum (∼ k ^−γ ) in Figure 7C. We also check the temperature weighted E _k in according to the chosen temperature (logT ∼ 7.05 K) for the Fe XXI 1354 Å line. A similar power law distribution pattern can be seen clearly, as shown in Figure 7D. We notice that Line A with larger Fe XXI 1354 Å line broadening shows a higher index (γ ∼ 5.34) as compared with the shallower line B (γ ∼ 2). This result indicates that large bulk flows can significantly contribute to the line broadening as shown by the velocity spectrum of Line A, in which the more turbulent flows remain in the low-k (or large scale) ranges.

In Figure 8, we plot the distribution and temporal evolution of synthetic Fe XXI 1354 Å emissions with a viewpoint such that the plasma sheet is viewed edge-on. At each MHD simulation cell on the x − y plane, we obtain the synthetic spectral lines assuming that the LOS is along the z − direction. The maximum intensity of the obtained Fe XXI 1354 Å spectra line at each MHD grid is shown in Figure 8A, in which the bright sheet structure and flare cusp regions can be clearly seen as predicted by the standard solar flare model. The non-thermal velocity is calculated and shown in Figure 8B.

Figure 8 shows that the region where there is strong non-thermal line broadening (up to ∼ 300 km s⁻¹) is dynamically evolving in both the CS and flare cusp regions. The region of high non-thermal broadening flows downwards along the CS from the primary reconnection X-point region to the flare cusp region. The white arrows in Figure 8B illustrate this process during about one characteristic timescale from t = 4.6t ₀ to 5.5t ₀. The downward-moving patch of high non-thermal broadening finally consolidates into the cusp region after t ∼ 5.5t ₀, and causes extended regions of high non-thermal broadening above the flare loop-top.

In order to understand the magnetic field topology and plasma flow properties in these high-v _nth regions, we plot out the horizontal momentum component (ρV _z ) and magnetic component (B _z ) in Figure 9. Because the emission intensity of the Fe XXI line is proportional to ρ ² and is primarily affected by horizontal flow V _z , the plasma momentum (ρV _z ) should be a sensitive factor to distinguish the effects of different structures on emission line profiles and line broadening. As shown in Figure 9A, the dashed lines and arrows highlight the heights characterized by enhanced v _nth , similarly to Figure 8B. It is clear that the position of the strong horizontal momentum component coincides with such heights. The strongest momentum regions are usually associated with well-developed magnetic flux-ropes. At time = 4.6t ₀, two of these flux-ropes are highlighted by plotting the chosen helical magnetic field lines where the strongest momentum component appears. At later times (4.9 ∼ 5.2t ₀), these two flux-ropes are associated with strong ρV _z components that move downwards to lower heights. In fact, the strong momentum components (ρV _z ) are due to the appearance of a guide field B _z in such a turbulent reconnection current sheet. As shown in the second row of Figure 9, the guide field B _z widely appears with a very turbulent behavior inside the whole current sheet. The typical relative strength of B _z to the total background (B) ranges from ∼5%–10%, and can be larger than ∼30% at particular positions and times (indicated by the dark-blue color in Figure 9). It is then not surprising that the well-developed flux-ropes with strong momentum components (ρV _z ) are generally found in high B _z regions. In addition, other structures all have significant contributions to the momentum component ρV _z , including the growing flux-ropes and very turbulent reconnection outflows with remarkable horizontal flow components along the z − direction. As shown in Figure 9, well-developed flux-ropes are relatively rare in the whole current sheet region. The fluctuations in ρV _z can be commonly found around the growing flux-ropes and turbulent reconnection outflow patches as well. For example, a large number of growing flux-ropes appear at the center region at time = 5.2t ₀, which is consistent with the high v _nth region above y = 0.7L ₀. Meanwhile, the nearby fully developed flux rope is located at a much lower altitude (y ∼ 0.65L ₀).

The downward high-v _nth features commonly exist for a longer duration in the reconnection process. Figure 10 displays the time–distance map (or “stack plot”), showing the distribution of the non-thermal velocity of Fe XXI 1354 Å line along at the CS center (x = 0) as a function of time. The high-v _nth emission flows away from the primary X-point regions nearby y ∼ 0.9, as can be clearly seen in Figure 10B. These high-v _nth patches are most like to originate from the reconnection X-point sites (around y ∼ 0.9L ₀), and spread with the reconnection outflows to the lower end of the CS. At different altitudes along the CS, there are two obvious regions of high non-thermal broadening: one is close to the primary reconnection X-point site and another is located above the flare loop-top. Above the flare loop-top region, the non-thermal velocity is strong but highly depends on the more complex plasma dynamics in this interface region. It could be roughly the same or decrease with height, as shown in Figure 10. We also notice that high non-thermal-broadenings also appear in upward flowing plasma near the upper boundary (y ∼ 1.0).

In general, both plasma turbulence and bulk flows can contribute to the non-thermal widening (e.g., Ciaravella and Raymond, 2008), and it is hard to quantitatively explore the turbulent strength based on the observational line profiles due to complex plasma environments along the LOS direction. On the other hand, the numerical models allow us to obtain detailed information about the local plasma flows and integrated effects on spectral line profiles. Therefore, it is interesting to examine the correlation between the non-thermal velocity and plasma turbulence over a sample of locations in both the CS and flare loop-top regions during a relatively long period from 4.6 ∼ 7.4 t ₀ in this model. We introduce a non-dimensional parameter V turb ≡ 1 n z ∑ v z i 2 / V A to describe the turbulence strength of plasma flows along the LOS. Here, v z i ≡ v z − v z ̄ is the turbulent fluctuation of each cell to the mean velocity v z ̄ , V _A is the characteristic Alfven speed, and n _z is the total cell number along the LOS. Figure 11 plots the two-dimensional population histogram of V _turb and non-thermal velocity (v _nth ). In Figure 11A, the background colors indicate the population of samples on the v _nth - V _turb map. An overarching feature of these plots is that the v _nth distribution is clearly proportional to the turbulence strength (V _turb ) because stronger plasma turbulent flows can naturally cause larger line-broadening, as discussed in the previous sections. However, the detailed distribution is more widely distributed with the increase of V _turb , and will not likely be fitted by using one simple linear line. Here, we annotate one typical growth direction using a red dashed line and a deviation direction using a black dashed line in Figure 11A, where the most dominant samples appear. The samples around the red line show higher non-thermal velocities with comparable turbulence strength compared with those around the black line.

The deviation among the black and red dashed lines in Figure 11A could be due to the variation of plasma flows and turbulence properties at different heights, especially in the CS region and flare cusp regions where the plasma density and plasma β environments are largely different. Therefore, we display all samples on the v _nth -V _turb map with their height information, in which colors indicate y − position of each sample in Figures 11B–D. We separate all samples into three groups (y ≤ 0.5L ₀, 0.5L ₀ < y < 0.65L ₀, and y ≥ 0.65L ₀) to make the difference more visible.

Figure 11B shows the first group of samples with the lowest heights below y ∼ 0.5L ₀, in which both the v _nth and V _turb are small. Because this group of samples is close to the more dense flare loop-top regions, where the turbulent flows are expected to be relatively lower than those in the above-the-looptop regions, v _nth is naturally small (e.g., < ∼ 100 km/s). The samples in the second group (Figure 11C) in the very turbulent flare cusp regions have stronger turbulence and larger non-thermal velocities.

Similar high-v _nth and high-V _turb can be seen in the CS region as shown in Figure 11D as well. Samples in both the CS and the above loop-top regions tend to follow the growth direction marked by the red dashed line in Figure 11A, except small abnormal patches with low v _nth and high V _turb that usually appear in the low end of the CS (y < ∼ 0.8L ₀) and upper the loop top region (y > ∼ 0.57L ₀). In fact, these abnormal samples that have departed from the red line direction may indicate strong bulk flows which will be discussed in the following section.

It is worth noticing that thermal broadening also contributes to the emission line broadening and has an impact on the corresponding v _nth distribution, especially when the plasma is hotter than the equilibrium temperature of Fe XXI 1354 Å at logT _eq ∼ 7.05 K. In the above investigations, there is a set of samples in our model that are slightly lower in temperature than the Fe XXI 1354 Å temperature (T _eq ) (e.g., see Figure 6A). Therefore, we investigate the v _nth distribution in different temperature environments. We maintain the plasma density the same as in the above analysis but scale the temperature to match the observations (e.g., logT 7.0 ∼ 7.2 K shown in Figure 2). In the non-dimensional MHD models, the characteristic temperature can be expressed by T 0 = B 0 2 β 0 / ( 2 μ ρ 0 R ̃ ) and Alfvén velocity is V A = B 0 2 / ( μ ρ 0 ) . Here, β ₀, R ̃ , and μ are background plasma β, gas constant, and magnetic permeability, respectively. Thus, the scaling of temperature (T ₀) with corresponding changes in characteristic velocity (V _A ) allows the MHD model to be scalable by maintaining the plasma β ₀ as the same. Therefore, we increased the simulated temperature to a factor of 1.4× higher than that in the original MHD model result, and the velocity is also increased by a factor of 1.4 accordingly. As a result, the samples in the model can cover a higher temperature range from logT ∼ 6.9 to ∼7.2K.

Figure 12A shows the population histogram of V _turb and v _nth in the scaled temperature case. The red and black dashed lines are exactly the same as in Figure 11A as well. The distribution feature of v _nth is consistent with the above results in the original MHD temperature case (T _MHD ), though the absolute values of v _nth are slightly larger by about 10% ∼ 20%. These larger v _nth are partially caused by the scaling process of MHD models, because the characteristic velocity increased by about ∼ 18 %, which could lead to larger Doppler shift velocities and corresponding larger line broadenings. The mean temperature of each sample is displayed by different colors in Figure 12B. It is interesting to note that the highest temperature (logT ∼ 7.2K) does not necessarily indicate the largest v _nth if the turbulence strength remains low (the dark-red samples), which suggests that the plasma perturbation behaviors have more crucial impacts on the non-thermal velocity distribution than the thermal broadening itself.

Figures 12C–E shows Fe XXI 1354 Å line profiles in two typical regions: CS region (S1, S3), and the above flare loop-top region (S2). Sample S1 represents the most dominant features with low turbulence strength and low non-thermal velocity. The line profile basically matches the Gaussian-type shape with minor line broadening due to the relatively weak turbulent flows. Sample S2 is at the high v _nth end with strong V _turb , where the strong plasma perturbations contribute to a remarkable line-broadening in both the 1.4× higher and slightly lower T _MHD cases. In addition, the plasma bulk flows also cause the second emission peak at ∼ 1354.4 Å which makes the line profile depart from a Gaussian distribution and causes a wider v _nth . Sample S3 indicates the abnormally high V _turb with low non-thermal velocity. The reason can be found in the line profiles which show two separate emission peaks: one is for the dominant Fe XXI 1354 Å with a smaller Doppler velocity and the other one is centered at ∼ 1355.4 Å due to the enormous plasma bulk flow, such as the well-developed flux-ropes as shown in Figure 9. In general, the turbulence strength (V _turb ) only refers to mean perturbation features, and it could be too simple to represent complex flow properties including randomly turbulent flows and large bulk flows. Therefore, the appearance of these abnormal samples around the black dashed line in Figure 12 suggests that the investigation of plasma turbulence based on the emission line broadening features must consider the different fine structures in the CS (e.g., plasma blobs or flux-ropes) and flare loop-top regions (e.g., macro plasma instabilities or SADs in Shen et al. (2022)).

In this section, we investigate the averaged Fe XXI lines profiles across the CS and compare them to the IRIS observations during the 2015–03–07 flare. We set the position of the simulated IRIS slits above the flare loop top region in the MHD model as shown by Figure 13A. The synthetic SDO/AIA 131 image in this plot is obtained by assuming that the LOS is along the z − direction. The synthetic bright flare loops, cusp region, and extending CS regions all match well the SDO/AIA images of the solar flare under study (Figure 13A) well. It is clear that the IRIS slits are located just above the bright flare arches in AIA 131 during this eruption event (also see Figure 1). Therefore, we chose a height of y ∼ 0.62L₀ in the MHD model to simulate the IRIS raster in the following analysis. As the reconnected magnetic flux is accumulated at the solar surface, the flare loop system gradually grows in both vertical and horizontal directions. Thus, the relative position of the IRIS slits to the flare cusp region gradually changes as a function of time. Because our MHD model is focused on the fast magnetic reconnection process in solar flares, it is reasonable to compare the predicted SDO/AIA emissions with the observed flare just after the impulsive phase when the post-flare loops are well-formed but the bursty magnetic energy release is still on-going. To simplify the analysis, we show a particular SDO/AIA observation at the starting time of IRIS observations (2015–03–07T23:00) in Figure 13B, and keep the simulated IRIS slits in the same position in this work.

Along each slit in Figure 13A, we average the Fe XXI 1354 Å line profiles over all MHD cells with I > I peak e . Here I is the maximum intensity of the Fe XXI line at each cell, and I _peak indicates the peak I across the CS on this slit. The non-thermal velocity of each slit is then calculated using the averaged line profile, assuming logT = 7.05K for the strongest Fe XXI emission. Figure 13C shows the non-thermal velocity variation with time on the chosen slits, No. 2, 4, and 6. An overall feature is that the v _nth variation tendency on these three slits is consistent with each other. For example, the local maximum v _nth peak (red shadowed region) appears initially at times 5.2 on slit two and then at 5.4t ₀ for slit 6. A similar tendency is also clear around time 6.8t ₀, as highlighted by the blue shadowed region. By comparing the v _nth trends above with Figure 10, one can see that the first v _nth peak (around times 5.3t ₀) is consistent with downwards moving structures with enhanced Fe XXI line widths between t = 4.6 to t = 5.5t ₀, as indicated by the yellow dashed line. In the later times (t = 5.5 to t = 7.4t ₀), the relatively weaker downward structures also caused the second local v _nth peak at 5.9t ₀ and possibly the third peak at 7.4t ₀ on Figure 13C.

In the ideal situation when the non-thermal velocity variation is due to the passing downward structures, the characteristic size of detectable downward patches can be estimated as follows: δ l + d slit cos θ ∼ v d × δ t . (6) Here δl is the length of moving patches, d _slit is the interval of the IRIS slits, θ is the intersection angle between slits and the solar surface, v _d is the moving speed of patches, and δt is the exposure time of each slit (or interval time between two slits). For instance, δt ∼ 30 s for the IRIS Fe XXI observations, d _slit is about 726 km, and θ ∼ 30°. Assuming the downward moving structures are in sub-Alfvénic speed (e.g., 300 km/s), a moving enhanced v _nth patch with the size of δl ≥∼ 8000 km should be well recognized. The above estimated δl equals ∼ 0.05 L 0 , which is reasonable and matches the modeling structures as shown in Figure 10.

Similar to the synthetic results in MHD models, we also see very similar v _nth evolution trends in the IRIS observations (Figure 13D). The red and blue shadowed regions show one example where the local maximum and minimum v _nth has been recognized first on slit 2, then on the following slits 4, and 6. Remarkably, the evolution of v _nth of the Fe XXI line is comparable between MHD model predictions and IRIS observations, which range from ∼ 20 km/s to 180 km/s during this particular period. We note that the time-resolution of the deduced v _nth variation from the IRIS is limited by the relatively long exposure time ( ∼ 30 s). Therefore, Figure 13D may unavoidably smooth out fine structures and cannot display the short-period perturbations that are visible in the synthetic modeling results.

Due to the similar variations in the MHD predictions and the IRIS observations described above, it is likely that the v _nth variation on IRIS Fe XXI lines is due to intermittent downwards moving structures along with the reconnection outflows as shown in the above models. In fact, we find very similar intermittent downwards-moving signals on SDO/AIA 131Å images. Figure 14 shows the time-distance map of AIA 131Å intensity above the flare loop top regions. In this plot, we count all emissions over 16 pixels around the CS region (the bright plasma sheet on AIA 131Å maps) at each height to increase the signal-to-noise ratio. Two typical downward features are fitted by using red and yellow dashed lines. These downflow features range from ∼ 250 km/s to ∼ 100 km/s, which is consistent with the model-predicted moving speed of turbulent structures as shown in Figure 10. We also overlay the non-thermal velocity trends on the AIA 131Å stack map as shown by the solid-dots lines. There is no clear correlation between the brightest AIA features and the highest v _nth because the AIA intensity is mainly dominated by plasma density in the CS region, while the v _nth values are more affected by the turbulent plasma flows.

We notice that the variation feature of non-thermal velocity above the flare loops (e.g., in Figures 11, 14) are consistent with the previous observational studies. The increasing behaviors on non-thermal velocity from the flare loops to higher altitudes have been reported in the past. For example, Doschek et al. (2014) found that the non-thermal velocity in the multimillion-degree regions increases with height above the bright coronal flare loops from Hinode/EIS observations. Our simulation results also match their explanation, in which the strong newly contracting hot and turbulent plasma mainly contributes to these high non-thermal motions.