The Dynamical Behavior of Reconnection-driven Termination Shocks in Solar Flares: Magnetohydrodynamic Simulations
Chengcai Shen, Xiangliang Kong, Fan Guo, John C. Raymond, Bin Chen
DDraft version December 5, 2018
Typeset using L A TEX manuscript style in AASTeX61
THE DYNAMICAL BEHAVIOR OF RECONNECTION-DRIVEN TERMINATION SHOCKS INSOLAR FLARES: MAGNETOHYDRODYNAMIC SIMULATIONS
Chengcai Shen, Xiangliang Kong, Fan Guo, John C. Raymond, and Bin Chen Harvard-Smithsonian Center for Astrophysics, 60, Garden Street, Cambridge, MA, 02138, USA Shandong Provincial Key Laboratory of Optical Astronomy and Solar-Terrestrial Environment, and Institute of SpaceSciences, Shandong University, Weihai, Shandong 264209, China Los Alamos National Laboratory, P.O. Box 1663, Los Alamos, NM 87545, USA Center for Solar-Terrestrial Research, New Jersey Institute of Technology, 323 Dr. Martin Luther King Blvd, Newark,NJ 07102, USA
ABSTRACTIn eruptive solar flares, termination shocks (TSs), formed when high-speed reconnection outflowscollide with closed dense flaring loops, are believed to be one of the possible candidates for plasmaheating and particle acceleration. In this work, we perform resistive magnetohydrodynamic simula-tions in a classic Kopp-Pneuman flare configuration to study the formation and evolution of TSs, andanalyze in detail the dynamic features of TSs and variations of the shock strength in space and time.This research focuses on the fast reconnection phase when plasmoids form and produce small-scalestructures inside the flare current sheet. It is found that the TS emerges once the downward outflowcolliding with closed magnetic loops becomes super-magnetosonic, and immediately becomes highlydynamical. The morphology of a TS can be flat, oblique, or curved depending on the detailed inter-actions between the outflows/plasmoids and the highly dynamic plasma in the looptop region. TheTS becomes weaker when a plasmoid is crossing through, or may even be destroyed by well developedplasmoids and then re-constructed above the plasmoids. We also perform detailed statistical analysison important physical quantities along and across the shock front. The density and temperatureratios range from 1 to 3 across the TS front, and the pressure ratio typically has larger values up to10. We show that weak guide fields do not strongly affect the Mach number and compression ratios,and the TS length becomes slightly larger in the case with thermal conduction. a r X i v : . [ a s t r o - ph . S R ] D ec Shen et al.
Keywords: magnetic reconnection — magnetohydrodynamics (MHD) — shock waves— Sun: flares
ASTEX Dynamical TS INTRODUCTIONIt has been widely accepted that a huge amount of magnetic energy (up to 10 erg) is violentlyreleased via magnetic reconnection in a typical solar eruption. The released magnetic energy isquickly converted into bulk plasma kinetic energy, plasma thermal energy, and energy contained inaccelerated energetic particles. However, the exact mechanism for accelerating the charged particlesremains unclear. There are several competing mechanisms for explaining the acceleration of chargedparticles in solar flares, including magnetic reconnection electric currents, turbulence/waves, andfast-mode shocks (see Miller et al. 1997; Zharkova et al. 2011, and references therein).A termination shock (TS) has been predicted in early solar flare models. In the standard flaremodel, also known as the CSHKP model (Carmichael 1964; Sturrock 1968; Hirayama 1974; Kopp &Pneuman 1976), fast reconnection outflow jets are produced as a result of magnetic reconnection.They collide with the newly closed magnetic loops and may form a fast-mode TS, if the outflowspeed exceeds the local fast-magnetosonic speed in the looptop region (Forbes & Malherbe 1986).The flare TS has been suggested as a possible particle accelerator by many authors (Masuda et al.1994; Shibata et al. 1995; Forbes & Acton 1996; Somov & Kosugi 1997; Tsuneta & Naito 1998; Mannet al. 2009; Guo et al. 2012; Kong et al. 2013; Li et al. 2013), and is widely adopted in standard flaremodel cartoons including the arguably most famous one by Shibata et al. (1995). It may also play animportant role in the heating of plasma in or above the post-reconnection flare loops (e.g., Masudaet al. 1994; Guidoni et al. 2015) .The TSs also have been investigated in numerical studies (e.g., Forbes & Malherbe 1986; Forbes1988; Workman et al. 2011; Takasao et al. 2015). A stationary fast shock is identified in magnetohy-drodynamic (MHD) numerical experiments with line-tied magnetic reconnection when the downward-directed reconnection jet encounters the obstacle formed by the closed loops (Forbes 1986). In general,the shock normal of the flare TS is nearly perpendicular to the magnetic field (i.e., the angle betweenthe upstream magnetic field and shock normal vector θ Bn > ◦ ). Forbes (1986) found the existenceof this shock with a compression ratio of 2.0 and a Mach number as high as 2.3. The numericalsimulations by Workman et al. (2011) have shown similar results. Forbes (1986) also pointed out Shen et al. that the transition from the supermagnetosonic flow region upstream of the shock to the nearly staticdownstream is complicated, and a deflection sheath downstream of the fast shock is necessary for theformation of fast shocks according to the MHD jump conditions. More recently, structures of TSshave received more attention in MHD simulations by Takasao et al. (2015) and Takasao & Shibata(2016). By including essential physics for solar flares such as magnetic reconnection, heat conduction,and chromospheric evaporation, their numerical model revealed that flare loops and the above-the-loop-top region are filled with various shocks and waves. They reported multiple shocks, includinghorizontal and oblique shocks above the looptop, and found the quasi-periodically oscillations in theabove-the-loop-top regions. This suggests that the structure of TSs is more complex than previouslyassumed and could significantly affect energetic electron acceleration and plasma heating, and theassociated observational signatures.Although TSs are often invoked in the standard solar eruption models, there are few solid observa-tional constraints because they are difficult to observe. One piece of evidence for TSs is the looptopradio emission, which shows spectroscopic features similar to solar type II radio bursts (an indicationof shocks in the corona) but with small frequency drift rates (Aurass et al. 2002; Aurass & Mann2004). If these radio bursts are associated with plasma radiation for which the emission frequency f ∝ n / e (where n e is the plasma density), a small frequency drift rate suggests a slow overall changeof plasma density during the shock evolution, which implies a quasi-standing shock wave locatedin the looptop region. Recently, using the Karl G. Jansky Very Large Array (VLA), Chen et al.(2015) presented radio spectroscopic imaging of an eruptive solar flare event. Their imaging resultsshowed that a TS, which appeared in the radio dynamic spectrum as a slow-drift type-II-burst-likefeature, is located in the looptop region in front of fast reconnection outflows. The correspondingRHESSI observation shows that the HXR loop-top source is nearly co-spatial with, but slightly be-low the shock front delineated by the radio source centroids made at different frequencies, agreeingwell with the scenario that the TS is responsible for accelerating high-energy electrons and producesHXR emission in the shock downstream region. Furthermore, radio spectroscopic observations havesuggested that the shock is unsteady in nature based on small but nonzero frequency drift in the ASTEX Dynamical TS MODELS2.1.
Overview
Our MHD model follows the classic two-ribbon flare configuration in two dimensions. It has beenused to study the evolution of small-scale structures (e.g., plasmoids) inside the reconnecting current
Shen et al. sheet, the evolution of magnetic loop and magnetic arch structures (Shen et al., 2011, 2013a), andthe dynamic features of TSs above flare loops (Chen et al., 2015). The simulation starts with acurrent sheet in mechanical and thermal equilibrium that separates two regions of magnetic fieldwith opposite polarity. By adding a small perturbation on the initial equilibrium current sheet,the system commences to evolve. Magnetic reconnection is initiated gradually from the perturbationregion, which produces two reconnection outflows and forms an arcade of newly-reconnected magneticloops anchored at the bottom boundary where the magnetic field has been set to be line-tied on thephotosphere.The governing resistive MHD equations are as following: ∂ρ∂t + ∇ · ( ρ v ) = 0 , (1) ∂ρ v ∂t + ∇ · ( ρ vv − BB + P ∗ ) = 0 , (2) ∂ B ∂t − ∇ × ( v × B ) = η m ∇ B , (3) ∂E∂t + ∇ · [( E + P ∗ ) v − B ( B · v )] = S, (4)where P ∗ is a diagonal tensor with components P ∗ = P + B / E isthe total energy density E = Pγ − ρv + B , (5) γ = 5 / S = µ η m j + ∇ || · κ ∇ || T , which includesOhmic dissipation and thermal conduction. Here µ , η m and κ are the magnetic permeability of freespace, magnetic diffusivity, and parallel component of Spitzer thermal conduction tensor.We normalize above MHD equations to non-dimensional forms using characteristic values in ourcalculations: B ∗ = B /B , ρ ∗ = ρ/ρ , p ∗ = p/ ( B /µ ) , v ∗ = v /V , T ∗ = ( β / T /T ), and J ∗ = J /J with V = B / √ µ ρ and β = 2 µ p /B . Here an ambient gas pressure p depends on β whichwill be assigned in the following parameter Table 1. We set L = 7 . × km, B = 0 .
004 Tesla, ρ = 1 . × − kg/m , T = 2 MK, t = 92 . V = 812 . J = 4 . × − Am − tomodel a typical solar flare in following simulations. ASTEX Dynamical TS
Code description
We use the Athena code to numerically solve the above MHD equations. Athena is a publiclyavailable grid-based code for astrophysical MHD (Stone et al. 2008). In this work, we include Ohmicresistivity in MHD equations, and omit gravity and optically thin cooling. The code is based onthe directionally unsplit high order Godunov method, and combines the corner transport upwind(CTU) and constrained transport (CT) methods. It provides superior performance for capturingshocks as well as contact and rotational discontinuities. The low numerical diffusion feature ofAthena code provides a significant advantage for resistive MHD simulations of processes like magneticreconnection. 2.3.
Initial and boundary conditions
The initial configuration includes a pre-existing Harris-type current sheet along y -direction withthe non-dimensional width w = 0 . B x ( x, y ) = 0 , (6) B y ( x, y ) = tanh ( x/w ) , (7) p ( x, y ) = (1 + β − B y ( x, y ) ) / , (8) ρ ( x, y ) = 2 p ( x, y ) /β . (9)In order to have the system evolve rapidly from the initial steady state to a bursty reconnectionphase, we introduce perturbation magnetic field B x and B y on this pre-existing current sheet asfollows: B x = 2 πL y A pert cos ( πxL x ) sin ( 2 π ( y − y c ) L y ) B , (10) B y = − πL x A pert sin ( πxL x ) cos ( 2 π ( y − y c ) L y ) B . (11)Here A pert = 0 .
01 is the non-dimensional perturbation strength, y c = 0 .
25 is the perturbation centeralong the y − axis. L x and L y are non-dimensional perturbation wavelength which are set to 2 and 3.5in order to minimize perturbations at boundaries. The initial current density, magnetic configurationand pressure distribution are plotted in Figure 1. Shen et al.
The boundary conditions are arranged in this following way. The right ( x = 1 . x = − . y = 2 .
0) of the simulation domain are open or free boundaries on which the plasmaand the magnetic flux are allowed to enter or exit freely, and the boundary at the bottom ( y = 0 . ∂B y ( x, , t ) ∂t = 0 , (12)and the plasma does not slip and is fixed to the bottom boundary, v ( x, , t ) = 0 , (13)respectively. For gas pressure and plasma density, we set ∂P ( x, y, t ) ∂y | y =0 = 0 , (14)and ∂ρ ( x, y, t ) ∂y | y =0 = 0 . (15)We also set J z vanish at the bottom boundary as well, namely J z = ( ∂B y ∂x − ∂B x ∂y ) | y =0 = 0 for t > . Then B x ( x, , t ) can be set according to the above condition.2.4. Parameter table
The parameters of our numerical simulation cases are listed in Table 1. In Cases 1 and 2, we fixthe ambient magnetic field strength to B and set different gas pressure to change the backgroundplasma β = p / ( B / µ ). In Cases 3, 4, and 5, we introduce a uniform B z in the whole calculationdomain to study the effect of magnetic guide field on TSs. In Case 6, we include anisotropic thermalconduction to compare with Case 7. We also perform simulations from different initial current sheetsin Cases 6 ∼ current sheet. Driven by same perturbations, these narrow B y ( x, y ) = sin ( πx/ w ) when | x | < = w , and B y ( x, y ) = sign ( x ) when | x | > w (see Shen et al. 2011, in detail). ASTEX Dynamical TS η m on the whole computational domain, which gives a constant magnetic Reynolds number R m = 10 . RESULTSDriven by the initial perturbation on magnetic fields defined in Equations (10) and (11), the initialcurrent sheet becomes thinner and more intense, and eventually develops magnetic islands. In thefollowing analysis, we focus on Case 1, and show the common features of dynamical terminationshocks. In section 3.1-3.3 we present detailed analysis on the properties and dynamic features of
Table 1.
Parameters and setup for all simulation casesCases Grids Background β Constant B z Thermal Conduction1 2048 2048 0 . .
02 0 No3 2048 2048 0 . . . . . Note —Magnetic Reynolds number R m = 10 . Cases 6 ∼ Shen et al.
Figure 1.
The evolution of current-density magnitude | J | , y component of velocity V y and pressure P . Atthe early phase ( t = 0 , t ), magnetic reconnection steadily takes place. After t = 90 t , rapid reconnectionappears as plasmoids develop inside the current sheet. Later, the reconnection becomes violent and morefine structures appear in this unsteady phase. the termination shock. In section 3.4, we will discuss effects of magnetic guide field, and presentsimulations including classical Spitzer thermal conduction (Spitzer 1962) in section 3.5. ASTEX Dynamical TS | J | , velocity in the y -direction V y , and pressure P . In the beginning, the magnetic field starts to diffuse at the X-point whereperturbation is introduced. The magnetic fields at two sides of the current sheet slowly move towardone another due to the Lorentz force attraction between the field lines of opposite polarity. A pair ofmagnetic outflows gradually form and closed magnetic loops accumulate below the reconnection sitedue to the reconnected magnetic flux. This process gradually squeezes the initial current sheet, andthe sheet narrows until it becomes intensified ( | J | is higher than 300 J after time 90 t ) and thinenough that the tearing mode grows and becomes nonlinear (e.g., see Furth et al. 1963; Loureiroet al. 2007; Ni et al. 2010). At this time, many magnetic islands appear inside the sheet andreconnection suddenly becomes violent. The reconnection outflow speed is significantly enhanced,reaching a maximum speed of 0.7 V . The fast reconnection results in the rapid growth of closedmagnetic loops, which would be observed as hot flaring loops (see Figure 1).We find that the TS develops when the downward plasma flow has a speed exceeding the localmaximum fast-magnetosonic speed and it is stopped by the magnetic loop. As shown in Figure 2, thereconnection downward outflow rapidly increases after time 89 t , and the outflow speed first exceedsthe local fast magnetosonic speed along the reconnecting current sheet. During this period, thesignificantly enhanced current density indicates the reconnection becomes more efficient and violentas well. We also monitor properties in flare loop-top regions and show how average magnetic fieldstrength, average density and average fast magnetosonic speed change with time in panels ( d ) − ( f ).It is clear that amount of magnetic flux and density has been accumulated as the rapid reconnectionwas taking place, and a strong magnetic loop has been formed after 89 t and becomes more denseafter 94 t . After 95 t , the loop-top density decreases as the inflow plasma gets thinner during laterphases of the fast magnetic reconnection. This relatively steady magnetic loop structure is importantfor the formation of TSs, as it becomes an obstacle that stops the downward super-magnetosonicoutflows. Note we do not include chromospheric evaporation in the current study, which wouldgreatly enhance the density of the post-flare loops and possibly, facilitate the shock formation as thelocal Alfv´en speed is inversely proportional to √ ρ .2 Shen et al. )0.250.500.75 (b) Maximum Downflow (V )12 (c) Maximum Mach number M>1 )345 (e) Loop-top ( )60 70 80 90 100Time (t )0.40.50.6 (f) Loop-top Magnetosonic Speed (V ) Figure 2.
The maximum current-density, downward outflow speed and Mach number ( M = v/ ( C s + V a ) / ) and average magnetic field strength, average density and average fast magnetosonic speedin loop-top regions within 0.1 L heights below the end point of downward outflows. Here v is plasma flowspeed, C s is the local sound speed, and V a is Alfv´en speed. The dashed black line indicates the time 89 t when the downward flow first exceeds the local fast magnetosonic speed. ASTEX Dynamical TS
Jump condition and Morphology
We first analyze the velocity distribution above the flare magnetic loops to identify the exactlocation of the TS front. The plasma velocity should be super-fast-magnetosonic on the upstream sideand become sub-fast-magnetosonic on the post-shock (or downstream) side, respectively. Therefore,we plot the Mach number defined by v/ (cid:112) C s + V a , the ratio of plasma velocity over the maximumfast-magnetosonic speed perpendicular to B in Figure 3(a) to show the surface of TS front. It is clearthat the downward reconnection flow Mach number steps from above 1.4 to less than 1.0 crossinga sharp boundary, which is the TS front. The TS front can also be easily identified by examiningthe velocity divergence ( ∇ · v ) map (Figure 3(b)): it appears as a sharp layer with negative values of ∇ · v , as strong compression is expected in the immediate vicinity of the shock. In Figure 3(b), theshock front is delineated in the ∇ · v map as a thin, red curve, and its location matches well with thesharp discontinuity in the Mach number map in Figure 3(a).We then represent primary variables in the co-moving frame with the TS front and compute thefast mode Mach number ( M F ). Based on the ∇ · v map, we can fit the TS front profile using 8 orderpolynomial fitting as shown by dotted line in Figure 3(b), and obtain the spatial position of the TSfront. We then estimate the TS normal velocity ( v T Sn ) using backward difference in time betweencurrent and previous frames. For each point on the TS front, we find its previous position (e.g., t = 95 . t ) along the TS normal direction, measure moving distance of the chosen point during thisperiod, and obtain v T Sn (see Figure 3(c)). At each position along the TS front curve, we choose a setof sampling cuts along the local TS normal direction as shown in Figure 3(b), and obtain the normalcomponent of primary variables along these cutting lines. Then the fast mode Mach number can becalculated by: M F = | v T Sn − v n | /C F . Here v T Sn is the normal TS velocity at each sampling cut, v n isthe normal velocity of plasma flows to the TS, and C F = [ ( C s + V a + (cid:113) ( C s + V a ) − C s V a cosθ Bn ] / is the local magnetosonic fast-mode wave speed, and θ Bn is the angle between the local magnetic fieldand the shock normal. We show distributions of M F along these chosen sampling cutting lines inFigure 3(d). Here the horizontal axis is for x − location of the TS, and the perpendicular axis is thedistance away from the TS front along each of cutting lines. It can be seen that M F sharply jumps4 Shen et al. between upstream and downstream regions. For example, along the green cutting line, M F is as highas ∼ ∼ ∇ · v map. Therefore, we can see the fast mode Mach number reduces to oneat the blue and orange cuts, and even less than one outside these two lines. We then can estimatethe expansion length of the TS front in x − direction, which ranges from ∼ − .
023 to ∼ .
02 in thecondition M F > ∇ · v ) and fast-mode Mach number ( M F ), (b) gaspressure ( p ), density( ρ ) and temperature( T ), (c) normal component B n and transverse component B t of magnetic field, (d) velocity components v n and v t , (e) mass flux component ρv n and magneticflux component B t v n − B n v t , (f) momentum flux p + B / − B n + ρv n and energy flux ( p + B / v n − B n ( B · v ) + ( ρe + ρv / B / v n .The transition of pressure between upstream and downstream regions is very sharp, and jumpsfrom 0.1 to around 0.38. The location of the shock front also can be identified by the minimum ∇ · v , which matches well with the pressure jump. Similar sharp jumps can also be found in density,temperature, and the transverse magnetic component B t and normal velocity competent V n profilesin Figure 4(b) ∼ (d). It is worth noting that we do not resolve the real dissipation scale of the shockfront (e.g., the order of ion skin depth) in these MHD numerical experiments because that scale isextremely narrow in high temperature plasmas such as the solar corona. Nevertheless, the transitionstill can be well-approximated as a discontinuous change between two sides of the shock front. Aswe can observe in Figure 4(a), the discontinuity is not infinite thin, with a thickness of roughly 2 ∼ B n (blue curve in Figure 4(c)) is small relative tothe transverse component B t , but it does not vanish either on the upstream or downstream sides.A horizontal flow v t , albeit very small, can also be clearly seen (orange curve in Figure 4(d)). That ASTEX Dynamical TS )0.580.600.620.640.660.68 y ( L ) (a) Time = 96.00 (t ) V / V A + C s )0.580.600.620.640.660.68 y ( L ) (b) Fitted TS frontTS normal= 30.5°TS normal= 92.6°TS normal=132.8° -156-117-78-3903978117156 V ( t ) )0.580.600.620.640.660.68 T S H e i g h t a t y - a x i s ( L ) (c) TS front at 95.9 t TS Front at 96.0 t V TSn at 96.0 t T S N o r m a l V e l o c i t y ( v ) )0.0080.0060.0040.0020.0000.0020.0040.0060.008 D i s t a n c e a c r o ss T S f r o n t ( L ) (d) DownstreamUpstreamFitted TS frontTS normal= 30.5°TS normal= 92.6°TS normal=132.8° 0.00.20.40.60.81.01.21.41.61.82.0 F a s t m o d e M a c h nu m b e r ( M F ) Figure 3. (a) Distribution of Mach number v/ ( C s + V a ) / , (b) Velocity divergence ∇ · v at time t = 96 t .The black contour lines in panels (a) and (b) show magnetic field lines, and the red curve with the minimum ∇ · v indicates the position of TS front in panel (b). Three straight cuts (blue, green and orange lines) showthe local TS normal directions in panel (b). Panel (c) shows TS heights at times 95.9 and 96.0 t , and thenormal velocity V T Sn at t = 96 t . Panel (d) shows the results across and along the TS at t = 96 t . Thehorizontal axis and three sampling cuts (blue, green and orange lines) are same as in panel (b). indicates the TS is nearly perpendicular in the particular position we picked (also see green line inFigure 3(b)), and the shock is dynamically evolving with weak horizontal flows.Along the cut, the momentum component and magnetic flux component are roughly constantbetween upstream and downstream regions. As can be seen in Figure 4(e), the relative changes for ρv n and B t V n − B n V t are roughly 1% and less than ∼
2% between the upper edge and lower edge ofthe TS region (gray shaded regions in Figure 4), respectively. The momentum flux and energy flux6
Shen et al. V a l o n g s a m p li n g li n e ( t ) (a) TS frontTime = 96.00 (t ) F a s t M o d e M a c h N u m b e r ( M F ) N o r m a li z e d P , , T (b) P (0.1p ) (1.4 )T (0.1 T )1.8e-031.9e-032.0e-032.1e-032.2e-032.3e-032.4e-032.5e-032.6e-03 - B n ( B ) (c) - B t ( B ) - V n ( V ) (d) - V t ( V ) )-1.00-0.98-0.96-0.94-0.92-0.90-0.88-0.86 v n ( . e - v ) (e) B t v n - B n v t ( . e - B v ) )1.001.021.041.061.081.10 p + B / - B n + v n ( . e - p ) (f) -1.100-1.080-1.060-1.040-1.020-1.000 ( p + B ) v n B n ( B v ) + ( e + v + B ) v n ( . e - p ) Figure 4.
Variation of primary variables across the shock front along the TS normal direction at x = 0in the co-moving frame: (a) ∇ · v and M F ; (b)pressure p , density ρ and temperature T ; (c) magnetic fieldcomponents B n and B t ; (d) velocity components v n and v t ; (e) mass flux and magnetic flux; and (f)momentum flux and energy flux. The red and blue shades indicate pre-shocked and post-shocked regions,respectively. The gray shade is for the TS front. parallel components are also conserved, and the relative changes between upstream and downstreamare still less than ∼
2% as shown in panel (f). It is noticed that there are unavoidable measurementerrors when estimating TS raising velocity and computing TS normal direction, which causes slight
ASTEX Dynamical TS y − axis in this paper) where the magnetic field onlyhas a horizontal component. In other words, the magnetic field direction both in the upstream anddownstream sides is parallel to the TS front at the symmetry center. However, the magnetic con-figuration could be asymmetric in many realistic environments, which may cause variations on theTS geometry. Instead of investigating the asymmetric configuration itself, we show here variableTS slopes when the system evolves to the asymmetrical case due to the numerical perturbations.In Figure 5, we display the shape of the TS illustrated by ∇ · v at four different times. The shockfronts show time-varying shapes including quasi-flat, horizontal and oblique components, curves, andoblique-flat shapes. In panel (b), a pair of sub-shocks can be seen below the horizontal TS front attime 96.30 t . It could appear as the reflection of plasma flows becomes strong below the TS, anddisappear once there are not reflected flows in the downstream regions. In the following analysis, wefocus on the upper shock structure which is directly driven by the reconnection downflow, if multiplesub-shock structures are present at the same x coordinate.3.2. Dynamics of TSs Shen et al.
Figure 5.
Divergence of fluid velocity at four different time frames showing the shapes of TS front. Thedotted lines indicate the edge of the shock front where primary variables jump crossing shock front. The TSshape is strongly modified by the intermittent reconnection outflow and loop-top structure.
In Figure 6(a), we plot the distribution of the flow Mach number ( v/ (cid:112) C s + V a >
1) along the y axis ( x = 0), which is the direction of the reconnection current sheet, as a function of time between t = 95 t and 100 t ( ∼ ∇ · v isminimum. During this period, the reconnection outflows are accelerated away from the X-point, andbecome super-magnetosonic in the exhaust region until they hit the top of the closed flaring loops.The acceleration of upward and downward moving plasmoids away from the X-point nearby y = 1 . v/ (cid:112) C s + V a > . < y < .
75. On the upstream side of the shock,
ASTEX Dynamical TS β close to the TS appears to vary with time. The local plasma β is clearly closeto the values in magnetic reconnection downflows, and significantly larger than the plasma β insideflare loops. As shown in Figure 6(b), the lower range of the local plasma β at the TS is around 100.This confirms that the TS forms above magnetic loops and the shock properties are not stronglyaffected by the plasma inside magnetic loops. y ( L ) (a) Minimum V V / C s + V a
95 96 97 98 99 100Time (t )0.00.51.01.52.0 y ( L ) (b) Minimum V P l a s m a Figure 6.
Distributions of (a) flow Mach number v/ (cid:112) C s + V a and (b) plasma β along the y-axis ( x = 0)as functions of time. The dashed lines show minimum ∇ · v , which indicate the height of TS fronts. In addition, we calculate the compression ratio between upstream and downstream, and show herehow the compression ratio changes spatially and temporally. In Figure 7, we plot gas pressure ratio,density ratio, temperature ratio, and M F along the shock front. We trace the TS front accordingto the minimum ∇ · v at each time, and obtain a set of sampling points at the TS front. For eachsampling point, we then calculate primary variables on upstream/downstream sides by 2D spatialinterpolation along the normal direction of TS front at this point. In Figure 7, the vertical axis shows0 Shen et al. the location of sampling points in the x-direction, and the horizontal axis is for time. During theplasmoid reconnection phase (95-100 t ), the density compression ratio intermittently increases dueto unsteady reconnection downward flows. The dominant density ratio ranges from 1 to 3 dependingon time and location. The maximum density ratio can be as high as 3.7 near a corner of two obliqueshock fronts where the downstream/upstream density could be slightly overestimated/underestimateddue to interpolation process. The intermittent nature of the quantities for both spatial distributionand temporal distribution is clear. The variation of pressure and temperature ratios are also clear,and it is also easy to see that maximum value can larger than 10 and 3, respectively. The distributionof fast mode Mach number M F also varies in space and time, and ranges from 1 ∼ / ( γ − / (Soward & Priest 1982; Forbes 1986) for Petschekreconnection with an inflow plasma of zero β . In the case of γ = 5 /
3, this gives M F ≈ .
73. In fact,the overall M F in our simulations agrees in general with their expected values as shown in Figure7(d). Furthermore, M F in our simulation also can occasionally exceed the predicted value due tothe occasion ”bursts” of the outflow speed, which is associated with the development of tearing (orplasmoid) instabilities.Along the shock front, the density compression dramatically varies with time. The spatial locationof the maximum density compression ratio is not always at the point directly below the currentsheet (i.e., x = 0), but varies in time from x = − . L to x = +0 . L . The highly variablenature of the density compression along the shock front is associated with the slope of the TSfront, the bursty reconnecting downward flows, and/or the detailed properties of plasmoids. Anotherinteresting feature is that the spatial position where maximum temperature ratio appears to departfrom either the position of maximum pressure ratio or density ratio at some particular time. Forexample, the maximum temperature jump appears at around x > . t , while the maximumdensity compression appears on another side x < x ∼
0) as well. This is not surprising since the magnetic field configuration
ASTEX Dynamical TS x ( L ) (a) p p o s t / p p r e x ( L ) (b) p o s t / p r e x ( L ) (c) T p o s t / T p r e
95 96 97 98 99 100Time (t )0.040.020.000.020.04 x ( L ) (d) F a s t m o d e M a c h N u m b e r M F Figure 7.
The compression ratio distribution across the TS front at different times. (a) Gas pressure ratiobetween post-shock and pre-shock sides; (b) and (c) are density ratio and temperature ratio; (d) is the fastmode Mach number M F . The threshold condition M F > . Shen et al.
In Figure 8, we display the normal direction of the shock front and examine the shock propertyduring its dynamical evolution. We measure the shock normal direction and magnetic field directionat each point along the shock surface, and show them in Figure 8(a) and (b). We plot the shock-normal angle, measured from the normal direction to the x-axis, and calculate the interior anglebetween TS normal direction and magnetic field direction in panel (c). It can be seen in panels (a)and (b) that the TS direction dynamically changes during time 95 t to 100 t . At the system center( x = 0), the shock front direction ranges from around 60 degrees to around 130 degrees. In panel(c), we can see that the interior angle is close to 90 degrees at the system center ( x = 0), wherethe TS can be regarded as a nearly perpendicular shock for most of the time. However, as shown inFigure 8(c), it is by no means a stable perpendicular shock, as usually depicted in the standard flarecartoons (references to e.g., Shibata et al. 1995), even at x = 0. Moreover, the TS becomes moreand more oblique towards the edge of shocked surface, where the interior angle gradually decreasesto <
30 degrees.It is convenient to obtain the TS length measured along the TS front profile on the ∇ · v map. InFigure 9, we trace shock front profiles extending to the left and right sides from the system center( x = 0), and define the end points according to the fast mode Mach number. The threshold value atthe end points are chosen either as 1.0 or 1.5. In Figure 9, the cyan line is for the threshold M F > . M F > .
5. Because the evolution of thelength is highly dynamical, the maximum length is larger than 0.09 L while the minimum length canbe as short as around 0.01 L at a particular time. Scaling it according to the characteristic length( L = 7 . × km) in Section 2, the TS length is in range of < . ∼ Effect of Plasmoids on Mach Number
As described above, we have shown that the collision between plasmoids and the closed magneticloops can cause dynamic evolution of the TS. The shock front can temporally change its height andshape. In this section, we look in detail at the shock strength when the collision takes place. Wechoose two typical episodes in the following analysis from Case 2 in which the background magnetic
ASTEX Dynamical TS x ( L ) (a) T S N o r m a l T S ( d e g ) x ( L ) (b) B D i r e c t i o n ( d e g )
95 96 97 98 99 100Time (t )0.020.000.02 x ( L ) (c) D i ff A n g l e ( d e g ) Figure 8. (a) Normal direction of the TS front at different times. The angles are measured from the x − axis direction to the normal direction; (b) Magnetic field direction at the shock front; (c) Interior anglebetween normal direction and magnetic field direction. The TS front is defined as in Figure 7. field is relatively stronger due to the initial lower plasma β . In this case, we can see more small scalemagnetic islands during the evolution, and choose two typical events to analyze the TS strength indetail.In Figure 10, we demonstrate how the Mach number of the TS changes when a small down-ward moving blob is merging into the closed magnetic loops. We color the downward flows with4 Shen et al.
95 96 97 98 99 100Time (t )0.010.020.030.040.050.060.070.080.09 L e n g t h o f T S F r o n t ( L ) M F > 1.0 M F > 1.5 Figure 9.
The length of the TS front. The cyan line is for threshold of M F > .
0, and orange line is for M F > . v/ (cid:112) C s + V A > y − axis in panel(d) for different times. At time 94.8 t , the shock front is located at the height y ∼ . L where themaximum Mach number is larger than 2.5. Meanwhile the downward moving magnetic island hasarrived to height y ∼ . L . We can see a lower Mach number at the core of this magnetic islanddue to higher pressure and magnetic field than the ambient downward flows. At time 94.9 t , thecenter of this magnetic island arrives at the shock front. In the Mach number plot, it manifests itselfas a slightly lower Mach number valley near the shock front. After this magnetic island has beentotally merged into the previous closed magnetic loops, the TS front appears at a higher altitude, y ∼ . L . We can see more clearly the evolution history for this collision process in panel (d). Inthis panel, the dashed line is for time 94.9 t when the magnetic island is crossing the shock front. Itis clear that the fast mode Mach number on the upstream side quickly decreases to around 1.5 fromabove 2.5. Once the magnetic island has merged into the magnetic loops, the Mach number risesagain to 2.2 as shown by the dark green line as the shock moved up to a higher altitude.Another collision event is shown in Figure 11. In this case, the incoming plasmoid has fullydeveloped for a relatively long time inside the current sheet and grows to a larger size. As this ASTEX Dynamical TS y ( L ) (a)t= 94.80 (t ) V / C s + V A Arriving Island V y ( L ) (b)t= 94.90 (t ) Arriving Island )0.500.550.60 y ( L ) (c)t= 95.10 (t ) 0.05 0.00 0.05x (L )1.0 1.9 2.7 -203 0 2030.4495 0.5000 0.5505 0.6010 0.6515y (L )0.00.51.01.52.02.5 V / C s + V A a l o n g y - a x i s (d)94.8t Figure 10.
Impacts of a small plasmoid on the termination shock. plasmoid accumulated magnetic flux and mass, it moved more slowly than the ambient outflow. Inpanels (a) and (b), we can see the plasma surrounding this magnetic island is sub-magnetosonic attimes 91.8 t and 91.9 t . At the same time, the TS front is at y ∼ .
54. Once this magnetic islandreaches the shock front after time 92.0 t , the TS front completely disappears, and the center region6 Shen et al. y ( L ) (a)t=91.8(t ) (b)t=91.9(t ) (c)t=92.0(t ) (d)t=92.1(t ) (e)t=92.2(t ) (f)t=92.3(t )0.05 0.00 0.05x (L )0.500.550.60 y ( L ) V / C s + V A -99 0 99 V )0.00.51.01.52.02.5 V / C s + V A a l o n g y - a x i s (g) )0.00.20.40.60.8 V e l o c i t y ( V ) (h) V at 91.8t C s V A V at 92.1t C s V A Figure 11.
The TS is destroyed by the well developed plasmoid, and then forms again above the plasmoid.Panels (a) ∼ (f) are Mach number and velocity divergence contours, and (g) shows Mach number distributionalong the y − axis; Panel (h) shows absolute velocity, local sound speed and local Alfv´en speed along the y − axis at two particular times using cyan lines(91.8 t ) and red lines(92.1 t ). above the closed loops turns to sub-magnetosonic. In this collision process, the TS has been totallydestroyed. At time 92.1 t , a weak shock front re-appears at higher altitude above this magneticisland. At time 92.3 t , we can see that a new TS has been restored as the plasmoid has mergedfully into the underlying flare loops and the super-magnetosonic reconnection outflows build up inthe upstream region again. The maximum Mach number then rises again to around 1.8. In panel ASTEX Dynamical TS y = 0 .
53, which is even lower than the shockfront before the collision. Therefore, the location of the TS front in y − direction can remarkable varyduring a short period.In this case, the variation of Mach number along the y − axis is shown in Figure 11(g). As the greendashed line shows, at time 92.0 t there is no TS front because the plasma is sub-magnetosonic. Wecan examine several characteristic speeds to explore why the shock front may disappear. In panel(h), we plot the plasma flow speed ( V ), sound speed( C s ), and Alfv´en speed ( V A ) along y − axis at twoparticular times, 91.8 t and 92.1 t . At the early time, 91.8 t , the shock front is at y = 0 .
538 wherethe flow speed drops from 0.37 to around 0.1. Until time 92.1 t , the magnetic island moves down tothe height y = 0 . V , lower than both the local Alfv´en speed and sound speed. Asshown by the red lines, the plasma speed is still significantly smaller than sound speed in surroundingregions. 3.4. Termination Shocks with Guide Fields
In this section, we compare the strength of the TS in Cases 3 through 5 with different strengthsof the guide field. The following simulations are performed by adding an initially uniform guidefield B z . The B z component might change strength during the evolution. Therefore, B z could affectthe local pressure equilibrium inside the current sheet significantly. The closed magnetic loops andcorresponding local characteristic speed could, in turn, vary with the changing guide field. In here,we consider three different initial values of the guide field B z = 0 . , .
25, and 0 . B (Cases 3, 4,5 in Table 1) in the simulation zone. The corresponding background plasma β is then equal to p ini / (0 . B ini + 0 . B z ), which are ∼ p ini = 0 .
05 and non-dimensional magnetic strength B ini = 1. We then analyze the TSproperties on the y − axis for above three cases.Figure 12(a) shows the maximum v/ (cid:112) C s + V A on the upstream side along y − axis at differenttimes. It is clear that the TS becomes weaker as the guide field increases to 0.5 (triangle points),and the corresponding maximum v/ (cid:112) C s + V A is less than ∼ . Shen et al. the strong guide field case ( B z = 0 . B z = 0 . , . ∼ . B z = 0 .
100 120 1401.01.52.02.53.0 M a x i m u m V / C s + V A (a)
100 120 1402.04.06.08.010.0 P r e ss u r e R a t i o (b)
100 120 140Time (t )1.01.52.02.53.03.5 D e n s i t y R a t i o (c)
100 120 140Time (t )1.01.52.02.53.03.54.0 T e m p e r a t u r e R a t i o (d) B z =0.05B z =0.25B z =0.50 Figure 12. (a) The maximum v/ (cid:113) C s + V A on the upstream side along the y − axis, and (b-d) compressionratios between upstream and downstream sides for these times with v/ (cid:113) C s + V A > y − axis. Theinitial guide fields are B z = 0.05, 0.25, and 0.5 B in the three cases (Cases 3-5 in Table 1). Effects of Thermal Conduction
ASTEX Dynamical TS ∼ . . ∼ . Shen et al. introduce a density variation along the TSs, which may be compared with future observational resultsof TSs made with radio spectral imaging .We perform statistical analysis during a short period and compare the shock compression ratiodistribution between two cases in Figure 15. We include a set of sampling points along the TS frontat each time within ∼ > ∼
35% in Case 7.Temperature ratios larger than 2 are slightly more frequent in thermal conduction case. IMPLICATIONSSimulations performed in this work revealed a dynamical nature of the flare TS driven by inter-mittent fast outflows of magnetic reconnection. Here we discuss briefly the implication of the aboveresults on the acceleration of particles at the TS and observations of TS signatures.As the upstream plasmoids and magnetic fluctuations interact with the magnetic loop in the loop-top region, the flare TS is likely to be turbulent and rippled, which may enhance the injection ratefor electron acceleration (Guo et al. 2012; Guo & Giacalone 2015). The typical compression ratioof 1.5–2.5 produced in the simulations is consistent with the value reported by Chen et al. (2015), C = n down /n up ≈ . f ( E ) ∝ E − δ with δ ∼ ASTEX Dynamical TS Figure 13.
Temperature and Mach number distribution in two simulations with (Case 6) and without(Case 7) thermal conduction. Red lines indicate TS fronts. for a large number of flare events at or near the limb, in which the intense footpoint HXR sourcesare occulted behind the limb and the weaker coronal HXR sources are revealed (Effenberger et al.2017). The strong variations in compression ratio in space and time should cause rapid variations inthe flux and spectrum of accelerated particles. Future particle acceleration studies need to take intoaccount the more realistic shock features.Magnetic configurations in post-shocked regions also show highly complex features in our numeri-cal experiments. The corresponding plasma compression and curved magnetic field lines can be seenin these regions as well due to the interaction between downflows and the loop top (as shown inFigure 10). These features suggest that energetic electrons have more opportunity to be trapped in2
Shen et al.
Figure 14.
Density and gas pressure distributions along the shock front (see red lines on Fig13(c)(d)) onthe post-shocked sides. Panel (a) is relative to Figure 13(c) and panel (b) is for Figure 13(d), respectively.The horizontal axis is the distance measured from the left-end of the shock front.
Figure 15.
Compression ratio distribution in two cases with and without thermal conduction.
ASTEX Dynamical TS ∼
10 MK plasmaand is routinely observed by the Interface Region Imaging Spectrograph (IRIS). They inserted anartificial termination shock into a Petschek-type reconnection geometry with plasmoid instability,and found that the termination shock leads to enhanced line intensity as well as blue- and/or red-shifted features (depending on the viewing geometry) in the Fe XXI line profile, which are mainlyassociated with the heated and compressed plasma in the downstream of the TS. They suggestedsuch signatures may be identified in the observed IRIS Fe XXI spectra. As our numerical simulationsshow, the dynamical TS causes the highly complex TS front and post-shocked plasma structures. Inthis case, the Fe XXI line profile may include more abundant fine features in either the line width oremission strength, which is worthwhile for further investigations in the future. SUMMARYIn this work, we performed a set of numerical experiments within the framework of the classicKopp-Pneuman flare configuration to investigate the formation and dynamics of TSs in solar flares.We focus on the region of the lower portion of the reconnection current sheet and the top of theclosed flare arcades where the TS is formed. We find that as the reconnection rate increases, thespeed of the reconnection outflow steadily grows and, at some point, exceeds the local magnetosonicspeed at the top of the flare arcades, thereby forming a TS at the looptop. We further show thatthe TS is highly dynamic, which responses promptly to the interaction between the intermittentreconnection flows in the shock upstream region and the obstacle in the downstream—the closed4
Shen et al. flare arcades. The TS front can be clearly delineated by the minimum of the velocity divergence,which allows us to analyze the physical parameters along the TS front in detail. To the contrary ofprevious assumptions for a steady-state, flat TS, however, our simulations show that the TS front canbe highly variable and asymmetric, displaying a variety of morphologies from flat, sloped, curved,to fragmentary or completely disrupted during the arrival of plasmoids of different sizes. We alsomake detailed measurements along the TS front on its compression ratio, fast mode Mach number,and inclination as it evolves. We find that, although they vary significantly with time and location,the fast mode Mach number and the density compression ratio ranges in M F ≈ C ≈ < . L to 0 . L ( < . ∼ ASTEX Dynamical TS ∼ ∼
3, and the pressure ratio hasnormally larger values ranging from 1 to 5, and even close to 10.4. Downward moving plasmoids can significantly reduce the strength of the termination shock.A well developed plasmoid may totally destroy the previous TS front, but a new one appears oncethe plasmoid merges into closed flare magnetic loops. This process leads to violent oscillation of TSheight. That is an important feature of intermittent magnetic reconnection.5. A background magnetic guide field causes both lower Mach numbers and compression ratios forTSs. A strong guide field obviously increases the local Alfv´en speed, which results in lower Machnumber. However, the impact of a weak guide field is relatively small. In cases with 5% or 25% guidefield, both the Mach number and compression ratio can still be larger than 2.The basic features of TSs in simulations with anisotropic thermal conduction are consistent withother cases that do not include it. Because the unsteady magnetic reconnection rate is not stronglyaffected by thermal conduction, it is not surprising that the maximum Mach number is around M ∼ Shen et al. . REFERENCES Aurass, H., Vr??nak, B., Mann, G., 2002, A&A,384, 273AAurass, H., Mann, G., 2004, ApJ, 615, 526ABhattacharjee, A., Huang, Y.-M., Yang, H., &Rogers, B. 2009, Physics of Plasmas, 16, 112102Blandford, R. and Eichler, D., 1987, PhysicsReports, 154, 1Beresnyak, A. 2017, ApJ, 834, 47 https://princetonuniversity.github.io/Athena-Cversion/ Carmichael, H. 1964, in The Physics of SolarFlares, ed. W. N. Hess (NASA SpecialPublication, Vol. 50; Washington, DC: NASA),451Chen, B., Bastian, T. S., Shen, C., et al. 2015, Sci,350, 1238Daughton, W., Roytershteyn, V., Karimabadi, H.,et al. 2011, Nature Physics, 7, 539Effenberger, F., Rubio da Costa, F., Oka, M., etal. 2017, ApJ, 835, 124Forbes, T. G. 1986, ApJ, 305, 553Forbes, T. G. 1988, SoPh, 117, 97
ASTEX Dynamical TS37