# Magnetohydrodynamic shocks in and above post-flare loops: two-dimensional simulation and a simplified model

## Abstract

Solar flares are an explosive phenomenon, where super-sonic flows and shocks are expected in and above the post-flare loops. To understand the dynamics of post-flare loops, a two-dimensional magnetohydrodynamic (2D MHD) simulation of a solar flare has been carried out. We found new shock structures in and above the post-flare loops, which were not resolved in the previous work by Yokoyama & Shibata (2001). To study the dynamics of flows along the reconnected magnetic field, kinematics and energetics of the plasma are investigated along selected field lines. It is found that shocks are crucial to determine the thermal and flow structures in the post-flare loops. On the basis of the 2D MHD simulation, we have developed a new post-flare loop model which we call the pseudo-2D MHD model. The model is based on the 1D MHD equations, where all the variables depend on one space dimension and all the three components of the magnetic and velocity fields are considered. Our pseudo-2D model includes many features of the multi-dimensional MHD processes related to magnetic reconnection (particularly MHD shocks), which the previous 1D hydrodynamic models are not able to include. We compare the shock formation and energetics of a specific field line in the 2D calculation with those in our pseudo-2D MHD model, and we found that they give similar results. This model will allow us to study the evolution of the post-flare loops in a wide parameter space without expensive computational cost and without neglecting important physics associated with magnetic reconnection.

## 1 Introduction

Solar flares are a phenomenon in which a huge amount of magnetic energy stored in the coronal field is rapidly released by magnetic reconnection. Where magnetic reconnection is a physical process in which a magnetic field in a highly conducting plasma changes its connectivity due to finite resistivity. The rapid energy release is accompanied by various plasma dynamics, like shock waves and the chromospheric evaporation. For more details about solar flares, the reader is referred to Shibata & Magara (2011).

Flares similar to solar flares have been observed from many stars and other astrophysical objects (e.g. Koyama et al., 1994; Tsuboi et al., 1998; Güdel, 2004). Since many of those flares and solar flares have many common features, solar flare models have been applied to those astrophysical flares (e.g. Hayashi et al., 1996; Shibata & Yokoyama, 2002; Machida & Matsumoto, 2003; Masada et al., 2010). Note that there should be differences in the plasma parameters between those flares and solar flares, like the plasma density, Alfven speed, and system size. Since these values determine the time scales of different processes, this can make the evolution of different flares distinct. Therefore, to comprehensively understand the flare physics in the universe, we need to explore various flares in a wide range of the plasma parameters.

There are many attempts to model solar and stellar flare loops. One-dimensional (1D) hydrodynamic models have been widely developed to study the thermal evolution and flows in a single loop. Depending on the assumed main energy transfer process, the models can be categorized into two branches: conduction-heating models (e.g. Nagai, 1980; Peres & Reale, 1993) and electron-beam-heating models (e.g. Fisher et al., 1985; Mariska et al., 1989). In general, the results of both models agree with observed emissions well. Hori et al. (1997, 1998) developed a pseudo-2D loop model based on 1D hydrodynamic calculations, and gave a simple description of observed soft X-ray emissions (see Warren (2006) for further development of this kind of models).

The first 2D magnetohydrodynamic (MHD) simulation of a solar flare based on a reconnection model with the heat conduction was performed by Yokoyama & Shibata (2001), and they studied the thermal evolution in the post-flare loops and the physics that determines the flare temperature. Shiota et al. (2005) performed MHD simulations of reconnection with the heat conduction to study the coronal mass ejections and associated giant arcade, where they do not include the chromosphere. Miyagoshi & Yokoyama (2004) carried out a MHD simulation of chromospheric evaporation jets formed as a result of reconnection between emerging magnetic flux and a coronal ambient field. Recently Longcope et al. (2009) and Longcope (2014) developed a theoretical model in which the dynamics of a reconnected field line is considered.

Shocks in and above the post-flare loops can be important for both the non-thermal particle acceleration and evolution of the thermal structure. The above-the-looptop hard X-ray source found by Masuda et al. (1995) has posed problems about the electron acceleration mechanism (for recent progress, see Krucker et al. (e.g. 2010) and Oka et al. (2015)). One of the promising scenarios is the acceleration by fast shocks which are expected to be formed above the soft X-ray post-flare loops (Somov & Kosugi, 1997; Tsuneta & Naito, 1998; Tanuma & Shibata, 2005). Recently acceleration from contracting plasmoids (Drake et al., 2006) and acceleration in plasmoids interacting with fast shocks (Nishizuka & Shibata, 2013) have been also proposed. It has been argued that a high-density region can be formed as a result of the shock-shock interaction at the top (Hori et al., 1997) and by the compression at termination shocks (Yokoyama & Shibata, 2001). Zenitani & Miyoshi (2011) carefully analyze shocks in a plasmoid in a ideal 2D MHD simulation. Despite the importance, the shock formation in and above the post-flare loops has not been yet investigated in detail.

How post-flare loops evolve in many astrophysical systems is one of our central interests. Since reconnection produces hot and high-speed plasma flows, it is expected that the thermal structure should be determined by a coupling among the plasma flows, shocks, and heat conduction. The heat conduction is essential to drive the hot high-speed evaporation flows. However, the calculation of the plasma dynamics with the heat conduction effects is computationally expensive. This makes the extensive parameter survey in the multi-dimensional MHD simulations difficult. The development of a simplified post-flare loop model based on a reconnection model is therefore desired.

In this paper, we study the shock formation and evolution of the thermal structure in and above the post-flare loops using MHD simulations. To understand the shock formation in 2D systems, a 2D MHD simulation of a solar flare has been carried out. We found new shock structures in and above the post-flare loops, which were not well resolved in the previous work by Yokoyama & Shibata (2001). To study the dynamics of flows along the reconnected magnetic field in detail, kinematics and energetics of the plasma are investigated along selected field lines. It is found that shocks in the post-flare loops are crucial to determine the thermal and flow structures in the post-flare loops. On the basis of the results of the 2D MHD simulation, we have developed a new post-flare loop model which we call the pseudo-2D MHD model. The model is based on the 1D MHD equations, and has been developed to model a 3D reconnected field line. Through a comparison, it is found that the shock formation and thermal evolution in the pseudo-2D MHD model are similar to those in the 2D MHD model.

This paper is structured as follows. In Section 2, the 2D MHD model of a solar flare is introduced. In Section 3, numerical results of the 2D MHD simulation are detailed. In Section 4, on the basis of the 2D MHD simulation, we develop our pseudo-2D MHD model of a post-flare loop. In Section 5, numerical results of the pseudo-2D MHD model are introduced and are compared with the 2D MHD simulation. Section 6 contains our summary and discussion.

## 2 2D MHD Model of a Solar Flare

### 2.1 Assumptions and Basic Equations

We performed a 2D MHD simulation of a solar flare similar to Yokoyama & Shibata (2001) simulations. We take a rectangular calculation domain in the - plane. The MHD equations in the following form are solved:

(1) | ||||

(2) | ||||

(3) | ||||

(4) | ||||

(5) | ||||

(6) | ||||

(7) | ||||

(8) | ||||

(9) |

where and . ( in cgs units) is the coefficient of the Spitzer-type conductivity, is the specific heat ratio (5/3 is used in this study), is the mean particle mass and is the Boltzmann constant. is the electric resistivity. is the conduction flux. The normalization units of our simulations are summarized in Table 1. The radiative cooling time is expected to be longer than the dynamical time in the post-flare loops. Since we mainly focus on the dynamical processes such as the shock formation and flows, we neglect the radiative cooling.

### 2.2 Initial and Boundary Conditions

Our model is identical to the Yokoyama & Shibata (2001) model, except for the domain size and boundary conditions. The initial condition is same as their model. The initial and boundary conditions are as follows. The domain is and , where =10 and . A schematic diagram of the initial and boundary conditions are shown in Figure 1. The initial magnetic field is assumed to be a force-free field and is given by

(10) | ||||

(11) | ||||

(12) |

where is the width of the initial current sheet. The gas pressure is assumed to be uniform (). The density distribution is described as

(13) |

where and are respectively the densities in the corona and chromosphere, and and are respectively the height and the width of the transition region between the corona and chromosphere. and are set to 1 and 0.2, respectively. The chromosphere is modeled as a dense and cool plasma. For simplicity, in this paper is set to . is the initial coronal temperature.

To allow the magnetic field to reconnect, we impose a localized resistivity in the form of

(14) |

where and . We localize the resistivity and fix it with time to realize a fast and quasi-steady magnetic reconnection with a single X-point (Ugai, 1992).

All the boundaries are symmetric. At the boundary at the sign of is reversed.

The numerical scheme is based on a Harten-Lax-van Leer (HLL) scheme developed by Miyoshi & Kusano (2005), HLLD (”D” stands for Discontinuities), which is a shock-capturing scheme. It has the second-order accuracy in space and time. The heat conduction term is solved with an implicit method similar to Yokoyama & Shibata (2001) method to reduce the calculation time. We modify their method to more accurately calculate the heat conduction flux. A detailed description of the implicit method is given in Appendix A. The calculated domain is resolved with grids.

## 3 Numerical Results of 2D MHD Model

### 3.1 Overview of Evolution

Figure 2 displays an overview of the time evolution of the 2D MHD simulation. The region where is also shown only for visual inspection. Reconnection takes place due to the localized resistivity and the Alfvenic collimated reconnection outflow is produced. The reconnected field lines pile up in the outflow region to form the growing loop system at the base of the corona. The temperature structure is smooth along the magnetic field owing to the heat conduction. The heat in the hot outflow is carried to the chromosphere along the magnetic field due to the heat conduction. As a result, the upper chromosphere is heated up and rapidly expands, leading to the formation of the hot dense upflows called the chromospheric evaporation. The evaporated plasma finally fills the reconnected field lines and forms the hot dense post-flare loops. We note that the weak disturbance which starts to propagate at the beginning of the simulation is generated because the initial condition is not in the thermal equilibrium between the chromosphere and corona. We confirmed that its effect is negligible for the dynamics of the post-flare loops.

One-dimensional plots parallel to the axis across (i.e. across the reconnection outflow) are shown in Figure 3. It is shown that a Petschek-type reconnection is established because of the localized resistivity, where two slow shocks emanate from the reconnection region. The slow shocks can be discerned as a pair of the discontinuities in Figure 3. We note that due to the heat conduction the slow shocks become isothermal slow shocks (Forbes et al., 1989; Yokoyama & Shibata, 1997; Chen et al., 1999; Seaton & Forbes, 2009).

The thermal structure of the post-flare loops is determined by a complex coupling among the plasma flows, shocks, and heat conduction. When the reconnection outflow collides with the loop system below, two oblique fast-mode shocks and sometimes a horizontal fast-mode shock (Mach disk) are formed at the top (see Figure 4), which indicates that the termination shock is a combination of two or three fast-mode shocks. Note that most of the regions where takes large negative values are slow or fast shocks. We found that the strength of the termination shock shows a quasi-periodic oscillation. We also found the shock reflection and Mach reflection in a concave magnetic field at the top (Figures 4 and 6). The relationship between the oscillation and the flow structure at the top will be discussed in our future papers.

The high-pressure plasma at the loop-top expands along the magnetic field to the foot-points, generating the high-speed downflows. Then the evaporation flows collide with the downflows, forming the dense regions in the post-flare loops (”humps” in Figure 5, named by Yokoyama & Shibata (2001)). After the collision, the fronts of the evaporation flows and downflows becomes slow shocks (see maps). The pair of the upward slow shocks finally collide with each other at the apex, forming another dense region (see Figure 6). Note that the ”blob” named by Yokoyama & Shibata (2001) is different from this high-density region. The blob is the high-density region in a concave magnetic field at the top.

Figure 7 displays the slow shocks in the post-flare loops mentioned above. The entropy is discontinuous at the two pairs of the discontinuities in the density map, but the temperature is smooth along the magnetic field, indicating that they are isothermal slow shocks. The plasma in a large fraction of the post-flare loops, as well as in the outflow region, is larger than unity (the contour indicates the level where ), meaning that the low- approximation cannot be applied to the reconnected magnetic field. This is a consequence of the shock heating and compression. The high-pressure post-flare loops are confined by the surrounding low- plasma. This is consistent with the observation by McKenzie (2013); they concluded that plasma is of the order of unity or larger than unity in the supra-arcade plasma in two flares analyzed.

### 3.2 Dynamics and Energetics along a Specific Field Line

We have seen the two-dimensional evolution of a solar flare. To study the dynamics of flows along the reconnected magnetic field in detail, kinematics and energetics of the plasma are investigated along selected field lines. We performed the same analysis for other field lines and confirmed that they give similar results.

We picked up a magnetic field line and measured the physical quantities along it. The field line used in the analysis is shown in Figure 8. Figure 9 displays the time-distance diagrams of the density, temperature, and along the field line whose foot-point is located at , where is the local sound speed. The distance is measured along the field line. Before the field line reconnects, the starting point of the distance is the intersection point between the field line and the top boundary. After reconnection, the starting point is the apex of the closed loop. The sign of is defined as positive and negative when a plasma flow travels to the apex and foot-points, respectively.

In Figure 9, we can discern the field line shrinking after reconnection. The heat released at the Petschek slow shocks is transferred to the chromosphere (see the temperature in Figure 9), leading to the generation of the chromospheric evaporation. The evaporation flows are seen as upflows from the chromosphere (see and maps). The reconnection outflow is significantly decelerated at . This is seen as the enhancement of the density and temperature and negative (compression), indicating that the kinetic energy of the outflow is converted to the thermal energy. After the termination, the high-pressure plasma at the top expands along the magnetic field, forming the downflows. The downflows collide with the chromospheric evaporation flows, leading to the formation of the humps (see also Figure 5). After the collision, the fronts of the evaporation flows and downflows become steep and finally become slow shocks. The pair of the upcoming shocks finally cross each other at the apex at , forming another high-density region (see also Figure 6).

We found slow shocks propagating along the field line, but there is no prominent shocks nor waves propagating back and forth from end to end. They are damped by the heat conduction (e.g. Ofman & Wang, 2002). Also note that the propagation speed of the shocks is largely decelerated by the evaporation flow (i.e. Doppler effect), which significantly affects the traveling time of the slow-mode waves/shocks (the local sound speed is , but the propagation speed is ). That is, a simple estimation by is not a good approximation of the traveling time, where and is the loop length and the sound speed in the loop, respectively. We observe no prominent signature of the standing slow-mode waves in the calculated time range (as for acoustic waves in the post-flare loops, see e.g. Nakariakov et al. (2004)).

Figure 10 shows the time evolution of the total, magnetic, internal (thermal) and kinetic energies integrated along two specific field lines. The top and bottom panels are for the field lines which are rooted at and , respectively. Considering that the cross-sectional area of the flux tube is inversely proportional to the magnetic field strength, we integrate the energies along a field line as follows:

(15) | ||||

(16) | ||||

(17) | ||||

(18) |

In this paper, all the energies are normalized by the initial total energy. After the field lines reconnect (), the magnetic energy is rapidly converted to the internal and kinetic energies. As for the field line rooted at , the maximum value of the kinetic energy is 0.35. Let us define the time when becomes as . If we compare the total, magnetic and internal energies at and those at , , , and . Note that the variation of the total energy remains within %. The field line rooted at also gives a similar result. The small variation in the total energy reflects the fact that the compressive and expanding motions by the surrounding plasma are localized in time and space (see Figure 9). Therefore, the compressive motion of the reconnected flux tube is found to be not significant with respect to the total energy variation.

## 4 Pseudo-2D MHD Model of a Post-Flare Loop

### 4.1 Physical Processes Considered

Through the analysis of the 2D MHD simulation, we found that the field-aligned plasma motions (particularly evaporation flows and slow shocks) and heat conduction seem to mainly determine the dynamics in the post-flare loops. Fast shocks are important for converting the kinetic energy of the reconnection outflow to the heat and for locally changing the cross-sectional area, but they do not largely change the total energy of the field lines (Figure 10). On the basis of the results, we aim to model the multi-dimensional reconnection and flare processes in a simplified MHD scheme.

What is important for conduction-heating-type flare modeling would be 1. to model a reconnected field line, 2. to model the reconnection inflow and outflow, 3. to include MHD waves (since MHD waves can carry a large amount of the released energy from the reconnection sites (Kigure et al., 2010)), 4. termination of the reconnection outflow and energy conversion of the kinetic energy of the outflow into the thermal energy, and 5. to include the heat conduction which is essential to determine the temperature of the flare loops.

Considering that the plasma motions are frozen-in a magnetic field, a 1D MHD model will be the simplest form among the possible MHD models. We developed a model based on the 1D MHD equations, where all the variables depend on one space dimension and all the three components of the magnetic and velocity fields are considered. We regard our model as a pseudo-2D MHD model. Note that the meaning of ”pseudo-2D” of our model is different from that of Hori et al. (1997) hydrodynamic model which consists of isolated and fixed semi-circular loops with different lengths and constant cross-section.

### 4.2 Assumptions and Basic Equations

We take a cartesian coordinate system in which all the variables are functions of . The reconnection outflow is in the -direction. is the ignorable coordinate (uniform in the -direction).

We mimic a reconnected magnetic field line by assuming a magnetic field with a sharp bend (see Figure 11). The magnetic field line shrinks with time and drives the flow perpendicular to the field which represents the reconnection outflow.

A schematic picture of our model is as follows. Figure 12 displays a comparison of our pseudo-2D MHD model with the 2D MHD model. is the height of the reconnection point, and is the location of the foot-point of the field line. They are linked by the relation of . The reconnection angle and the initial plasma beta are treated as free parameters. If a guide-field (the -component of the magnetic field) is included, the guide-field angle will be an additional parameter.

The shrinking motion of the reconnected field line will stop when it collides with the magnetic loops piling up below. To model this process, the final configuration of the field line is given in the model and the outflow is decelerated when the field line approaches the final state. The termination process is modeled by adding a damping term to the equation of motion perpendicular to the moving field line. We let the damping term work only when the reconnected field line comes close to the final configuration.

According to Figure 10, the total energy of a field line is conserved within several 10 %. On the basis of this result, we hypothesize that the total energy in a reconnected flux tube is conserved. The cross-sectional area in our model is assumed to be uniform and constant in time.

The basic equations are as follows:

(19) | ||||

(20) | ||||

(21) | ||||

(22) | ||||

(23) | ||||

(24) | ||||

(25) | ||||

(26) | ||||

(27) |

where , , , , and . All the physical quantities are only dependent of and . Note that we include a damping term in the momentum equations that slows down only the x and y-components of the velocity perpendicular to the magnetic field. The detailed description of the damping term is given below. Note that the total energy is conserved along a field line, and that the kinetic energy reduced by the damping term is converted into the thermal energy. The normalization units are the same as those in the 2D MHD model (See Table 1). Note that this model can treat slow-mode, fast-mode and Alfven waves.

Considering the symmetry, we only solve the domain within . At the boundary of the mirror symmetric boundary conditions are imposed. The boundary of is free.

The initial conditions are as follows. The free parameters that determine the initial magnetic field are the plasma beta , and the angles and (see Figure 12). The guide field effect is detailed in Appendix B. The initial magnetic field is given by

(28) | ||||

(29) | ||||

(30) | ||||

(31) |

The gas pressure is assumed to be uniform (). The domain is divided into the two layers, namely chromosphere and corona:

(32) |

where is the location so that . Where is a function that describes the configuration of the field line.

The initial and final magnetic field configurations are described as follows. See also Figure 13. The reconnection point (the location of the sharp bend) is assumed to be at . The initial condition of a magnetic field is written as

(33) |

for . The final state of the magnetic field is approximated by a quadratic function of

(34) |

which has the same slope of the tangent line with the equation (33) at .

We virtually consider the travel distance of the reconnected field line in the -direction. When the field line approaches the final state, the damping term is applied only to the - and -components of the velocity perpendicular to the field line. We define the distance in the -direction between the field line at and the field line in the final state as

(35) |

The damping term only works when the field line approaches to the final state:

(36) |

where , is the outflow speed in the y-direction, and is a free parameter that denotes a typical braking distance. To prevent the field line from shrinking further even after the travel time , we increase the damping coefficient by a factor of 100 after the time . Note that the term arising from the damping terms is not included in the energy equation under the assumption that the total energy along a field line is conserved. The kinetic energy decreased by the damping term is converted only into the thermal energy.

The reconnection angle and plasma beta are important parameters to determine the total released magnetic energy and energy conversion rate (say, reconnection rate). To choose a physically acceptable parameter set, we utilize the analytical approach by Falle et al. (1998). A detailed description is presented in Appendix B.

The height of the reconnection point is assumed to be , and the reconnection angle is . The domain size is therefore . The width of the transition region and the typical damping distance are respectively and .

The numerical scheme of the pseudo-2D MHD model is based on the Vögler et al. (2005): the fourth-order space-centered difference for spatial derivative and an explicit four-step Runge-Kutta time integration. We explicitly solve the heat conduction term. Using the current computational resources, it is not difficult to explicitly solve the heat conduction in our 1D calculations. The domain is resolved by 640 grids.

## 5 Numerical Results of Pseudo-2D MHD Model

### 5.1 Dynamics and Energetics

Figure 14 demonstrates the time evolution of the magnetic field structure. The field line retracts and sweeps up the plasma like a slingshot. The reconnection outflow is decelerated when the field line approaches the assumed final configuration. In the following, we performed the same analysis as for the 2D MHD model.

One-dimensional plots shown in Figure 15 demonstrate the formation of a pair of two slow shocks attached to the reconnection outflow, that is, the Petschek-type slow shocks. As well as the 2D MHD model (Figure 3), the Petschek-type shocks are isothermal shocks due to the heat conduction.

Figure 16 compares the field-aligned motions in the pseudo-2D MHD model and 2D MHD model. As well as in the 2D MHD model, the heat generated at the Petschek slow shocks is transferred to the chromosphere by the heat conduction, forming the chromospheric evaporation. The humps and high-density regions at the top are found to be formed in the same way as in the 2D MHD model: the downflow-evaporation collision and shock-shock interaction, respectively. The damping of the slow shocks are also observed.

The energy evolution is examined in the pseudo-2D MHD model. When we integrate the energies, we drop off the term in the integrand because in the pseudo-2D MHD model the variation of the cross-sectional area is not considered. Figure 17 is the same as Figure 10 but for the pseudo-2D MHD model. It is shown that , , and , which is similar to the results of the 2D MHD model.

### 5.2 Dependence on Parameters

Previously a formula that determines the post-flare loop temperature was derived under the assumption that the energy input to a loop balances with the conduction cooling rate (Fisher & Hawley, 1990). The formula is given by

(37) |

where is the volumetric heating rate and is the half length of the magnetic field line of a post-flare loop. The heating rate by magnetic reconnection is determined by the Poynting flux: . Using this, the temperature can be written as

(38) |

by assuming that . This scaling law was derived by Yokoyama & Shibata (1998, 2001). We checked whether the scaling law based on a magnetic reconnection model holds in our pseudo-2D MHD model.

Figure 18 shows the numerically-obtained -, - and - relations. denotes the temperature in the reconnection outflow. denotes the maximum temperature after the pair of the two slow shocks generated by the chromospheric evaporation flows collides at the apex. As shown in Figure 18, it is found that the temperature in the loop obeys the scaling law well.

## 6 Summary and Discussion

In this paper we investigated the flow structure, shock formation and thermal evolution in the post-flare loops using MHD simulations. On the basis of the 2D simulation result, we have developed a new post-flare loop model (the pseudo-2D MHD model). We compare the flow structure, shock formation and energetics of a specific field line in the 2D calculation with those in the pseudo-2D MHD model, and then we find that they give similar results. Here we summarize the results and compare our results with previous studies.

Performing a 2D MHD simulation, we found new shock structures. The termination shock consists of two oblique fast-mode shocks and sometimes a horizontal fast-mode shock (Figure 4). A hump is formed as a result of the collision of the downflow and the chromospheric evaporation flow. After the collision, the fronts of the evaporation flow and downflow become slow shocks, and then the hump appears as a dense region behind the two slow shocks (Figure 5). The upward slow shock finally interact with the slow shock coming from the other side at the top, forming the high-density region. Note that the high-density region is separated from the blob in Yokoyama & Shibata (2001) and is formed below it (Figure 6).

We found that the strength of the termination shock in the 2D MHD model shows a quasi-periodic oscillation, which is a multi-dimensional feature. In addition, the shock reflection and Mach reflection are sometimes observed in a concave magnetic structure, which could be important for understanding the heating in the loop-top. These complicated structures at the top will be detailed in our future papers.

We observe no prominent shocks nor waves propagating back and forth from end to end. In addition, no prominent standing slow-mode waves are observed within the calculated time range. By performing 2D MHD simulations without the heat conduction, we confirmed that without the heat conduction the slow shocks formed at the fronts of the downflows from the top propagate back and forth from end to end. Thus the propagation of the slow shocks in the post-flare loops is found to be significantly affected by the heat conduction and evaporation flows. The slow shocks are damped by the heat conduction. The propagation speed of the slow shocks is reduced by the evaporation flow (Doppler effect), which makes the shock propagation time longer. Therefore it is essential to consider the flows resulting from reconnection (particularly downflows from the top) for understanding the behavior of magneto-acoustic waves in the post-flare loops.

In fully 3D situations, reconnection can be intermittent in space and time, which could affect the shock structures found in this study. 3D component reconnection, where reconnecting magnetic field lines are not perfectly anti-parallel, can result in the reconnection outflow jet with a speed insufficient for the formation of the fast shock above the loop top. However, in the case that the reconnection outflow speed exceeds the fast mode phase speed in the outflow region, we expect the formation of the termination shock found as in our 2D MHD model. We also expect that the slow shocks presented in the 2D MHD model will be formed in 3D, since the field-aligned plasma flows which are essential to form slow shocks are well treated in our MHD flare modeling scheme.

From the 2D simulation, we found that the field-aligned plasma motions (particularly evaporation flows and slow shocks) and heat conduction mainly determine the dynamics in the post-flare loops. Considering this, we construct the pseudo-2D MHD model which is basically a 1D MHD model. The pseudo-2D MHD model is compared with the 2D MHD model, and we found that the dynamics (particularly flow structure and shock formation) and the energetics are similar between the two models. The scaling law for the temperature based on a reconnection model is also examined, and it is found that the scaling law holds in the pseudo-2D MHD model. These facts indicate that our pseudo-2D MHD model captures important features of a reconnection model of solar flares.

1D hydrodynamic models, like Mariska et al. (1989) and Hori et al. (1997, 1998), have been used for the post-flare loop modeling. The models are useful to study the thermal evolution and flows in the flare loops, but the energy input (in many cases the heat input) must be done by hand. Our model includes many features of the multi-dimensional MHD processes related to magnetic reconnection, like the heating by the Petschek slow shocks and conversion of the kinetic energy of the reconnection outflow to the heat. Another important point is that our model can treat MHD waves and shocks generated in the post-flare loops, which could be important to understand the flow structure in the loops.

We note that all the previous 1D hydrodynamic models lack the strong downflows from the top (for example, see Figure 5 in Hori et al. (1997)). We showed that the downflows play important roles in forming the humps and slow shocks (Figures 5 and 6). Our pseudo-2D MHD model naturally produces the downflows from the top, allowing us to study the dynamic evolution of the thermal structure.

A theoretical model in which field lines shorten after localized 3D reconnection based on the thin flux tube approximation was proposed by Longcope et al. (2009). The model assumes that the plasma has always low-. Our 2D MHD simulation demonstrates that, because of the shock heating and compression, the plasma becomes larger than unity not only in the reconnection outflow but also in a large part of the post-flare loops (see Figure 7). In such regions, we need to consider the effects of the gas pressure to understand the flow and shock structures. Our pseudo-2D MHD model can treat the high- plasma. Longcope et al. (2009) model could be valid in the situation in which the shock heating does not break the low- assumption (e.g. in the situation in which the guide-field is much stronger than the reconnection field). Regarding this issue, see also Appendix B.

The high- condition could lead to disordering of the magnetic field when turbulence is important so that the geometrical assumptions made for our pseudo-2D MHD model are not met. We found no evidence that turbulence is important in the current sheet in our 2D MHD model. Our model will be useful to model the reconnected field lines in such a situation in which turbulence is not important.

As a result of the approximations made for the pseudo-2D MHD model, some multi-dimensional processes such as termination shock formation and turbulence cannot be modeled. However, the pseudo-2D MHD model are able to treat plasma flows and waves/shocks (slow-mode waves/shocks and Alfven waves) along the magnetic field, which could be useful to understand the plasma motions in the post-flare loops.

Our pseudo-2D MHD model requires much smaller computational cost than other multi-dimensional MHD models. This model will allow us to study the evolution of the flare loops in a wide parameter space without expensive computational cost. Also, it will be much easier to include detailed physics like the non-equilibrium ionization effect (e.g. Imada et al., 2011). These will be our future work.

## Appendix A Numerical Method of Heat Conduction in 2D MHD Model

We modify the time step splitting method to calculate the heat conduction flux by Yokoyama & Shibata (2001). Using the MHD energy flux and heat conduction energy flux , we can write the energy equation as follows:

(A1) |

where is the total energy. The discretized form is

(A2) |

First, we calculate the MHD part

(A3) |

where the superscript denotes the results of the MHD step. Then we calculate the heat conduction step

(A4) |

The procedure mentioned above is the same as Yokoyama & Shibata (2001). The difference appears in the expression of the heat flux formula. The heat flux formulae in Yokoyama & Shibata (2001) are

(A5) | |||||

(A6) | |||||

(A7) | |||||

(A8) |

where the subscripts and denote and . The time and space discretization are, for example,

(A9) |

We modify this to

(A10) |

The difference of two formulae appears in the operator matrix of the implicit scheme used. The number of the non-zero components in one row is 5 in Equation (A9), and 9 in Equation (A10). Geometrically, in Equation (A10) we use all the neighboring 8 grids around each point at step to calculate the heat flux (Figure 19). This method gives more accurate results than the previous method, particularly in the grid points where a magnetic field is largely bend and oblique to coordinate.

## Appendix B Analytical Approach to the Riemann Problem: Possible Solutions

We shall describe the exact solutions of the symmetric MHD Riemann problems that can appear in our pseudo-2D MHD model. In a symmetric MHD Riemann problem, the initial states in the left and right hand side are assumed to be as follows.

(B1) | |||||

(B2) | |||||

(B3) | |||||

(B4) | |||||

(B5) | |||||

(B6) |

where the subscripts and represent the left hand side and the right hand side, respectively.

Let us consider the situation without the guide field component (). Since the magnetic field vectors in the initial states are in the - plane, the rotational discontinuities will not appear in the solutions. Due to the initial discontinuity, one contact discontinuity, two slow shocks, and two fast rarefaction waves are generated.

From the jump condition across the contact discontinuity, the tangential magnetic field, , should be equal to zero in the region between two slow shocks. Therefore, it is necessary to have either the switch-off slow shock (SS) or the switch-off fast rarefaction wave (FRW). In the former case (see the left panels of Figure 20), the exact solution consists of the switch-off SS and FRW. On the other hand, in the latter case (the middle panels of Figure 20), the exact solution consists of the pure hydro shock (HS) and switch-off FRW.

There are two parameters in the symmetric MHD Riemann problems: plasma beta and the reconnection angle . A sophisticated procedure to derive the exact solutions of more general MHD Riemann problems was developed by Falle et al. (1998). We utilize the procedure to derive the exact solutions of our problems. The result is summarized in the phase diagram of Figure 21, where and are the key parameters. Here the heat conduction is neglected. The - space is divided into two regions: switch-off slow shock regime and pure hydro shock regime. The parameter set of the example introduced in this paper is within the switch-off slow shock region. If we include the heat conduction effect, the pure hydro shock region will become slightly wider.

Figure 22 displays the normalized reconnection rate in the - diagram, where the normalized reconnection rate is defined by . Where is the outflow speed, , and . The maximum reconnection rate observed in solar flares and simulations is of the order of 0.1. To take the reconnection rate similar to 0.1 in our model, we need to choose a reasonable parameter set from this diagram.

If a guide field exists, the rotational discontinuities (RD) will appear in addition to SS and FRW. In this case it is necessary to have either the switch-off SS, switch-off FRW, or RD because of the boundary condition at the contact discontinuity (). Among them, only the solution with RD (the right panels of Figure 20) meets the requirement that the continuity condition of at the contact discontinuity should be satisfied (Petschek & Thorne, 1967). Longcope (2014) claimed that pure hydro shocks can be formed in a reconnected flux tube even when a guide field exists. However, considering the analysis here, pure hydro shocks will not appear when a guide field exists.

Quantity | Unit | Value |
---|---|---|

Length | 3,000 km | |

Velocity | 170 km s | |

Time | 18 s | |

Temperature | K | |

Density | (10 cm ) | |

Pressure | 0.47 dyn cm | |

Magnetic field strength | 3.4 G |

### Footnotes

- affiliation: Kwasan and Hida Observatories, Kyoto University, Yamashina, Kyoto 607-8471, Japan
- affiliation: Institute of Space and Astronautical Science (ISAS), Japan Aerospace Exploration Agency (JAXA), 3-1-1 Yoshinodai, Sagamihara, Kanagawa 252-5210, Japan
- affiliation: Kwasan and Hida Observatories, Kyoto University, Yamashina, Kyoto 607-8471, Japan
- affiliation: Kwasan and Hida Observatories, Kyoto University, Yamashina, Kyoto 607-8471, Japan

### References

- Chen, P. F., Fang, C., Tang, Y. H., & Ding, M. D. 1999, ApJ, 513, 516
- Drake, J. F., Swisdak, M., Che, H., & Shay, M. A. 2006, Nature, 443, 553
- Falle, S. A. E. G., Komissarov, S. S., & Joarder, P. 1998, MNRAS, 297, 265
- Fisher, G. H., Canfield, R. C., & McClymont, A. N. 1985, ApJ, 289, 414
- Fisher, G. H., & Hawley, S. L. 1990, ApJ, 357, 243
- Forbes, T. G., Malherbe, J. M., & Priest, E. R. 1989, Sol. Phys., 120, 285
- Güdel, M. 2004, A&A Rev., 12, 71
- Hayashi, M. R., Shibata, K., & Matsumoto, R. 1996, ApJ, 468, L37
- Hori, K., Yokoyama, T., Kosugi, T., & Shibata, K. 1997, ApJ, 489, 426
- Hori, K., Yokoyama, T., Kosugi, T., & Shibata, K. 1998, ApJ, 500, 492
- Imada, S., Murakami, I., Watanabe, T., Hara, H., & Shimizu, T. 2011, ApJ, 742, 70
- Kigure, H., Takahashi, K., Shibata, K., Yokoyama, T., & Nozawa, S. 2010, PASJ, 62, 993
- Koyama, K., Maeda, Y., Ozaki, M., et al. 1994, PASJ, 46, L125
- Krucker, S., Hudson, H. S., Glesener, L., et al. 2010, ApJ, 714, 1108
- Longcope, D. W., Guidoni, S. E., & Linton, M. G. 2009, ApJ, 690, L18
- Longcope, D. W. 2014, ApJ, 795, 10
- Machida, M., & Matsumoto, R. 2003, ApJ, 585, 429
- Mariska, J. T., Emslie, A. G., & Li, P. 1989, ApJ, 341, 1067
- Masada, Y., Nagataki, S., Shibata, K., & Terasawa, T. 2010, PASJ, 62, 1093
- Masuda, S., Kosugi, T., Hara, H., et al. 1995, PASJ, 47, 677
- McKenzie, D. E. 2013, ApJ, 766, 39
- Miyoshi, T., & Kusano, K. 2005, Journal of Computational Physics, 208, 315
- Miyagoshi, T., & Yokoyama, T. 2004, ApJ, 614, 1042
- Nagai, F. 1980, Sol. Phys., 68, 351
- Nakariakov, V. M., Tsiklauri, D., Kelly, A., Arber, T. D., & Aschwanden, M. J. 2004, A&A, 414, L25
- Nishizuka, N., & Shibata, K. 2013, Physical Review Letters, 110, 051101
- Ofman, L., & Wang, T. 2002, ApJ, 580, L85
- Oka, M., Krucker, S., Hudson, H. S., & Saint-Hilaire, P. 2015, ApJ, 799, 129
- Peres, G., & Reale, F. 1993, A&A, 275, L13
- Petschek, H. E., & Thorne, R. M. 1967, ApJ, 147, 1157
- Seaton, D. B., & Forbes, T. G. 2009, ApJ, 701, 348
- Shibata, K., & Yokoyama, T. 2002, ApJ, 577, 422
- Shibata, K., & Magara, T. 2011, Living Reviews in Solar Physics, 8, 6
- Shiota, D., Isobe, H., Chen, P. F., et al. 2005, ApJ, 634, 663
- Somov, B. V., & Kosugi, T. 1997, ApJ, 485, 859
- Tanuma, S., & Shibata, K. 2005, ApJ, 628, L77
- Tsuboi, Y., Koyama, K., Murakami, H., et al. 1998, ApJ, 503, 894
- Tsuneta, S., & Naito, T. 1998, ApJ, 495, L67
- Ugai, M. 1992, Physics of Fluids B, 4, 2953
- Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335
- Warren, H. P. 2006, ApJ, 637, 522
- Yokoyama, T., & Shibata, K. 1997, ApJ, 474, L61
- Yokoyama, T., & Shibata, K. 1998, ApJ, 494, L113
- Yokoyama, T., & Shibata, K. 2001, ApJ, 549, 1160
- Zenitani, S., & Miyoshi, T. 2011, Physics of Plasmas, 18, 022105