# Simulations Reveal Fast Mode Shocks in Magnetic Reconnection Outflows

###### Abstract

Magnetic reconnection is commonly perceived to drive flow and particle acceleration in flares of solar, stellar, and astrophysical disk coronae but the relative roles of different acceleration mechanisms in a given reconnection environment are not well understood. While outflow fast mode shocks have been predicted analytically, we show for the first time via direct numerical simulations that such shocks do indeed occur in the outflows of fast reconnection when an obstacle is present. These shocks are are distinct from the slow mode Petschek inflow shocks. If Fermi acceleration of electrons operates in the weak fast shocks the associated compression ratios will induce a Fermi acceleration particle spectrum that is significantly steeper than strong fast shocks commonly studied, but consistent with the demands of solar flares. While this is not the only acceleration mechanism operating in a reconnection environment, it is plausibly a ubiquitous one.

PACS codes: 96.60.Iv, 52.35.Vd, 95.30.Qd, 52.35.Tc, 96.60.qe, 98.54.Cm

July 17, 2019

## I Introduction

Magnetic reconnection is a frequently studied process in plasma astrophysics whereby magnetic energy is converted into some combination of thermal, non-thermal, and flow kinetic energy. Reconnection underlies models of solar and stellar flares and coronae of astrophysical accretion disks. The process has also been studied in laboratory plasma experiments yamadakulsrud ().

Most work has focused on understanding the reconnection rate Birn et al. (2001) and under what circumstances fast reconnection occurs. Several distinct paradigms for the onset of fast reconnection in nature have been studied and may operate in complementary environments. The first occurs when the ion inertial length is larger than the magnetohydrodynamic (MHD) Sweet-Parker resistive length, facilitating current driven plasma instabilities that enhance dissipation near the X-point biskamp00 (); kulsrud01 (); uzdensky03 (). A second occurs in the context of turbulent MHD, where the global reconnection rate is enhanced by the contemporaneous action of many small scale sites loureiro07 (); kowal09 (); eyink11 (). A third path to fast reconnection occurs for a long enough laminar MHD current sheet that becomes tearing mode unstable and again multiple small scale reconnection sites acting together enhance the rate. bhy2009 ().

There has been less synthesized progress understanding how the converted magnetic energy partitions and on the accelerated particle spectra although different mechanisms have been investigatedromanovalovelace92 (); blackman1 (); blackman96 (); dalpino (); kowal (); blackman97 (); tsuneta97 (). Because of the variety of plasma reconnection conditions, it is best to think of reconnection as an acceleration environment rather than a single mechanism. Different processes can operate with different strengths depending on the conditions, and whether reconnection is fast or slow.

The basic question we address here, using numerical simulations, is whether fast shocks arise generically in reconnection outflows. Analytic predictions and phenomenological scenarios suggest that weakly compressive fast shocks may form downstream of reconnection sites blackman1 (); blackman97 (); tsuneta97 (); tsuneta98 (), such as in coronal flares if supersonic flows from the X-point impinge on a previously reconnected loop top from above. However, there has yet to be a study of this specific such shock formation numerically. Note that the fast outflow shocks we are focusing on are distinct from slow mode inflow shocks of the Petschek model petschek (). Even weak fast shocks in reconnection outflows would be important because of their potential ubiquity and role for shock-Fermi acceleration bell2 (); JandE (). A fast downstream shock in solar flares blackman1 (); tsuneta98 (); blackman2 () also benefits from magnetic dissipation near the X-point that may provide the needed injection electrons that seed the Fermi process.

Here we present 2-D magnetohydrodynamic (MHD) simulations of initial Harris current sheet configurations that are unstable to magnetic reconnection, and determine the conditions for which the outflows are weakly super fast-magnetosonic. By injecting a dense plasma obstacle into the downstream, we then directly show the formation of weak outflow shocks. From the numerically computed compression ratios, we discuss the implications for electron energy spectra from shock Fermi acceleration and discuss how the results may apply to solar flares.

## Ii Methods

For our simulation we use the MHD code ATHENA athenamethod (); athenamethod2 (), a cartesian, time explicit, unsplit, Godunov, code parallelized by MPI for compressible MHD. The mass conservation, momentum, and magnetic induction equation to be used are given by

(1) | |||||

(2) | |||||

where , is fluid pressure, is the magnetic field, is the mass density, is the velocity field, and the magnetic diffusivity equals the resistivity, i.e. . In the units used in ATHENA the magnetic permeability is unity and one converts to cgs units by dividing by .

In Eq. (II) , and we write as

(4) |

where is a uniform background resistivity and is an enhanced resistivity. The quantity determines the radial scale of the enhanced resistivity about the simulation origin surrounding the X-point. The functional form of the gaussian fall off of this enhancement is arbitrary and is simply chosen to provide the option to consider the effect of enhanced dissipation near the X-point for , the importance of which we will discuss later.

We set up initial isothermal equilibria for all cases but then distinguish cases in which perturbations are evolved adiabatically or isothermally. These cases are referred to as ”adiabatic” and ”isothermal” respectively. Note that by adiabatic we mean no heat loss from the system. Neither the adiabatic or isothermal cases are isentropic. As the isothermal and adiabatic simulations bracket a wide range of cooling efficiencies for the perturbed configurations motivates our inclusion of both cases.

For both adiabatic and isothermal simulations we initialize as

(5) |

where is the sound speed and is set equal to .5 initially and is initialized as discussed later in this section.

For our isothermal simulations no energy equation is solved and Eq. (5) becomes a fixed definition. For our adiabatic simulations we replace Eq. (5) with the equation for energy evolution given by

(6) |

where

(7) |

and evolves as with and evolves as evolves according to Eq. (6).

Our initial conditions follow those in the subset of the MHD GEM challenge simulationsBirn et al. (2001) with only resistivity and no viscosity. Technically, in MHD the gyro-radii and plasma particle inertial lengths are infinitesimal, and the non-ideal nature of the plasma is scalable entirely with the Lundquist or Magnetic Reynolds numbers given by where is the Alfvén speed.

It is tempting to normalize lengths to the ion inertial length and our time scale to the inverse ion gyrofrequency to see what kind of scales our choices of number density and magnetic field, and temperature would imply for a real system as well as to compare to the MHD and particle in cell simulations of the GEM ChallengeBirn et al. (2001). In our simulations, is the inflow axis and is the axis of the current sheet. We employ a rectangular domain with and , where and in arbitrary units. Upstream conditions of a particular solar flare tsuneta (), , and G, corresponds to an ion inertial length of cm, ion gyroradius of cm, and Alfvén speed so this indicates that that for this flare we (and the authors of other simulationsBirn et al. (2001)) would be assuming MHD is valid for scales cm. This is misleading however, since the presence or absence of shocks once in the MHD regime does not depend on the actual normalization of the density or the magnetic field. One need only be concerned that the MHD approximation is valid on scales for which a comparison to a real system is being made and to check the correpsonding Lundquist number.

Resolutions for the runs in Table 1 range from to . Most simulations were run at . The primary effect of resolution was to change the steepness of the flow profiles near the x point but not the final magnitude of the outflow at the simulation edges. Outflow boundary conditions were used at all edges for both fluid and magnetic quantities. The mass and magnetic field were not resupplied as the simulations evolved. Perturbation around a stable configuration causes the commencement of reconnection, which draws material toward the y-axis.

ATHENA maintains a zero divergence of the magnetic field if the initial divergence is zero. To ensure in ATHENA, we initialize the magnetic field via a vector potential such that . For our unperturbed Harris configuration we used

(8) |

where and for all simulations which then returned

(9) |

The density profile for our initial Harris configuration is

(10) |

where and for all simulations. Inclusion of a constant background density was used to prevent a large an Alfvén speed in low density regions that would result in an unreasonably small time step. The overall pressure balance and dynamics were left unchanged by the background density.

Using , makes , and and represents an initial equilibrium around which we peturb. For the cases we call ”adiabatic” we set the initial unperturbed in Eq. (6) to the isothermal equilibrium pressure, but then set . The deviation from equilibrium then evolves adiabatically according to Eq. (6).

To initiate reconnection, we use a vector potential perturbation

(11) |

where and for all simulations. The perturbed magnetic field is then

(12) |

Table 1 summarizes selected runs. The A, I, C, and E prefixes indicate adiabatic, isothermal, uniform resistivity (), and enhanced resistivity () respectively. Simulations IS and IL are equivalent to IE, but with half and twice the domain size respectively.

The prefix W indicates simulations in which a dense plasma ”wall” of zero temperature is injected at one end of the domain perpendicular to the direction between and to facilitate a shock when the outflow is super-magnetosonic. The wall is intended simply as a generic obstacle. In a solar flare a natural obstacle arises in standard geometries as the downward directed component of the outflow from the X-point impinges upon the closed soft X-ray magnetic loop of previously reconnected field linestsuneta ().

The wall is injected only once per simulation at time , but we insert the wall at different times for different runs to test the sensitivity to . The sound speed of wall density is set to zero to prevent the wall from expanding due to pressure gradients. The injected wall must supply a density contrast at least as high as that across any shock it facilitates. The wall density for the isothermal simulations was twice that of the adiabatic simulations, consistent with the fact that the mass conservation jump condition demands a higher compression ratio across isothermal shocks. The five wall simulations with enhanced resistivity led to stronger outflows than the two wall simulations with .

## Iii Results

Runs IC, IE, AC, & AE were unstable to reconnection and resulted in accelerated outflows along the x = 0 line. For these wall-free runs, Table 1 shows that at the final time , the uniform () vs. enhanced () resistivity cases produce Mach numbers (given by ) of 0.5 and 1.1 for the adiabatic simulations, and 1.1 and 2.9 for the isothermal simulations. Figure (1) and Figure (2) show the resulting flow profiles along the axis at times 42 and 50 and demonstrate that the flows had reached a steady state by . No flows were found to occur for test runs in which the perturbation to the magnetic field was not included.

The adiabatic simulations require to produce super-magnetosonic outflows but both adiabatic and isothermal cases produce higher outflow speeds with . These trends are expected: First, isothermal flows imply a release of thermal energy from the system so that the conversion of magnetic to bulk ram pressure in the momentum equation is more efficient for isothermal flows, leading to a higher Mach number. Heat is also released at isothermal shocks implying higher compression ratios compared to the adiabatic cases. Second, when the central resistivity is enhanced, Petschek type ”fast” reconnection ensues biskamp00 (); kulsrud01 () albeit for the artificial incompressible limit, a centrally enhanced resistivity facilitates the rapid rotation of field lines through the plasma near the X-point required for the Petschek configurationkulsrud01 (). The inflow magnetic field exerts a tension force parallel to the outflow that helps accelerate the latter. Slow shocks (not to be confused with the outflow fast shocks) appear in the inflow near the reconnection region for all our simulation cases where enhanced resistivity is used. This is shown in Figure (6) and Figure (7). Finally, lower mach number flows are expected in adiabatic simulations as the temperature of the plasma is increased which corresponds to an increase in the sound speed.

Although our present simulations are strictly MHD, since the ion-inertial length depends only on the plasma density, an enhanced resistivity would be expected for fully kinetic versions of our simulated parameter regime based on a comparison between the Sweet-Parker thickness and the ion inertial length using such parameters: For example, at of a typical (e.g. isothermal) simulation we calculate the Sweet-Parker thickness () where Using , inflow (and ) and we find in our units. Plasmas with (such as pre-flare solar loops krucker07 ()), are unstable to enhanced dissipation near the X-point biskamp00 (); kulsrud01 (); uzdensky03 (); uzdensky ().

Fig. 3 shows the vertical flow profile at for an isothermal EOS with varying combinations of and without the wall. Isothermal simulations were used as the lack of heating allowed for a clearer study of the effects of resistivity on the flow profile. The peak outflow speed toward the edge of the box and the compression ratios are quite insensitive to for , and insensitive to for .

Our main result is the formation of shocks in the wall simulations (those with W in Table 1), across which we measured compression ratios at along the y-axis. Fig. 4 shows snapshots of the vertical flow along at for adiabatic cases in which the wall was inserted at times, (cases AEW0, AEW30, & AW43 in Table 1) along with an adiabatic case with wall inserted at (ACW30). Only in cases with , did a shock form, and in all such cases maintained a compression ratio of . The final strength of the shock is relatively insensitive to the time the wall was inserted and the shock location converged by in all cases.

## Iv Potential Implications for Shock-Fermi Acceleration and Flares

In the limit that particle velocities are large compared to shock velocities and both are non-relativistic, the analytically computed accelerated particle energy spectrum from shock-Fermi acceleration is bell2 () , where . For our adiabatic shocks , so . (Note that a different proposed first order Fermi acceleration mechanism that also gives was proposed in dalpino (). There, in the context of relativistic particle acceleration in active galactic nuclei, Fermi acceleration is argued to occur as the oncoming reconnection inflow streams flow toward each other.)

For a solar flare, the index at emission sites is related to the photon number spectral index via , (where is the electron number spectral index) based on emission from electron-ion interactions at magnetic loop foot-pointsthe98 (). Electron energy spectral indices inferred from typical hard phases of flares krucker07 (); gb08 () are , consistent with the adiabatic prediction of from our simulations and consistent with scenarios that invoke fast shocks in reconnection down-flowsblackman97 (); tsuneta98 (); blackman2 ().

A typical flow travel time between the inferred X-point height and the loop top is seconds tsuneta (), and the time scales for conductive losses tsuneta () and bremsstrahlung cooling blackman2 () for keV electrons are comparable or within an order of magnitude less. Any cooling that does incur would make inferred at the emission site an upper limit to at the acceleration site. Note that the shock itself would be expected to be adiabatic because the shock thickness, determined by ion gyro-radii, is much thinner than the cooling scale.

Having chosen initial parameters to match the GEM challengeBirn et al. (2001) for consistency, we note that the distance bewteen the X-point and loop top obstacle in a typical flare is times larger than our simulation box given the box density if we were to artificially scale dimensions into units of the ion-inertial length based on our chosen density. However, our results are strictly MHD and so the fact that the real solar flare is much larger in such dimensionless units makes our conclusions that shocks are expected even more robustly relevant for the solar flare scales since the existence of shocks does not depend on the density as an independent parameter once in the MHD regime.

In addition to the mere presence of shocks and the expected power law index, it is noteworthy that the fast reconnection shocks in our simulations extend across the entire outflow width with very little deviation in the compression ratio across along the shock (see Fig. 6). All of the inflow then passes through the shock on its way out via mass conservation. The efficiency of electron acceleration at the shock would then not be a question of geometric limitation but a question of how efficiently the specific type of fast shocks accelerate electrons. This contrasts, for example, direct electric field electron acceleration along the x-line, which encounters only a limited electron throughput. A detailed study of the microphysics of shock acceleration at specifically weak fast shocks is an ongoing PIC-simulation research topic that is beyond the scope of the present work. The extent to which electrons must be ”pre-injected” above a certain energy threshold JandE () to be accelerated at such shocks is not entirely clear. More on this point below.

Observationally, solar flares have different populations of electrons and different paradigms have been proposed for the potential role that down stream fast shock acceleration might play blackman1 (); blackman2 (); tsuneta97 (); tsuneta98 (). Some flares have more non-thermal electrons by fraction than others alexander97 (); emslie03 (). In compact-impulsive flares, which have a large percentage of non-thermal electrons, there is a knee separating two power laws in the energy spectrum dulk92 () with the two power laws seemingly requiring different mechanisms blackman2 (): For the particles within the higher energy steeper power law, lower energy particles lag behind the higher energy particles, whereas, within the lower energy power law electrons, the high energy particles lag the lower energy particles asch95 (); asch97 (). In the shock reprocessing paradigm blackman2 (), the lower energy power law is taken to be the result of stochastic fermi acceleration in turbulent outflows from the reconnection X-point and and the highest energy particles are drawn from those which subsequently get injected into the shock as the flow passes through. The ”obstacle” causing the shock is a reconnected loop in the standard inverse-y picture of solar flares but the fraction of particles passing through the shock in the outflow region then depends the convexity of the loop; any pre-injection acceleration that might be needed for the particles to engage in the shock-Fermi acceleration could be accomplished by the stochastic acceleration region. In principle this can accommodateblackman2 () the needed rate of electron throughput in typical flares masuda94 ().

As we stated in the introduction, reconnection is best thought of as an acceleration environment and the relative importance of various mechanisms being possibly different in different flares. The discussion of the previous paragraph is one example of how such shocks may fit into such a paradigm. However, the main point of our paper here is to highlight the robust appearance of such shocks in simulations and their generic presence, regardless of whatever else may also be operating to produce observed phenomenology.

## V Conclusions

From 2-D MHD simulations, we have found that super-magneotsonic outflows are a generic feature of fast magnetic reconnection simulations and produce shocks upon encountering a plasma wall. We used the same initial parameters as those of the GEM challenge MHD simulationsBirn et al. (2001) to study this question. We find the fast shocks in our adiabatic cases (the most relevant for solar flares) to be weakly compressive with . Our results are generally consistent with previous analytic predictions of the existence of weak fast mode shocks in reconnection outflows blackman1 () and consistent the use of such shocks for particle acceleration in solar flares blackman97 (); tsuneta98 (); blackman2 (). The potential importance of such fast mode shocks and the complexity of magnetic reconnection warranted our explicit demonstration of their formation in a reconnection numerical simulation. It is also important to emphasize that the weak fast mode downstream shocks are distinct from the Petschek slow mode upstream shocks.

If Fermi acceleration is operative at such shocks, it will produce steeper particle energy spectra than maximally compressive adiabatic shocks. The resulting steeper particle spectral indices are not inconsistent with inferred distributions of solar flare electrons. An important next step is to directly study the microphysics of particle acceleration at these weakly compressive fast mode reconnection outflow shocks; most work has focused on particle acceleration at strong fast mode shocks.

We have not considered a guide field along the -axis in the third dimension and depending on on the spatial gradient of of such a field the Mach numbers of the outflow may change. The influence of a guide field would likely depend on its initial spatial gradient as to whether its net effect is to enhance or reduce the outflow pressure gradient. For example, if the guide field is very stronger at the X-point and weaker farther out, it could increase the outflow Mach number. If the field is uniform it would not produce a strong force. This is another topic for future work.

## References

- (1) M. Yamada, R. Kulsrud, & H. Ji Reviews of Modern Physics, 82, 603 (2010)
- Birn et al. (2001) J. F. Birn, M. A. Shay, B. N. Rogers, R. E. Denton, M. Hesse, M. Kuznetsova, Z. W. Ma, A. Bhattacharjee, A. Otto, & P. L . Pritchett, Journal of Geophysical Research, 106, 3715 (2001)
- (3) D. Biskamp, Magnetic Reconnection in Plasmas(Cambridge University Press, Cambridge, UK, 2000).
- (4) R. M. Kulsrud, E, P, & S, 417, 53 (2001)
- (5) D. Uzdensky ApJ 587 450 (2003)
- (6) N.F. Loureiro, A. A. Schekochihin, & S.C. Cowley,. Physics of Plasmas, 14, 100703 (2007)
- (7) G. Kowal, A. Lazarian, E.T. Vishniac, & K. Otmianowska-Mazur ApJ, 700, 63 (2009)
- (8) G.L. Eyink, A. Lazarian, & E.T. Vishniac, arXiv:1103.1882 (2011)
- (9) A. Bhattacharjee, ,Y.-M. Huang, H. Yang, & B. Rogers, Physics of Plasmas, 16, 112102 (2009)
- (10) M.M. Romanova, & R.V.E. Lovelace, Astron. Astrophys., 262, 26 (1992)
- (11) E. G. Blackman & G. B. Field, Phys Rev Letters, 73, 3097 (1994)
- (12) E.G. Blackman ApJ, 456, L87 (1996)
- (13) de Gouveia dal Pino, E. M., & Lazarian, A. 2005, Astron. Astrophys., 441, 845
- (14) G. Kowal, E.M. de Gouveia Dal Pino & A. Lazarian, arXiv:1103.2984 (2011)
- (15) E.G. Blackman,, ApJL 484 L79 (1997)
- (16) S. Tsuneta, S. Masuda, T. Kosugi, J. Sato, ApJ, 478, 787 (1997)
- (17) S. Tsuneta & T. Naito, ApJ, 495, L67 (1998)
- (18) H.E. Petschek, in Physics of Solar Flares edited by W. N. Ness, NASA SP-50, p. 425, (1964)
- (19) A.R. Bell, MNRAS 182, 443 (1978)
- (20) F. Jones & D. Ellison, Space Science Reviews, 58, 259 (1991)
- (21) R. Selkowitz & E. G. Blackman, MNRAS, 379, 43 (2007)
- (22) J. M. Stone, T. A. Gardiner, P. Teuben, J. F. Hawley, &J. B. Simon, ApJ Supplement, 178, 137 (2008)
- (23) K. Beckwith & J. M. Stone, ApJ Supplement, 193, 6B, (2011)
- (24) S. Tsuneta, ApJ, 456, 840 (1996)
- (25) S. Krucker, E.P. Kontar, S. Christe, S., & R.P. Lin, ApJL, 663, L109 (2007)
- (26) T.N. Larosa, R.L. Moore, J.A. Miller, S.N. Shore ApJ, 467, 454 (1996)
- (27) D.A. Uzdensky, Physical Review Letters 99, 261101 (2007)
- (28) E. Tandberg-Hanssen, & A.G. Emslie, Cambridge and New York, Cambridge University Press (1988)
- (29) P.C. Grigis & A.O. Benz, A. O. ApJ 683 1180 (2008)
- (30) D. Alexander & T.R. Metcalf ApJ, 489 442 (1997)
- (31) A.G. Emslie, ApJ, 595, L119 (2003)
- (32) G.A. Dulk, A.L. Kiplinger, R.M. Winglee R. M., 1992, ApJ, 389, 756
- (33) de Gouveia dal Pino, E. M., & Lazarian, A. 2005, Astron. & Astrophys, 441, 845
- (34) M.J. Aschwanden, R.A. Schwartz, D.M. Alt, ApJ, 447, 923 (1995)
- (35) M.J. Aschwanden, R.M. Bynum, T. Kosugi, H.S. Hudson, R.A. Schwartz, ApJ, 487 936 (1997)
- (36) S. Masuda.,T. Kosugi , H. Hara, S. Tsuneta, Y. Ogawara, Nat, 371, 495 (1994)

Sim. # | res. | |||||||
---|---|---|---|---|---|---|---|---|

AC | .01 | 0.0 | NA | NA | 75 | 0.46 | NA | |

AE | .01 | 10.0 | 0.25 | NA | 75 | 1.07 | NA | |

AEW0 | .008 | 10.0 | 0.25 | 0.0 | 50 | 1.85 | 1.9/2.0/2.2 | |

AEW30 | .008 | 10.0 | 0.25 | 30.0 | 50 | 1.97 | 1.8/1.8/2.1 | |

AEW43.5 | .008 | 10.0 | 0.25 | 43.5 | 50 | 1.79 | NA/2.1/2.1 | |

ACW30 | .008 | 0.0 | 0.25 | 30.0 | 50 | 0.35 | 0.0/0.0/0.0 | |

IC | .01 | 0.0 | NA | NA | 75 | 1.08 | NA | |

IE | .01 | 10.0 | 0.25 | NA | 75 | 2.91 | NA | |

IEW0 | .008 | 10.0 | 0.25 | 0.0 | 50 | 2.55 | 6.0/5.8/5.5 | |

IEW30 | .008 | 10.0 | 0.25 | 30.0 | 50 | 2.6 | 7.0/6.8/5.7 | |

IL | .008 | 10 | 0.25 | NA | 50 | 2.86 | NA | |

IS | .008 | 10 | 0.25 | NA | 50 | 2.85 | NA | |

IDCO | .016 | 0 | 0 | NA | 50 | 1.1 | NA | |

IDC | .016 | 10 | 0.25 | NA | 50 | 2.92 | NA | |

IDE | .008 | 20 | 0.5 | NA | 50 | 2.92 | NA | |

IDEC | .016 | 20 | 0.5 | NA | 50 | 2.83 | NA | |

IHEC | .008 | 2.5 | .125 | NA | 50 | 2.83 | NA |