Analysis of fast turbulent reconnection with self-consistent determination of turbulence timescale

TL;DR

This paper uses a Reynolds-averaged turbulence model to simulate magnetic reconnection, with a self-consistent turbulence timescale derived from the evolution of turbulent quantities, and compares the reconnection rates across different turbulence regimes.

Contribution

It introduces a self-consistent turbulence timescale calculation into magnetic reconnection modeling, improving upon algebraic models and enabling analysis of various turbulence regimes.

Findings

Reconnection rate varies with initial turbulent energy.

Self-consistent turbulence timescale influences reconnection dynamics.

Model aligns with resistive MHD rates in certain regimes.

Abstract

We present results of Reynolds-averaged turbulence model simulation on the problem of magnetic reconnection. In the model, in addition to the mean density, momentum, magnetic field, and energy equations, the evolution equations of the turbulent cross-helicity $W$, turbulent energy $K$ and its dissipation rate $\varepsilon$ are simultaneously solved to calculate the rate of magnetic reconnection for a Harris-type current sheet. In contrast to previous works based on algebraic modeling, the turbulence timescale is self-determined by the nonlinear evolutions of $K$ and $\varepsilon$, their ratio being a timescale. We compare the reconnection rate produced by our mean-field model to the resistive non-turbulent MHD rate. To test whether different regimes of reconnection are produced, we vary the initial strength of turbulent energy and study the effect on the amount of magnetic flux reconnected in time.

Equations

$$f=F+f^{\prime },F=⟨f⟩$$

$$f$$

$$F$$

$$f^{\prime }$$

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times (\mathbf{V}\times \mathbf{B}+\mathcal{E})+\eta {\nabla }^2\mathbf{B},$$

$$\mathcal{E}\equiv ⟨{\mathbf{v}}^{\prime }\times {\mathbf{b}}^{\prime }⟩.$$

$$\frac{D{\mathcal{E}}_i}{Dt}=$$

$$-⟨v_{\ell }^{\prime }{v}_k^{\prime }+b_{\ell }^{\prime }{b}_k^{\prime }/(\mu \overline{\rho })⟩\mu {ϵ}_{ijk}{\partial }_{\ell }{B}_j$$

$$+⟨v_{\ell }^{\prime }{b}_k^{\prime }+b_{\ell }^{\prime }{v}_k^{\prime }⟩{ϵ}_{ijk}{\partial }_{\ell }{V}_j+H.T.,$$

$$⟨{ϵ}_{ijk}{B}_{\ell }({\partial }_{\ell }{f}_j^{\prime })f_k^{\prime }⟩=\frac{1}{3}{\delta }_{\ell i}{B}_{\ell }⟨{ϵ}_{mjk}({\partial }_m{f}_j^{\prime })f_k^{\prime }⟩,$$

$$⟨f_i^{\prime }{g}_j^{\prime }⟩=\frac{1}{3}{\delta }_{ij}⟨f_{\ell }^{\prime }{g}_{\ell }^{\prime }⟩$$

$$\frac{D{\mathcal{E}}_i}{Dt}=$$

$$-\frac{1}{3}⟨v_k^{\prime }{v}_k^{\prime }+b_k^{\prime }{b}_k^{\prime }/(\mu \overline{\rho })⟩\mu {ϵ}_{i\ell j}{\partial }_{\ell }{B}_j$$

$$+\frac{2}{3}⟨v_k^{\prime }{b}_k^{\prime }⟩{ϵ}_{i\ell j}{\partial }_{\ell }{V}_j+H.T.$$

$$⟨{\mathbf{v}}^{\prime }\times {\mathbf{b}}^{\prime }⟩=\alpha \mathbf{B}-\beta \mu \mathbf{J}+\gamma \mathbf{\Omega }.$$

$$\alpha$$

$$\beta$$

$$\gamma$$

$$H$$

$$K$$

$$W$$

$$\alpha =\int \text{mbox}d\mathbf{k}\underset{-\infty }{\overset{t}{\int }}\text{mbox}d{\tau }^{\prime }G(k,\mathbf{x};\tau ,{\tau }^{\prime },t)$$

$$\beta =\int \text{mbox}d\mathbf{k}\underset{-\infty }{\overset{t}{\int }}\text{mbox}d{\tau }^{\prime }G(k,\mathbf{x};\tau ,{\tau }^{\prime },t)$$

$$\gamma =\int \text{mbox}d\mathbf{k}\underset{-\infty }{\overset{t}{\int }}\text{mbox}d{\tau }^{\prime }G(k,\mathbf{x};\tau ,{\tau }^{\prime },t)$$

$${\nu }_{\mathrm{T}}={\nu }_{\mathrm{T}0},$$

$${\nu }_{\mathrm{T}}\sim u\ell \sim u^2\tau ,$$

$${\nu }_{\mathrm{T}}=\frac{7}{15}\int d\mathbf{k}{\int }_{-\infty }^{t}d{\tau }_{1}G(\mathbf{k},\mathbf{x};\tau ,{\tau }_1)Q(\mathbf{k},\mathbf{x};\tau ,{\tau }_1,t),$$

$${\nu }_{\mathrm{T}}={\int }_{k_0}^{\infty }dk\tau (k)E(k),$$

$$E(k)=C_{\mathrm{K}}{\epsilon }^{2/3}{k}^{-5/3},$$

$$K={\int }_{k_0}^{\infty }E(k)dk={\int }_{k_0}^{\infty }{C}_{\mathrm{K}}{\epsilon }^{2/3}{k}^{-5/3}dk\sim {\epsilon }^{2/3}{k}_0^{-2/3}.$$

$${\tau }_{\mathrm{E}}\sim \ell /u\sim K/\epsilon .$$

$${\tau }_{\mathrm{A}}\sim (k{V}_{\mathrm{A}}{)}^{-1},$$

$${\tau }_{\mathrm{AE}}\sim ({\tau }_{\mathrm{E}}/{\tau }_{\mathrm{A}}){\tau }_{\mathrm{E}}.$$

$$E(k)=C_{\mathrm{IK}}(\epsilon V_{\mathrm{A}}{)}^{1/2}{k}^{-3/2},$$

$$K={\int }_{k_0}^{\infty }E(k)dk={\int }_{k_0}^{\infty }{C}_{\mathrm{IK}}(\epsilon V_{\mathrm{A}}{)}^{1/2}{k}^{-3/2}dk\sim (\epsilon V_{\mathrm{A}}{)}^{1/2}{k}_0^{-1/2},$$

$${\tau }_{\mathrm{A}}\simeq (k_0{V}_{\mathrm{A}}{)}^{-1}\simeq {\left(\frac{\epsilon V_{\mathrm{A}}}{K^2}{V}_{\mathrm{A}}\right)}^{-1}=\frac{K}{V_{\mathrm{A}}^2}\frac{K}{\epsilon }=\frac{K}{V_{\mathrm{A}}^2}{\tau }_{\mathrm{E}}.$$

$${\tau }_{\mathrm{A}}{V}_{\mathrm{A}}^2\sim {\tau }_{\mathrm{E}}K.$$

Full text

Turbulence – Magnetic reconnection – Subgrid-scale

††preprint: To be submitted to Phys. Plasmas††preprint: April 2019

I Introduction

The conversion of stored magnetic energy into different forms by magnetic reconnection is a key dynamical process for astrophysics and laboratory plasmas. In the sun, magnetic reconnection is believed to be the triggering process of solar flares as well as the mechanism responsible for heating the solar corona. Based on linear theory, the rate of energy conversion is given by the inverse magnetic Reynolds number $R_{m}$ (Sweet–Parker scaling) which is large for astrophysical plasmas. It turns out that the related reconnection rate is very small and can not explain the speed at which dynamical events in the solar corona take place. Several processes for collisionless plasmas have been proposed to increase the energy conversion rate.Schindler (2007); Büchner (2007); Silin, Büchner, and Vaivads (2005); Büchner and Elkina (2005); Terasawa (1983); Birn et al. (2001) For MHD, the plasmoid instability can increase the rate of magnetic reconnection.Loureiro, Schekochihin, and Cowley (2007); Uzdensky, Loureiro, and Schekochihin (2010); Bhattacharjee et al. (2009); Huang and Bhattacharjee (2015) However, turbulence is ubiquitous in astrophysical plasmas and may play a very important role for reaching fast magnetic reconnection.Matthaeus and Lamkin (1985); Lazarian and Vishniac (1999); Kowal et al. (2009) Lazarian and Vishniac (1999) proposed that the effective thickness of the diffusion layer is largely increased by the stochastic movement of the magnetic field lines. The argument has been shown not to hold with 3-dimensional simulations of stratified turbulence.Jabbari et al. (2016) A possible reason is that turbulence should be determined simultaneously with the evolution of mean-fields and not be externally imposed. In such situations, the statistical description of turbulence in mean-field theory is adequate/appropriate.

The range of astrophysical parameters for realistic numerical simulations is not reachable in a foreseeable future and turbulence has to be modeled. Even though several approaches to turbulence model exist,Krause and Rädler (1980); Sagaut (2001) we concentrate our work on a mean-field turbulence model. The model chosen is a Reynolds-averaged turbulence model where turbulence is self-generated and -sustained by inhomogeneities of the mean-fields.Yokoi, Higashimori, and Hoshino (2013) Through the model, mean-field MHD equations are solved together with transport equations modeling subgrid turbulent motions. Statistical quantities such as the turbulent cross-helicity (cross-correlation between magnetic-field and velocity-field fluctuations), the turbulent energy and its dissipation rate are the building blocks of the model. The dynamical balance between the turbulent energy and cross-helicity has been presented to produce fast reconnection.Yokoi, Higashimori, and Hoshino (2013); Yokoi (2019) Such mean-field model has already been confirmed in the context of magnetic reconnection through filtering procedure of high resolution simulations of plasmoid instability in MHD.Widmer, Büchner, and Yokoi (2016a) Other confirmations such as the expression of the electromotive force by direct numerical simulations (DNS)Yokoi and Balarac (2011) or the reproducibility of the Alfvénicity of the solar wind,Yokoi and Hamba (2007) strengthens the validity of the present mean-field approach of turbulence.

An important point in turbulence modeling is the determination of the turbulence timescale which is related to the closure scheme of the turbulence model under consideration. The original turbulence model for magnetic reconnection had been proposed in order that the timescale of turbulence may be self-consistently determined through the turbulence dynamics.Yokoi, Higashimori, and Hoshino (2013) However, as a first step, a simplified algebraic model for the turbulence timescale was numerically solved.Higashimori, Yokoi, and Hoshino (2013) There different turbulence timescales used as an external parameter led to different results. Depending on the initial value of the turbulence timescale, two regimes of energy conversion have been found to be laminar and fast turbulent while a third one was determined to be a turbulent diffusion not favorable to a situation of fast reconnection.Higashimori, Yokoi, and Hoshino (2013) These results have been confirmed and extended to force-free and Harris-type current sheets with out-of-plan guide magnetic field.Widmer, Büchner, and Yokoi (2016b) However, as noted in the previous work,Yokoi, Higashimori, and Hoshino (2013); Higashimori, Yokoi, and Hoshino (2013) the timescale of turbulence should be determined directly by the nonlinear dynamics of turbulence. Its spatiotemporal evolution should be determined either by the transport equation of timescale or by some quantities representing it. For instance, in the Kolmogorov picture of hydrodynamic turbulence, the ratio of the turbulent energy with its dissipation rate corresponds to a timescale. Solving simultaneously transport equations for the turbulent energy and its dissipation rate by numerical methods gives a spatiotemporal distribution of the turbulence timescale. In this paper, we solve such transport equations and present results confirming that fast reconnection is obtained by the model. On the other hand, the turbulent diffusive regime is not recovered even for large initial amplitude of turbulence.

II Modeling turbulence transport and timescale

II.1 Mean-Field Turbulence Model

In the Reynolds-averaged turbulence model, any instantaneous physical quantity $f$ is divided into a mean $F=\left<f\right>$ and a fluctuation around it, $f^{\prime}$, as

[TABLE]

with

[TABLE]

where $\boldsymbol{v}$ is the velocity, \mbox{\boldmath\omega}(=\nabla\times\boldsymbol{v}) the vorticity, $\boldsymbol{b}$ the magnetic field, $\boldsymbol{j}(=\nabla\times\boldsymbol{b})$ the electric-current density, $\rho$ the density, $p$ the pressure, $h$ the internal energy. Here, $\left<{\cdots}\right>$ is the ensemble average satisfying the Reynolds’ rules: $\langle{\langle{f}\rangle}\rangle=\langle{f}\rangle$, $\langle{f^{\prime}}\rangle=0$, $\langle{f\langle{f}\rangle}\rangle=\langle{f}\rangle\langle{f}\rangle$.Schmidt (2015); Sagaut (2001) Under the ensemble averaging, the mean-field induction equation reads

[TABLE]

where $\boldsymbol{\cal{E}}$ is the turbulent electromotive force (EMF) defined by

[TABLE]

The EMF solely represents the direct effects of turbulence on the mean induction equation. From the equations of the fluctuation fields, the evolution equation for $\boldsymbol{\cal{E}}$ is written as

[TABLE]

where $D/DT(\equiv\partial/\partial t+\boldsymbol{V}\cdot\nabla)$ is the Lagrange or material derivative, and $\rm{H.T.}$ stands for higher order terms whose detailed expressions are suppressed here.Yokoi (2013, 2018a)

If we approximate, for the sake of simplicity, that the inhomogeneity of the fluctuations $f^{\prime}=(\boldsymbol{v}^{\prime},\boldsymbol{b}^{\prime})$ along the mean magnetic field can be represented by the curls of $f^{\prime}$ as

[TABLE]

with $g=(\boldsymbol{v}^{\prime},\boldsymbol{b}^{\prime})$,

[TABLE]

This suggests that the turbulent electromotive force $\cal{E}$ is estimated as

[TABLE]

Here, the transport coefficients $\alpha$, $\beta$, and $\gamma$ are determined by the statistical properties of turbulence. Equation (7) suggests that they are expressed as

[TABLE]

where $\tau_{s}$ with $s=(\alpha,\beta,\gamma)$ is the timescales of turbulence, and the turbulent magnetohydrodynamic (MHD) energy $K$, the turbulent residual helicity $H$, and the turbulent cross helicity $W$ are defined by

[TABLE]

As we have just seen above, even with the simplest manipulation of the fluctuation equations we get a basic expression of the EMF [Eq. (8)] with its transport coefficients $\alpha$, $\beta$, and $\gamma$ [Eq. (9)]. In order to obtain more suitable expressions, we have to adopt more elaborated analysis of the MHD turbulence. With the aid of a closure theory of inhomogeneous MHD turbulence,Yoshizawa (1990); Yokoi (2013, 2018b) the theoretical analytical expression of $\boldsymbol{\cal{E}}$ [Eq. (4)] was investigated. On the basis of the theoretical results, the model expression for the EMF was proposed. In the model, the $\alpha$, $\beta$ and $\gamma$ coefficients are expressed in terms of Green functions which depend themselves on time.Yoshizawa (1990); Yokoi (2013, 2018b) A set of elaborated expressions for the transport coefficients of turbulence takes the form:

[TABLE]

where $G$ is the Green’s function of turbulence and $Q_{vv}$, $Q_{bb}$, $H_{vv}$, $H_{bb}$, $Q_{vb}$ and $Q_{bv}$ are the spectral functions in the wavenumber space. The derivation and physical meaning of Eq. (11) can be found in(Yokoi, 2013) As we will see in the following subsection, Eq. (11) is a natural generalization of the simplest expressions given by Eq. (9).

In view of Eqs. (9) and (11), it is obvious that we have to adopt appropriate timescales of turbulence in order to properly model the turbulent transport coefficients $\alpha$, $\beta$, and $\gamma$.

II.2 Turbulence transport and timescale

Depending on the level of the turbulence closure scheme, there are several ways to express the effective transport coefficients due to turbulence. Let us see this point, for example, in the case of the eddy or turbulent viscosity. The following arguments are based on a rather simple picture of homogeneous isotropic turbulence, but will provide a perspective on the relationship of the timescale of turbulence and turbulent transport.

The simplest way of describing the turbulent momentum transport may be just to treat the turbulent viscosity as a parameter. In this level of description, the eddy viscosity is given as

[TABLE]

with $\nu_{{\rm{T}0}}$ being a parameter, and the value of $\nu_{{\rm{T}}0}$ is determined in a heuristic manner.

The next higher level of description might be given by the so-called mixing-length modeling, where the transport coefficient is expressed by the characteristic velocity and length scales of turbulence, $u$ and $\ell$, as

[TABLE]

where $\tau$ is the characteristic timescale of turbulence and $\ell$ is called the mixing length, which is chosen in a heuristic manner based on the physical arguments. In this sense, the choice of $\ell$ is not directly based on the turbulence dynamics.

In the present work, we use the results by a more elaborated closure theory of inhomogeneous turbulence, the two-scale direct-interaction approximation (TSDIA). In this framework, the transport coefficients are expressed in terms of the Green’s and spectral functions. In this formulation, the eddy viscosity $\nu_{\rm{T}}$ is expressed as

[TABLE]

where $G({\bf{k}},{\bf{x}};\tau,\tau_{1})$, and $Q({\bf{k}},{\bf{x}};\tau,\tau_{1},t)$ are the Green’s and energy spectral functions, respectively. Since the Green’s function represents the weight how much the past states affect the present one, it basically gives information of timescale of turbulence. In a simpler case where the spectral function $Q$ is not affected by the past time, the time integral of the Green’s function is separately calculated. Then Eq. (14) is reduced to the spectral expression of the eddy viscosity

[TABLE]

where $k_{0}=(2\pi/\ell)$ is the wavenumber of the energy-containing eddy or mixing length $\ell$, $\tau$ is the timescale of turbulence, and $E(k)$ is the energy spectrum, showing how much energy is contained in the scale represented by wavenumber $k$. Note that if $\tau$ is represented by motions of the largest or energy-containing scales, and does not depend on scale, Eq. (15) reads to the simple mixing-length expression for the eddy viscosity [Eq. (13)]. From these arguments, we see that Eq. (14) is a much more elaborated expression capturing more complicated spatiotemporal properties of turbulence.

As we observed in Eqs. (12)-(14), whatever levels of description might be, the timescale of turbulence is directly related to the turbulent transport properties, and is one of the most important ingredients of turbulence modeling. There are several timescales in turbulence dynamics; the eddy distortion or turnover time, the spectral transfer time, the timescale associated with the decay of the triple correlations, etc. In the Kolomogorov’s picture for the homogeneous hydrodynamic turbulence, all those timescales are the same, and we substantially have only one timescale.

We assume the Kolmogorov’s scaling for the energy spectrum:

[TABLE]

where $C_{\rm{K}}$ is the Kolmogorov constant of the magnitude of $O(1)$ and $\varepsilon$ is the energy dissipation rate, which is equivalent to the energy flux or transfer rate in the wave-number space in the Kolmogorov picture of turbulence.

Substituting Eq. (16) into Eq. (15), we estimate the turbulent energy $K$ as

[TABLE]

With this estimate, the eddy turnover time $\tau_{\rm{E}}$ is expressed in terms of the turbulent energy ($K\sim u^{2}$) and its dissipation rate as

[TABLE]

Note that Eq. (18) is also obtained directly from the Kolmogorov’s four-fifth law $u\sim\delta u\sim(4/5)(\varepsilon\ell)^{1/3}$. In the Kolomogorov’s picture, the energy dissipation rate $\varepsilon$ is equivalent to the energy transfer rate: how much the turbulent energy fluxes from larger to smaller scales.

In the presence of the magnetic field, in addition to the eddy turnover time $\tau_{\rm{E}}$, the timescale associated with the Alfvén-wave propagation along the magnetic field may play an important role in turbulence evolution. This timescale, the Alfvén time $\tau_{\rm{A}}$, is defined as

[TABLE]

where $V_{\rm{A}}[=|\boldsymbol{b}|/(\mu\rho)^{1/2}]$ is the Alfvén speed of the magnetic field $\boldsymbol{b}$, and $k$ is the wavenumber ($\mu$: magnetic permeability, $\rho$: density).

In the Alfvén-wave turbulence, the energy transfer is caused solely by the interaction of two Alfvén-wave packets propagating in opposite directions along the magnetic field. Each interaction time is the Alfvén time $\tau_{\rm{A}}$, which is typically short compared with the eddy turnover time ($\tau_{\rm{A}}\ll\tau_{\rm{E}}$). In order to produce the same amount of energy transfer as in the hydrodynamic turbulence, where the eddy turnover time is $\tau_{\rm{E}}$, many ($\tau_{\rm{E}}/\tau_{\rm{A}}$ times) Alfvén-wave interaction events are needed. This suggests that the energy transfer time for the scale represented by the wavenumber $k$, $\tau_{\rm{AE}}$, in the Alfvén wave turbulence is longer than the hydrodynamic counterpart $\tau_{\rm{E}}$ as

[TABLE]

We assume the spectrum for the inertial range in Alfvén-wave turbulence is the Iroshnikov–Kraichnan type:

[TABLE]

where $C_{\rm{IK}}$ is the Iroshnikov–Kraichnan constant of the magnitude of $O(1)$. With this spectrum, the turbulent energy in Alfvénic turbulence is estimated as

[TABLE]

where $k_{0}$ is the characteristic wavenumber in the energy containing region or largest scale of turbulent motion. With this estimate, the Alfvén time can be connected to the eddy turnover time in terms of the turbulent energy and the mean magnetic-field energy as

[TABLE]

This should be compared with Eq. (20). Equation (23) suggests that the relative magnitude of the Alfvén time to the eddy turnover time is proportional to the ratio of the turbulent energy to the mean magnetic field

[TABLE]

If the turbulent energy $K$ is comparable with the mean magnetic energy represented by $V_{\rm{A}}^{2}$, the Alfvén time is comparable with the eddy turnover time $\tau_{\rm{E}}$:

[TABLE]

II.3 Synthesized Timescale

If the energy transfer is caused by both the processes of eddy turnover and Alfvén-wave interactions, we need to construct a synthesized timescale $\tau_{\rm{S}}$ from the eddy turnover time $\tau_{\rm{E}}$ and the Alfvén energy transfer time $\tau_{\rm{AE}}$. The representative one may be the timescale based on average:

[TABLE]

with the power index $p$, for which $p=1$ or $2$ is typically adopted. If we adopt $p=1$, corresponding to the harmonic average, it gives a synthesized timescaleYokoi et al. (2008):

[TABLE]

where use has been made of the relationship between $\tau_{\rm{A}}$ and $\tau_{\rm{AE}}$ [Eq. (20)] in the second estimate, and that between $\tau_{\rm{A}}$ and $\tau_{\rm{E}}$ [Eq. (24)] in the third estimate.

Equation (27) shows that in the weak magnetic-field limit ($V_{\rm{A}}^{2}\ll K$), the synthesized timescale recovers the usual eddy turnover time as

[TABLE]

On the other hand, in the strong magnetic-field limit ($V_{\rm{A}}^{2}\gg K$), the synthesized timescale is reduced to the Alfvén time as

[TABLE]

In Eq. (25) we saw that the Alfvén time is comparable to the eddy turnover time: $\tau_{\rm{A}}\sim\tau_{\rm{E}}$ if the turbulent energy is comparable to the mean magnetic field energy, $K\sim V_{\rm{A}}^{2}$. In this case, it follows from Eqs. (20), (23) and (25) that all the timescales; the synthesized, eddy turnover, Alfvén times, and Alfvén energy transfer time are comparable:

[TABLE]

This means that the synthesized energy transfer rate time $\tau_{\rm{S}}$ can be expressed by the eddy turnover time $\tau_{\rm{E}}(=K/\varepsilon)$, which is eventually the same as the Alfvén time $\tau_{\rm{A}}$.

In the simulation with a turbulence model (Sec. III.2), the mean-field equations are simultaneously solved with the turbulent-field equations. By solving the dissipation-rate equation ($\varepsilon$ equation) as well as the turbulent energy equation ($K$ equation), we obtain the spatiotemporal evolution of the timescale of turbulence by Eq. (27) or more simply by Eq. (18) or Eq. (30). This procedure provides a possibility of more self-consistent estimate of the turbulent transport than any turbulence model in which the timescale $\tau$ is just given as a parameter. Under the assumption that the eddy turnover and the Alfv n time are about the same and that the time integrals of the Green’s function can be separately computed from the turbulent field propagators, the transport coefficients $\alpha$, $\beta$ and $\gamma$ [Eq. (11)] are modeled asYoshizawa (1984); Yokoi and Hoshino (2011); Yokoi et al. (2016)

[TABLE]

where $C_{\beta}$, $C_{\gamma_{t}}$ and $C_{\alpha}$ are model constants and $\tau$ the timescale of turbulence evolution. As a simplified approach, we consider that the timescale of turbulence can be dynamically determined as

[TABLE]

The derivation and physical origins of the model expression Eq. (8) as well as Eqs. (31a)-(31c) and the timescale of turbulence can be found in Yokoi (2013).(Yokoi, 2013)

III Basic equations and numerical implementation

III.1 Mean-field equations

The GOEMHD3 code Skála et al. (2015) is used to solve the dimensionless mean-field MHD equations normalised to a mass density $\rho_{0}$, magnetic field strength $B_{0}$ and a typical length scale $L_{0}$ taken to be the half-width of the Harris-type current sheet.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\overline{\rho}$, $\boldsymbol{V}$, $\boldsymbol{B}$ and $\boldsymbol{J}$ are the mean mass density, fluid velocity, magnetic field and current density, respectively. The mean current density is calculated as $\boldsymbol{J}=\nabla\times\boldsymbol{B}/\mu$. The symbol $\overline{h}$ represents the internal energy of the mean-fields and is related to the thermal pressure by the equation of state $P=2h^{\gamma}_{0}$ for adiabatic conditions ($\gamma_{0}=C_{V}/C_{p}=5/3$ is the ratio of specific heats).Widmer, Büchner, and Yokoi (2016a) The resistivity is given by $\eta=10^{-4}$ and the parameter $\chi$ describes the amount of local smoothness used to avoid numerical instabilities. It is defined as $\chi=\chi_{0}+\chi_{\rm{loc}}$ where $\chi_{0}=10^{-3}$ is uniform and constant while $\chi_{\rm{loc}}\ll\chi_{0}$ is activated locally if necessary.

The turbulent electromotive force $\langle{{\boldsymbol{v}}^{\prime}\times{\boldsymbol{b}}}^{\prime}\rangle$ as well as $\varepsilon$ represents the feedback of MHD turbulence in the mean induction equation. As a first step, the turbulent stress-tensor $\langle{\boldsymbol{v}^{\prime}\otimes\boldsymbol{v}^{\prime}-\boldsymbol{b}^{\prime}\otimes\boldsymbol{b}^{\prime}}\rangle$ is neglected in the momentum equation (34) and only the turbulent electromotive force relates turbulence to the mean fields.

III.2 Turbulence equations

In addition to the system of mean-field MHD Eqs. (33)-(36), the governing equation for the turbulent energy $K$, its dissipation rate $\varepsilon$ and the turbulent cross-helicity $W$ are implemented in the GOEMHD3 code and solved as

[TABLE]

[TABLE]

[TABLE]

where $\boldsymbol{\cal{E}}$ is the turbulent electromotive force [Eq. (8)] and $C_{\varepsilon n}(n=1-3)$ model constants of order $\mathcal{O}(1)$ and are chosen as $C_{\varepsilon 1}=1.4$, $C_{\varepsilon 2}=1.9$ and $C_{\varepsilon 3}=1.0$.Yoshizawa (1990) The symbol $P_{K}$ is the turbulent energy production mechanism given as

[TABLE]

The purpose of solving Eq. (39) in addition to Eqs. (37) and (38) and the mean-field MHD Eqs. (33)-(36) is to determine the timescale of turbulence accordingly to the turbulence dynamics. In such a way, we can test whether turbulence produces different regimes of reconnection as obtained in previous works using a constant-$\tau$ approximation.Higashimori, Yokoi, and Hoshino (2013); Widmer, Büchner, and Yokoi (2016b) Since such regimes of reconnection were obtained independently of the initial current sheet configuration used, we only consider a Harris-type current sheet. The turbulent helicity related term $\alpha$ is neglected in the present work because there is no symmetry breakage capable to generate it in a Harris-type current sheet without guide magnetic field.Widmer, Büchner, and Yokoi (2016a)

The geometry is defined by the unit vectors $\bf{e}_{y}$ (across the current sheets), $\bf{e}_{z}$ (along the current sheets) and $\bf{e}_{x}$ (in-plan direction). The resolution is $4\times 2048\times 2048$ for a box size of $0.4L_{x}/L_{0}\times 80L_{y}/L_{0}\times 80L_{z}/L_{0}$ where the normalizing length $L_{0}$ is the current sheet half-width. To use periodic boundary conditions, a pair of Harris-type current sheet are initialized as

[TABLE]

for an asymptotic value of the magnetic field $B_{0}$=1. Both current sheets are located at $\pm 10L_{y}/L_{0}$. Reconnection is triggered by a divergence free perturbation

[TABLE]

with $b_{p}/b_{0}=10^{-3}$. The initial zero flow condition ($\boldsymbol{V}=0$) gives $W_{0}=0$ as an initial condition for the cross-helicity. Hence, the initial balance for turbulence is given by

[TABLE]

Here $P_{K_{0}}$ is the initial turbulence energy production given by

[TABLE]

Equation (43) provides an initial condition for the turbulent dissipation rate. It depends only on the initial value of the turbulence energy $K_{0}$ which is varied. The initial strength of the turbulence energy $K_{0}$ and its dissipation rate $\varepsilon_{0}$ are given as

[TABLE]

[TABLE]

where $\mu_{0}\boldsymbol{J}_{0}=\nabla\times{\boldsymbol{B}}_{0}$.

The turbulent timescale $\tau$ is regularized to avoid numerical errors asMaget et al. (2018)

[TABLE]

where $M$ and $m$ are the maximum and minimum values that the timescale of turbulence $\tau$ can attain. Since the turbulent timescale is normalised to the Alfv n time $\tau_{A}$ in our simulations, both limiters avoid unphysical values of $\tau$. Figure 2 represents Eq. (47) for different values of $M$ and $m$. In previous work, the algebraic timescale was considered through the parameter $C_{t}=\tau_{0}/\tau$ where $\tau_{0}=C_{\beta}^{-1/2}|\boldsymbol{J}|_{t=0,z=0}^{-1}$. We choose $m=0.06$ and $M=10$ to let $\tau$ the possibility to span reconnection regimes from laminar to turbulent as well as turbulent diffusion obtained within the range $C_{t}\in[0.05,3]$.Higashimori, Yokoi, and Hoshino (2013); Widmer, Büchner, and Yokoi (2016b)

IV Fast Reconnection

IV.1 Reconnection rate

The simplest theory of reconnection is the Sweet–Parker theory based on the conservation of magnetic flux and matter.Parker (1957) The model, developed for resistive-MHD (hereafter $\eta$-MHD), assumes that two regions exist around the current sheet. Magnetic reconnection is then estimated at the time the current sheet has reached a steady state, i.e., when the diffusive term ($\eta$) of the induction equation balances the convective term ($\boldsymbol{V}\times\boldsymbol{B}$). The estimated reconnection rate ($M_{\rm{SP}}$) is found to be proportional to the inverse of the square root of the magnetic Reynolds number

[TABLE]

with the resistivity $\eta$, the typical length scale $L$, and the plasma velocity $V$. This simple form is not a good approximation for large-magnetic-Reynolds-number astrophysical plasma. Using our mean-field turbulence model, the magnetic reconnection rate can be approximated by replacing $\eta\rightarrow\eta+\beta$ providing a turbulent magnetic Reynolds number $(Rm^{\rm{T}})^{-1/2}\equiv\sqrt{(\eta+\beta)/(VL)}$ for the Sweet–Parker model. Such a simple estimation shows that for collisionless astrophysical (large-magnetic-Reynolds-numbers) plasmas, turbulence might enhance the magnetic reconnection rate by a factor $\sqrt{1+\beta/\eta}$ above the $\eta$-MHD estimation.

Figure 1 depicts the spatial distributions of the mean current density $\boldsymbol{J}$, the mean vorticity $\boldsymbol{\Omega}$, the turbulence cross-helicity $W$, the turbulence energy $K$, its dissipation rate $\varepsilon$, and the turbulence timescale $\tau$ estimated from the evolution of $K$ and $\varepsilon$. We found that the turbulence energy $K(\propto\beta)$ is acting as a localized anomalous resistivity located at and around the ‘X’-point. Also, the turbulence energy dissipation rate $\varepsilon$ is finite at and around the diffusion region where the turbulence energy $K$ is maximum. The location where $\varepsilon$ is finite also represents the region in which the large scale magnetic field energy is transported to smaller scales where it can be more effectively dissipated. The timescale of turbulence $\tau$ is therefore maximum near the diffusion region where magnetic reconnection takes place. The intensity of $\tau$ is in the range $[1.2;1.4]$ corresponding to the regimes of fast turbulent reconnection obtained for a constant turbulence timescale.Higashimori, Yokoi, and Hoshino (2013); Widmer, Büchner, and Yokoi (2016b)

In previous works, a regime of slow energy conversion similar to the $\eta$-MHD rate (laminar reconnection) and that of turbulent diffusion with slower rate of energy conversion (turbulent diffusion) were obtained as well as the regime of fast reconnection (turbulent reconnection).Higashimori, Yokoi, and Hoshino (2013); Widmer, Büchner, and Yokoi (2016b) These regimes were consequences of variation of an adjustable timescale parameter $\tau$ or the initial turbulence energy intensity $K_{0}$ as being independent parameters. In this work, only the initial intensity of the turbulence energy $K_{\rm{init}}$ is varied because the counterpart of the energy dissipation rate, $\varepsilon_{0}$, is directly related to $K_{0}$ by Eq. (46). The model constants $C_{\beta}$, $C_{\gamma}$, $C_{n}(n=1-3)$ obtained from the turbulence modeling are not varied. Also, the current-density intensity is determined by the Harris-type current sheet initialization.

Figure 3 presents the time evolution of the reconnection rate as $K_{\rm{init}}$ is varied for a resistivity $\eta=10^{-4}$. Turbulence produces fast reconnection in comparison to the $\eta$-MHD regime.

Because the turbulence diffusivity ($\beta$) is localized, the Petschek model of reconnection describing a localized diffusion region due to standing shock wavesPetschek (1964) is a more appropriate estimation of the magnetic reconnection rate. According to the Petschek model, the reconnection rate ($M_{\rm{P}}$) is expected to be proportional to the inverse of the logarithm of the magnetic Reynolds number as

[TABLE]

Utilizing Eq. (49), the turbulent reconnection rate ($M_{\rm{T}}$) may be estimated as

[TABLE]

with the intensity of $K$, $\varepsilon$, and $\tau=K/\varepsilon$ being taken at the reconnection peak (Fig. 1).

Figure 4 represents Eq. (50) for various $K_{0}$. A turbulent energy intensity of $K=0.15$ enhances the $\eta-$MHD reconnection rate by a factor of about $3.2$ independently of the initial turbulence intensity $K_{0}$. Figure 5 shows time history of the turbulent reconnection rate to the $\eta-$MHD ratio. The ratio $M_{T}/M_{P}$ is about $3.2$ at the time both regime reach their saturation state. Turbulence enhances the reconnection rate for large magnetic Reynolds number above the $\eta-$MHD value as a Petschek-like reconnection.

In marked contrast to the previous results obtained from the simulations with a constant timescale parameter $\tau$, only fast turbulent reconnection is obtained. In fact, the other two regimes of energy conversion, laminar reconnection and turbulent diffusion, are not reproduced in the present work, where the timescale of turbulence is self-consistently solved within the mean-field MHD equations. In this sense, the other two regimes, laminar reconnection and turbulent diffusion, are artifacts arising from the turbulence timescale as a constant parameter. Figure 6 depicts the production $P_{A,i}$, transport $T_{A}$, advective $\boldsymbol{V}\cdot\nabla A$ and dissipation $\epsilon_{A}$ terms in Eqs. (37)-(39) at $t/t_{A}=120$. The reconnection rate initiates its saturation regime at this time. Here $A$ refers to $K$, $W$, or $\varepsilon$ and $i=1,2$. The turbulent energy production term $P_{K,1}=C_{\beta}\tau K\boldsymbol{J}^{2}$ is located at the current sheet center and is mostly balanced by its dissipation term $\varepsilon_{K}$, a similar behavior is obtained for the turbulent energy dissipation rate $\varepsilon$. The evolution of turbulence reaches a steady state resulting in a saturation of the reconnected magnetic flux amount.

V Remark on Strong Compressibility Limit

The strong localisation of the turbulent energy diffusive-like term $\beta$ around the reconnection region produces a Petschek-like reconnection process.Petschek (1964) In such a situation, slow shock-waves contributing to the fast Petschek-like reconnection are usually associated with large density variance at the reconnection region. We thus consider a contribution of the density variance in the electromotive force $\cal{E}$ (Eq. (8))Yokoi (2018b, a)

[TABLE]

where $\chi_{\bar{\rho}}$ is the transport coefficient associated with the density variance $\langle{\rho^{\prime}{}^{2}}\rangle$ representing genuine compressibility effects. It is important to stress that the last term in Eq. (51) has no counterparts in the incompressible case. The derivation of the electromotive force in a strong compressible regime as well as the physical origins of the $\chi_{\bar{\rho}}$ related term can be found in Yokoi 2018.Yokoi (2018b, a) There, it is shown that the $\bf{B}\times\nabla\bar{\rho}$ direction contributes to the turbulent electromotive force, enhancing the turbulent energy production across the reconnection slow shock front. Such a phenomena might increase magnetic reconnection. A simple consideration of a Sweet–Parker current sheet for the electromotive force Eq. (51), neglecting the $\alpha$ related term, lead to the following reconnection rate

[TABLE]

the hats denoting normalised variables . Equation (52) is represented in Fig. 7 for $\hat{\chi}_{\bar{\rho}}=0$, $0.25$. The red dot represents the values of $\beta$ and $\gamma$ at $t=120\tau_{A}$ and the white line represents the limit of turbulence saturation for fast reconnection $|\gamma|/\beta=1$.Yokoi and Hoshino (2011) We see that this simple estimate underestimate the reconnection rate obtained through our full-compressible MHD numerical simulations for $\chi_{\bar{\rho}}=0$ but is in a similar range when $\chi_{\bar{\rho}}=0.1$. Including the density variance effect in the turbulent electromotive force $\cal{E}$ enhances the reconnection rate. Further numerical validations based either on Direct Numerical Simulations (DNSs) or turbulent mean-fields simulations evolving evolution equations for turbulent transport coefficients, as done previously,Widmer, Büchner, and Yokoi (2016a); Higashimori, Yokoi, and Hoshino (2013); Widmer, Büchner, and Yokoi (2016b) are required.

VI Conclusions

In this work, we have performed numerical simulations of a turbulence model for magnetic reconnection with an equation of the turbulent energy dissipation rate ($\varepsilon$ equation) implemented.Yokoi, Higashimori, and Hoshino (2013) In the simulations, where the timescale of turbulence is evaluated from the turbulent energy $K$ and its dissipation rate $\varepsilon$, we obtained only the regime of fast reconnection among the three regimes of the turbulent reconnection. From the viewpoint of the mean-field turbulence model, this is quite a natural result. The turbulence timescale is determined by the nonlinear dynamics through the coupling of the evolutions of the turbulence energy $K$ and its dissipation rate $\varepsilon$. The contribution of the production and dissipation terms in the evolution equations of $K$ and $\varepsilon$ are competing. They mostly balance close to the time of reconnection rate saturation, turbulence production and dissipation reach a steady state and cannot enhance the reconnection process further.

We found that the magnitude of the turbulence energy $K$ is of the same order as the counterpart of the turbulent dissipation rate $\varepsilon$, and that the timescale provided by the ratio $\tau=K/\varepsilon$ is about twice unity in the diffusion region. The value of $\tau$ corresponds to the time needed for an Alfv n wave to cross the initial current sheet width ($\tau_{\rm{A}}\approx\tau_{\rm{AE}}$). In such situations, the eddy-turnover time can be used to estimate the timescale of turbulence as $\tau_{\rm{E}}=K/\varepsilon$. As this result, the rate of energy conversion of our mean-field turbulence model is enhanced above the rate of energy conversion of $\eta$-MHD.

Acknowledgements.

One of the authors (FW) acknowledges the International Max-Planck Research School (IMPRS) at the University of Göttingen as well as the CRC 963 project A15. JB thanks the Max-Planck–Princeton Center for Plasma Physics for its support. A part of this work was performed during the period when one of the authors (NY) visited Goslar and Göttingen in August 2017, and under the support of the JSPS Grants-in-Aid for Scientific Research 18H01212.