Current-sheet Oscillations Caused by Kelvin-Helmholtz Instability at the Loop Top of Solar Flares

Current sheets (CSs), long stretching structures of magnetic reconnection above solar flare loops, are usually observed to oscillate, their origins, however, are still puzzled at present. Based on a high-resolution 2.5-dimensional MHD simulation of magnetic reconnection, we explore the formation mechanism of the CS oscillations. We find that large-amplitude transverse waves are excited by the Kelvin-Helmholtz instability (KHI) at the highly turbulent cusp-shaped region. The perturbations propagate upward along the CS with a phase speed close to local Alfv\'{e}n speed thus resulting in the CS oscillations we observe. Though the perturbations damp after propagating for a long distance, the CS oscillations are still detectable. In terms of detected CS oscillations, with a combination of differential emission measure technique, we propose a new method for measuring the magnetic field strength of the CSs and its distribution in height.


INTRODUCTION
Solar flares are one of the most energetic phenomena in the solar atmosphere and often appear as a sudden emission enhancement over the whole electromagnetic spectrum. In the past decades, their energy release mechanism, spatial structures, and dynamical properties have been widely studied. The well-known standard flare model summarizes several key features of flares, including two parallel bright ribbons, cusp-shaped loop-top structure, elongated current sheet (CS), and the erupting magnetized plasmoid (CSHKP ;Carmichael 1964;Sturrock 1966;Hirayama 1974;Kopp & Pneuman 1976;Shibata et al. 1995;Lin & Forbes 2000). These features are all closely related to the fundamental energy release process, i.e., magnetic reconnection.
Above the flare loops, the thin and long stretching or ray-like high-temperature structures, typically of 10-20 million Kelvins, are suggestive of model-predicted CSs (see Ciaravella et al. 2002;Webb et al. 2003;Ko et al. 2003;Liu et al. 2010;Patsourakos & Vourlidas 2011). For a few events, it is observed that the thin and long CSs even extend to the high corona for serval solar radii and connect to the erupting coronal mass ejections (e.g., Cheng et al. 2018). Because of the large Lundquist number of the corona, plenty of small-scale plasmoid structures are generated by tearing mode instability and then propagate upward and downward along the CSs (e.g., Shen et al. 2011;Mei et al. 2012;Lin et al. 2015;Lee et al. 2021).
Interestingly, the CSs usually also present a transverse oscillation. Chen et al. (2010) reported propagating largescale waves along a helmet streamer CS with a period of one hour and wave-length of several solar radii. Similar waves with a shorter period of 30 min were also captured by Ling et al. (2014). Samanta et al. (2019) analyzed a similar but more slowly propagating wave (∼ 19 km/s) along the ray-like CSs and proposed that they might be vortex shedding in the corona. Verwichte et al. (2005) found that the supra-arcade downflows (SADs), probably corresponding to the reconnection outflows, as well as the bright CS-like supra-arcade fans (SAFs), also oscillated but the periods were usually much shorter, mostly < 5 min. The first evidence for the CS oscillations was reported by Li et al. (2016), who found that it had a period of ∼ 11 min and propagated away from the sun with a speed of ∼ 200 km/s. Figure  1 exhibits a similar CS oscillation event captured by the Atmospheric Imaging Assembly (AIA) on board the Solar Dynamic Observatory (SDO) (Pesnell et al. 2012;Lemen et al. 2012) on 2011 Oct. 22. One can clearly see that the ray-like structures above the flare loops, supposed to be the CSs, sway in the plane of the sky. The CS oscillations might be related to the fine processes at the cusp-shaped flare loop region where turbulent flows and plenty of shocks take place (see McKenzie 2013;Takasao et al. 2015;Shen et al. 2018;Cai et al. 2019;Ye et al. 2020). Takahashi et al. (2017) reported that, under high-Lundquist-number conditions, the quasi-periodic oscillations were driven by the horizontal motions of termination shocks (TSs) with oblique fronts at the flare loop top and the buffer region at the CME bottom. Recently, Xie et al. (2021) performed a similar simulation and found that the oscillations are caused by Rayleigh-Taylor instability (RTI) in the buffer region with the break of symmetry and that the oscillations propagated towards the sun and caused prominent displacements of the CS structure.
In this letter, we further explore the origins of the CS oscillations based on a high-resolution 2.5-dimensional magnetic reconnection model. We found that the Kelvin-Helmholtz instability (KHI) of the TS tail flow initializes the development of asymmetry flows and thus gives rise to the CS oscillations. Furthermore, we analyze the propagation properties of perturbations causing the CS oscillations and confirm the method we propose for estimating the magnetic field strength of the CS and its distribution in height. Section 2 briefly introduces our numerical model. The main results are presented in Section 3, which are followed by a summary and discussion.

NUMERICAL MODEL
Our numerical model is based on the resistive MHD equations including gravity, anisotropic thermal condition, radiation cooling, and background heating: where, P * = p + B 2 /2, e = p/ (γ − 1) + ρu 2 /2 + B 2 /2, γ = 5/3, n e and n i are respectively the number density of electrons and ions, and other variables are denoted by standard notations. The conductivity parallel with magnetic field is determined by κ = κ 0 T 2.5 , where κ 0 = 6.67 × 10 5 erg · s −1 cm −1 K −3.5 . The gravity acceleration is calculated by g = −gê y , where g = g 0 / (1 + y/R ) 2 , R denotes the solar radii, and g 0 = 2.7390 × 10 4 cm/s 2 . We use a widelyused optically thin radiation cooling function Λ (T ) (see Klimchuk et al. 2008;Ye et al. 2020;Shen et al. 2022). The background heating, H = n i n e Λ (T cor ), is supposed to maintain the initial energy balance in corona and also keep balancing the cooling effect in the bottom chromosphere region (y < 0.2), where T cor is the initial coronal temperature.
In this paper, all variables are normalized according to constant units identical to Wang et al. (2021). Physical units of main variables are listed in Table 1. The initial temperature of the gravitationally stratified atmosphere is set as

Units
where, T cor = 0.1, T chr = 0.002, h chr = 0.12, and w tr = 0.02. The initial pressure profile, calculated by p (y) = p 0 exp − y 0 g/T dy , is set to balance gravity (see also Ye et al. 2020), where p 0 = p ref exp g/T dy and p ref = 0.08 is the pressure at y r = 0.22. The initial force-free magnetic field, identical to Wang et al. (2021), forms a vertical CS resembling that in the CSHKP model. To trigger the reconnection, we use a localized anomalous resistivity near y = 0.5 and x = 0, which is also identical to Wang et al. (2021). The initial velocity is set to zero everywhere. The initial static equilibrium of the atmosphere is exactly satisfied in the corona region (y¿0.2), but is slightly perturbed in the transition region (0.1¡y¡0.2) where the thermal conduction term causes a localized downward energy flow due to the temperature gradient. However, this initial perturbation in the transition region is ignorable compared with the dominating reconnection process (see the animation of Figure 2). The boundary conditions are arranged as follows. The left (x = −0.5) and right (x = 0.5) are free boundaries, the top (y = 2) is no-inflow boundary, and the bottom (y = 0) is symmetric boundary.
The above system is simulated with the Athena++ code (Stone et al. 2020). We use the HLLD Riemann solver (Miyoshi & Kusano 2005), the 2-order piecewise linear method (PLM), and the 2-order van Leer predictor-corrector scheme to solve the conservation part of Eq. 1. The resistivity, thermal conduction, gravity, radiation cooling, and heating terms are calculated by the explicit operator splitting method. The 2-order RKL2 super-time-stepping algorithm is adopted to reduce computational costs (Meyer et al. 2014). Uniform Cartesian grids are used on both directions, namely, 1920 and 3840 grids on x and y directions, respectively. The pixel scale is ∆L = ∆x = ∆y = 26 km. The maximum simulation time in our simulation is t max = 12, which corresponds to 23 minutes in physical time.

Overview
The magnetic reconnection is initiated at a small region that has an anomalous resistivity. At t = 5, a flare loop system and an erupting plasmoid above the principal X-point appear (see also Yokoyama & Shibata 2001;Takasao et al. 2015). Under the small background resistivity, plasmoid instability dominates the reconnection and the CS is fragmented into various magnetic islands of different scales (see also Bhattacharjee et al. 2009;Shen et al. 2011;Mei et al. 2012;Ni et al. 2016;Ye et al. 2019;Kong et al. 2020;Zhao & Keppens 2020). The reconnection also shows a quasi-periodic characteristic. At t = 7.5, a relatively large magnetic island forms in the CS, it then moves downward and collides with the flare loop top (see the animation of Figure 2). As the island annihilates there gradually, the kinetic and magnetic energies are further released (see Wang et al. 2021). At the same time, the cusp-shaped loop top becomes highly turbulent. Following this big island, several relatively small islands are generated and also enter the loop top region. After t = 9, large-size islands are rarely generated and the oscillations of CS gradually appear. At t = 9.1, under the TS, the y-asymmetry of the flows starts to grow (see the 1st column of Figure 2). The tip of the cusp-shaped loop top show prominent oscillation after t = 9.7 (see Figure 2). The CS bottom sways back and forth, and the transverse waves are excited and propagate upward at the same time. To compare with the observations, following the method used by Xie et al. (2019) and Ye et al. (2020), we synthesize the AIA 94Å images by I 94 = n 2 e f (T ) dl = n 2 e f (T ) L los , where f (T ) is the temperature response function, and L los = 10 9 cm is the scale of the line-of-sight on z-direction. The wave-like swinging of the CS can be clearly observed in the synthesized AIA 94Å images (See Figures 2(d1)-(d6)).  Figure 3). At t = 9.1, under the TS front (the downstream region), the downflow speed is still strong (see Figure  3(d1) and (e1)). Once the tail flow is blocked by the relatively stable flare loops, two backflows form on both sides (e.g., Figure 3(c1) and (d1), also see Takasao & Shibata (2016)). Consequently, two shear layers form below the TS, where the KHI is triggered.

Initialization of CS Oscillations
Theoretically, if the magnetic field is strong enough, the magnetic tension along the shear layer can stabilize the perturbations and thus suppress the KHI (Jones et al. 1997). The key criterion of the KHI in magnetized plasmas is the local Alfvénic Mach number of the velocity transition in shear layers (e.g., Ryu et al. 2000), which is defined as √ ρ is the projected Alfvén speed, and B is the magnetic strength along the shear layer. To be specific, the KHI is stabilized by the magnetic tension for M Akh < 2, and if M Akh > 4, the magnetic field is too weak and the evolution of the KHI is almost fully hydrodynamical (see Ryu et al. 2000).
To estimate M Akh of two shear layers, we set a slit perpendicular to the shear layers (see Figures 3(b1)-(d1)) and extract the profiles of the parallel velocity (here u y ) and C A along the slit (Figure 3(e1)). U 0 is the difference between the speeds of sheared flows (e.g., the velocity difference between the orange and purple dots in Figure 3(e1)). Because C A varies along the slit, the average value across the shear layers is used when calculating M Akh . As shown by Figure  3(e1), the values of M Akh on the two shear layers (left and right) are 3.26 and 2.80, respectively, which means that the two layers are both unstable to the KHI and their evolutions are not symmetric about y-axis at t = 9.1. Subsequently, the plasma under the TS starts to stir significantly and the flows show different behaviors on both sides of y-axis. At t = 9.46, the previously y-symmetric TS front becomes oblique (see Figure 3(d2)). Both shear layers are still unstable to the KHI (see Figure 3(e2)). The left one starts to rotate after this moment, and, at t = 9.66, it evolves into a horizontal layer that still satisfies the condition of the KHI (see Figures 3(c3)-(e3)). Thereafter, similar behaviors of shear layers are observed and the oscillation amplitude of the CS bottom also grows.
We now examine the condition of the RTI in the loop-top region. The RTI can be switched on if the perturbation wave along the density interface satisfies (Hillier 2016;Carlyle & Hillier 2017;Xie et al. 2021) where k and ω denote respectively the wave-number and frequency of the perturbation, g ⊥ is the acceleration of gravity perpendicular to the interface, B is the magnetic field strength parallel with the interface, and ρ u and ρ l are the upper and lower density, respectively. Equation 3 means that the perturbation wave-length λ = 2π/k should be larger than the critical length λ c = 4πB 2 /g ⊥ (ρ u − ρ l ) for the sake of initiating the RTI. The density interface where the RTI might grow satisfiesê y · ∇ρ > 0 to guarantee ρ u − ρ l > 0 and its normal vector isn rt = ∇ρ/ |∇ρ|. We can numerically calculate the local value of λ c after defining B 2 ≡ B 2 x + B 2 y − (B ·n rt ) 2 , g ⊥ ≡ |g ·n rt |, and ρ u − ρ l ≡ |∇ρ| ∆L. Because the perturbation wave-length should be smaller than the scale of the loop-top region, we locate all the positions satisfying λ c < 0.1 and find that they show a highly scattered distribution (see the white dots in Figure 3(a1)-(a3)). Moreover, moving with turbulent flows, these isolated positions vary dynamically in time (see the animation of Figure 3), which means that their lifetime might be too short to initiate the RTI. In contrast, the main shear layers triggering the KHI typically extend ∼ 0.1L 0 continuously in space. As the development of KHI, the shapes of the shear layers vary gradually in time but their spatial scales are approximately maintained (see Figure 3). Meanwhile, in our previous simulations in which gravity is not included and the RTI effect is ignorable, the loop-top oscillations can still be observed (see Wang et al. 2021). Therefore, we can conclude that the KHI dominates the symmetry breaking of the loop-top region and thus the initialization of CS oscillations.

Propagation properties of CS Oscillations
We further quantitatively analyze the propagation properties of the CS oscillation. Based on the method introduced in Appendix A, we determine the central positions x cs (y) of the CS (see Figures 2(b1)-(b6) and 4(b)). The distancetime plot of the positions clearly shows some wave-like structures (Figure 4(a)). During the first three oscillations, both the amplitude and the period grow (see Figure 4(c)), and the maximum oscillation amplitude reaches ∼ 0.01. Though damping as propagating upward, the oscillation can still be clearly detected at a high altitude (see Figure  4(a)).
To estimate the propagation speed of the perturbations, we take six contours in Figure 4(a) and denote them by y i p (t), where i = 1, 2, · · · , 6 (see Figure 4(a)). These contour lines track the propagation of the perturbations causing a zero-displacement at the x-direction. Note that, the propagation curves contain small noises caused by the motion of magnetic islands. We fit the contour lines by the power function to clean the noises and then take their time derivatives to obtain the propagation speeds, namely, V i p (y). The speeds calculated from six different tracks are very similar and we take their averaged value V i p to evaluate their mean feature (see Figure 5(a)). The speed increases from 0.89 to 2.17 with height, corresponding to the real values from 388 to 946 km/s. Furthermore, we calculate the local Alfvén speed at the CS outer boundary (V Acs ). It is found that the variations of V i p with height are largely similar to the profile of the time-average Alfvén speed V Acs (see Figure 5(a)). Considering the similarity between V Acs and V i p , we propose a novel way to determine the magnetic field strength of the CS outer boundary. Observationally, one can use the averaged propagation speeds of several CS oscillations, V i p , as a proxy of the local Alfvén speed. With a combination of the plasma density estimated by differential emission measure (DEM), the magnetic field strength of the CS boundary is then derived as Here, we verify this method using the numerical data ( Figure 5). The time-averaged value ρ cse are used to simulate the density derived by the DEM analysis (Figure 5(b1) and (b2)). We find that the profiles of B i p we evaluated are basically in agreement with the distribution of B Acs in height (see Figure 5(c)).

SUMMARY AND DISCUSSION
In this letter, we study the formation mechanism of the CS oscillations using MHD simulation. The fully-developed reconnection presents a quasi-periodic feature. Various downflow magnetic islands collide with the flare loop top and result in a turbulent cusp-shaped structure. Following the collision of a relatively large island at the loop top, fewer and smaller islands are generated. The KHI is then switched on in the interface of two sheared flows, from the persistent reconnection downflows. The asymmetric flows, caused by the asymmetric KHI, finally drive large-amplitude swaying of the flare loop-top and thus the oscillations of the CS.
The CS oscillations, though damping with distance, are still detectable when they propagate to a higher height. More importantly, the CS oscillations are proven to propagate with the local Alfvén speed. It thus can be used for evaluating the magnetic field strength of the outer boundary of the CS once we derive the local plasma density, for example, through the DEM technique. This method is well verified by our numerical data and is hopeful to be comparable with other methods (e.g., Chen et al. 2011;Nakariakov & Ofman 2001;Liu et al. 2011;Tian et al. 2012;Yang et al. 2020).
However, note that this method only provides a zero-order approximation of magnetic field strength. Observationally, the main challenge is how to accurately measure the propagation speed of the CS oscillations, which are largely limited by observational tempo-spatial resolution. Our result shows that the magnetic field at the outer boundary of the CS increases with the height, which is different from real observations and recent simulations of the flux rope eruption that use decaying magnetic fields in height (e.g., Chen et al. 2020). The main reason is that we take advantage of a uniform background magnetic field along y-direction and the appearance of upward moving magnetic islands during the CS oscillations.
Differing from Takahashi et al. (2017), in which the oscillations above the flare loops were also reproduced but for the TS fronts, here, we mainly focus on the fundamental mechanism and properties of large-amplitude oscillations of the whole CS. The mechanism for the CS oscillations revealed here is different from that in Xie et al. (2021) who considered the effects of CME but did not include thermal conduction, radiation cooling, and background heating. They contributed the CS oscillations to the RTI at the bottom of CME. Although excluding the flux rope eruption, our model provides much more fine structures at the loop top region. We find that, the positions where the RTI may be initiated are spatially scattered and temporally varied and thus, compared with the RTI, the KHI plays a more important role in giving rise to the CS oscillations.

ACKNOWLEDGMENTS
We would like to thank the anonymous referee for valuable suggestions. This research is supported by the Natural Science Foundation of China grants 11722325, 11733003, 11790303, and 11790300.

A. THE CS STRUCTURE
We take four typical slits across the CS at t = 9.1 to exhibits the detailed CS structures (see Figure A1(e)). The red dashed lines mark the position with the maximum temperature gradients (|∂T /∂x|) on both sides of the CS, which approximately correspond to the discontinuity fronts of the slow-mode shocks (also see Mei et al. 2012). We take their average x-coordinate as the central position x cs (y) of the CS (see the blue dashed lines in Figure A1).
Between the slow shocks, the polarity of magnetic field reverses rapidly in space (see Figures A1(d1)-(d4)). However, differing from standard CS models (e.g., the Harris-sheet or the Petschek-sheet), the CS structure shows complex details. Near the principal X-point (y = 1.1), the CS is thinnest and shows a Gauss-type distribution. At y = 0.8 and y = 1.8, on both sides of the CS, the slow-shock fronts cause rapid decreases of magnetic field strength and thus induce two localized strong currents, between which a lower current plateau forms. J z near x = 0 can be either negative (Figure A1(c1), also see Mei et al. (2012)) or positive (Figure A1(c4)). Near the center of a magnetic island (y = 1.356), the density and current increase significantly (Figure A1(a3) and (c3)), while the temperature and magnetic field profiles are similar with other slits (Figure A1(b3) and (d3)).
For all four positions, however, one can see that the region of the CS is well enclosed by the red dashed lines and is approximately symmetric about the blue dashed lines. Therefore, we determine the central position of the CS by the blue dashed line.