Magnetic reconnection-driven turbulence and turbulent reconnection acceleration

The MHD-PIC module embedded in PLUTO code can simulate the evolution of fluids by solving the set of MHD equations and simultaneously integrate the motion of CR particles using conventional PIC techniques. Particle back-reaction on the fluid is included in the form of momentum–energy feedback, and the CR-induced Hall term in Ohm’s law was introduced by Bai et al. (2015). The structure of this module mainly included three parts: solution of the MHD equations, integral of the motion equation of particles, and stepping time scheme. In general, this module, dealing with microphysics such as the resistance and Hall effect in the MHD part, is appropriate to capture the dynamical evolution of a plasma consisting of a thermal fluid and a non-thermal component represented by relativistic charged particles. Since the gyroradius of CRs (typically representing energetic protons) greatly exceeds the ion skin depth of the plasma, this module extends the range of applicability to much larger spatial (and longer temporal) scales.

The MHD-PIC module adopts finite-volume numerical schemes to solve the set of MHD equations on a computational grid, while the injected computational particles represent clouds of physical particles that are close to each other in phase space (Mignone et al. 2018). The divergence-free condition is enforced using either constrained transport or divergence-cleaning approaches. In addition, particles are treated kinetically either using conventional PIC time-reversible integrators or simple Runge-Kutta schemes. We notice that this module has been successfully tested for Bell instability and shock acceleration, and the results are in good agreement with previous studies based on Athena code (Bai et al. 2015). In the parallel computing tests, the CPU utilization efficiency can be stable at over 0.8 (Mignone et al. 2018 for more details).

To study turbulent magnetic reconnection processes, we perform numerical simulations employing the MHD-PIC module (Mignone et al., 2018) mentioned above to solve the dimensionless ideal MHD equations as follows

These equations are the continuity, momentum, energy, induction, and solenoidal condition, respectively. Here, $\rho$ is the mass density, $\bm{m}=\rho\bm{v}_{\rm g}$ the momentum density, $I$ unit tensor, $e$ specific internal energy, $\bm{v}_{\rm g}$ the gas velocity, $p$ the gas pressure, $\bm{B}$ magnetic field, and $E_{\rm t}=\rho e+\bm{m}^{2}/2\rho+\bm{B}^{2}/2$ the total energy density. The evolution of the magnetic field is governed by Faraday’s law of $\bm{E}=-\bm{v}_{\rm g}\times\bm{B}$. It should be stressed that we do not include resistant dissipation in our simulation.

Our numerical simulations are performed in a 3D domain with physical dimensions $1.0L\times 0.5L\times 1.0L$, where the dimensionless length scale is set to $L=C/\omega_{\rm p}=10^{4}$ with the light speed of $C=10^{4}V_{\rm A}$ and the plasma frequency of $\omega_{\rm p}=\sqrt{4\pi\rho q^{2}/m_{\rm p}}$ ($q$, $m_{\rm p}$ are a charge and mass of the proton, respectively). The configuration of the initial magnetic field is considered as a Harris type by (Harris 1962)

where $w$ is the initial width of the current sheet, and $B_{0}$ is the magnetic field strength controlled by parameter $\sigma=B_{0}^{2}/(4\pi\rho)$. The magnetic fields are antiparallel to the $X$-direction, resulting in a discontinuity placed on the $X-Z$ plane. By counteracting the Lorenz-force term with a thermal pressure gradient, we can obtain an initial equilibrium condition

where the plasma parameter is set as $\beta=0.01$ (Puzzoni et al., 2021). As a result, the total pressure remains invariable throughout the whole current sheet.

The evolution of CR particles can be described by the following equations

where $\bm{x}_{\rm p}$ and $\bm{v}_{\rm p}$ indicate the spatial position and velocity of a charged particle (proton for our scenario), respectively. The parameter $\alpha_{\rm p}$ in Equation (9) denotes the particle charge to mass ratio. When numerically solving Equations (8) and (9) by the Boris integrator (Boris & Shanny 1970), the electric and magnetic fields $\bm{E}$ and $\bm{B}$ are computed from the magnetized fluid by a cubic spline interpolation approach.

We uniformly inject $10^{6}$ particles into the whole box space and initialize them with a Maxwellian distribution of the thermal velocity of $0.1V_{\rm A}$ in each direction. The test particles are evolved together with the fluid using the Boris algorithm (Mignone et al. 2018). We use periodic boundaries in the $X$ and $Z$ directions and reflective conditions in the $Y$ direction. To numerically solve Equations (1) to (5), we choose the HLL Riemann solver with the Characteristic Tracing Contact (CT-Contact) electromagnetic fields averaging scheme and the second-order piecewise linear reconstruction and 2rd-order Runge Kutta time stepping. Except for the part on comparative studies of the resolution, we set the numerical resolution to $512~{}\times~{}256~{}\times~{}512$ through this paper, leading to the grid size of $\delta L\simeq 20$. In our simulations, we set some fixed parameters such as the uniform background plasma density of $\rho=1.0$, light speed of $C=10^{4}$, and Alfvén velocity of $V_{\rm A}=1.0$. The variable parameters are listed in Table 1 where we set the initial velocity perturbation with the amplitude of $V_{\rm eps}\sim 0.1V_{\rm A}$ to promote the generation of turbulence, that is, to shorten the simulation time. When the spectral energy distributions of the test particles reach a statistically steady state, we terminate the simulation at the final integration time of $t=7\times 10^{4}\omega_{\rm p}^{-1}$.

We run four models with the initial parameters listed in Table 1 by adding the stochastic noise of velocity fluctuations to launch a magnetic reconnection process. We first explore how the current sheet evolves during the spontaneously driven turbulent reconnection process. As an example, Figure 1 shows the current density slices on the $Y-Z$ plane at the fixed $X=0$ for the selected three snapshots $t=4~{}000,10~{}000$ and $70~{}000$ in units of the code. We find that the current density with noise structure was uniformly distributed within a thin current sheet at the early stage of evolution and then begins to form high-density clumps ($t<4~{}000$). When the fluid evolves to $t\sim 10~{}000$, we see prominent tearing of the current sheet. There are some magnetic islands-like structures similar to 2D magnetic reconnection simulation (e.g., Loureiro et al. 2007; Mignone et al. 2018; Puzzoni et al. 2021). This is because the initial current sheet with noise structure provides an instability incentive, which develops various fluid instabilities (e.g., Kelvin-Helmholtz and tearing instabilities) that make the current sheet fragment. Kowal et al. (2020) found that tearing instability dominates turbulence generation in the early stages of the evolution, and due to the suppression of tearing instabilities by magnetic field components perpendicular to the current sheet, the dominant mechanism in the later stages shifts to Kelvin-Helmholtz instabilities. Although quantitative analysis of instability is beyond the scope of our paper, we think that similar mechanisms, i.e., tearing and Kelvin-Helmholtz instabilities, should work at the hybrid scales between MHD and plasma scales, which needs to be further confirmed. As the fluid further evolves, the current sheet begins to thicken and form a significant large-scale magnetic flux rope structure at $t=70~{}000$. In our simulation, the setup of reflective boundary conditions in the $Y$-direction will cause the magnetic flux to reflect at the boundaries. The rebounded magnetic flux enters the reconnection region and intertwines with the magnetic flux in the reconnection region, thus forming a large-scale magnetic flux rope structure. This is different from the case of outflow boundary conditions applied in the direction of the perpendicular current sheet plane, in which such large-scale interactions of magnetic fluxes are not allowed (see Kowal et al. 2017).

The structure of the current sheet can also be observed from the 3D perspective in Figure 2 which shows a 3D visualization of the logarithm of the current density obtained from M2 at $t=70~{}000$ in units of the code. It is noticed that the current sheet on the $Y-Z$ plane has obvious tearing, while only thickening on the $X-Y$ plane. The reason is that the magnetic field direction is set along the $X$-axis direction. In addition, due to magnetic tension, the tearing of the current sheet is suppressed on the $X-Y$ plane.

Figure 3 shows the change of magnetic and kinetic energies within the reconnection region over the evolution time. As shown in the top panel of Figure 3, the continuously decreased magnetic energy with magnetic reconnection evolution demonstrates the presence of the release of magnetic energy by reconnecting. For model M2, the reduction of magnetic energy always lags behind other models, because its low initial noise added in velocity fluctuations makes its turbulent reconnection slightly slower than other models. Compared with models M3 and M4, the reduction in magnetic energy from model M4 is more significant, especially in the later stage ($t\gtrsim 40~{}000$). This is because the guiding field makes the magnetic field lines on both sides of the current sheet no longer strictly parallel so that they have a certain angle, which tends to enhance the reconnection efficiency. Compared with models M1 and M4, the release of magnetic energy in model M1 is more efficient. Although, this is consistent with conventional thinking: a large $\sigma$ value, i.e., a strong magnetic field, leads to the release of more reconnected magnetic energy. However, the more fundamental reason is that instabilities of initial noise excitation make the current sheet thicker (seems to touch near the boundaries) and tear more in the case of a strong magnetic field, as shown in the right upper panel of Figure 1.

The evolution of fluid kinetic energy over time first experiences a short plateau period then increases rapidly to form a sharp peak, and finally gradually tends to a quasi-steady state. In the initial stage of the beginning of reconnection, the release efficiency of magnetic energy is low and stochastic noise needs time to excite instabilities that induce the generation of turbulence. Subsequently, the magnetic energy is rapidly released, and the kinetic energy is suddenly increased. When the increase and dissipation of fluid kinetic energy reach a balance, the fluid kinetic energy reaches a peak. With the development of magnetic reconnection, a part of the dissipated magnetic energy is converted into the kinetic energy of the fluid. As shown in the bottom panel of Figure 3, at the M1 peak position, the magnetic energy decreases by about 0.2, while the fluid kinetic energy increases by about 0.13, and the remaining energy is converted into the kinetic energy of test particles. Throughout evolution, the average increase in kinetic energy hardly exceeds 0.1, while the continuous release of magnetic energy can reach about 0.6 ($=1-0.4$). Furthermore, since there are a few parts of particles in the current sheet at the initial stage of evolution, most of the magnetic energy released by the reconnecting is converted into kinetic energy. As the current sheet expands over time, the number of particles in the current sheet will increase, resulting in the direct conversion of magnetic energy into test particle kinetic energy.

Comparing these four models, we find that model M2 has a lag, which is consistent with the evolution of magnetic energy, but the lag phenomenon is more prominent. Model M4 has a higher peak value and a faster increase in kinetic energy than M3. The reason is that it has a significantly faster decrease in magnetic energy. Model M1 presents a more efficient increase of kinetic energy because of its stronger magnetic field. Both deformations of the current sheet and the formation of magnetic flux ropes can lead to changes in the strength of the outflow from the local reconnection region and interact with the developing turbulence. Turbulence increases when reconnection promotes this interaction. In turn, developing turbulence affects the reconnection process, leading to an increase in reconnection rates. With the evolution of time, the region around the current sheet is also successfully turbulent, providing an opportunity for the particles escaping the current sheet or the local particles to interact with turbulence.

We show the power spectra at the different evolution snapshots in Figure 4, including the magnetic fields (left column), velocities (middle), and densities (right) for the model M1 (upper panels) and M3 (lower). These snapshots correspond to $t=10~{}000$ to $70~{}000$ at intervals of $10~{}000$, where for comparison we also include the power spectrum at the initial stage of evolution of $t=2~{}000$. Since turbulence just begins to develop, the power spectrum fluctuates greatly for magnetic fields, velocities, and densities. With the evolution of time, the reconnection rate gradually increases, and the fully evolved power spectrum of MHD turbulence gradually converges. In the case of the velocity and density, they all well present the scaling of $k^{-5/3}$, which is the same as the Kolmogorov (1941) spectrum. This demonstrates that reconnection-driven turbulence also provides a scaling of $-5/3$ similar to that of Kolmogorov and GS95. However, the magnetic field presents the scaling of $k^{-2.2}$ as shown in the left panels, which is steeper than the Kolmogorov spectrum. This numerical experimental result implies that magnetic and velocity vector fields do not couple well. In our work, we also calculate the power spectra from models M2 and M4 (not shown) and find no significant differences between them.

Given the consistency of the velocity spectra with the generally accepted GS95 model, we analyze the anisotropic property of velocity in depth. We plot in Figure 5 the contour maps of the second-order structure function of velocity in a local reference frame for four models explored at the final snapshot ($t=70~{}000$). We can see that in all of these four models, it presents a well-anisotropic property. The semi-major axis is along the $\ell_{\parallel}$ direction, which is parallel to the local mean magnetic field. In addition, the models M1 and M4 show slightly stronger anisotropy than M2 and M3. This may be because the two formers correspond to a strong mean magnetic field and an additional guiding field, respectively, both of which have stronger magnetic field effects.

To quantify the anisotropy of reconnection-driven turbulence, we also show the anisotropy scaling of velocity in Figure 6, arising from the early ($t=30~{}000$) and late ($t=70~{}000$) snapshots. For the early snapshot (upper panel), the models M1 and M4 are very similar, both of which will satisfy the GS95 scale-dependent anisotropy of $\ell_{\parallel}\propto\ell_{\perp}^{2/3}$ in the small-scale region, while the model M2 tends to $\ell_{\parallel}\propto\ell_{\perp}$ and M3 is in between. This indicates that a strong mean magnetic field and guiding field will significantly promote the development of reconnection-driven turbulence. As for the late snapshot (lower panel), model M1 is similar to M4, while model M2 to M3. In the large-scale region, the anisotropic scaling of the models M2 and M3 are slightly smaller than that of M1 and M4, which is consistent with the results shown in Figure 5. We found that when self-generated turbulence reaches a statistically steady state (terminated at $t=70~{}000$), four models all are presenting the GS95 scale-dependent anisotropy in the inertial range of turbulence cascade.

In this section, we explore the evolution of CRs over time in reconnection-driven turbulence, including the gyroradius, momentum diffusion coefficient, and spectral energy distribution. The gyroradius of particles evolving with time is shown in Figure 7 containing the average gyroradius of all particles (upper panel) for four models and four accelerated particles randomly selected (lower panel) in model M2. We can see that for all four models the gyroradius increase gradually with time, and the difference is that the growth rate of models M1 and M4 is higher than that of the other two models. This should be due to the effect that both a high-$\sigma$ parameter for M1 and a non-zero guide field for M4 can improve the reconnection rate to accelerate particles more efficiently. Moreover, the growth rate of the gyroradius of model M1 (high-$\sigma$ setup) is significantly higher than that of M4. This implies that the guide field is not the only factor that improves the reconnection rate. For the phenomenon of the lowest growth rate of model M2, it can attribute to the addition of low initial noise, which makes turbulent reconnection slightly slower.

The lower panel of Figure 7 plots the time evolution of the gyroradius of selected four particles: two of them inside the initial current sheet (red and purple curves), one around (blue curve), and away from (green curve) the current sheet. For particles within the current sheet, sufficient energy is obtained in the early stages, and its gyroradius is increased by about 3-4 orders of magnitude, while the acceleration time is delayed for the particle around and away from the current sheet. It is not difficult to understand that as magnetic reconnection develops, the current sheet becomes wider and wider, and more particles are involved in acceleration.

Figure 8 shows the evolution of the momentum diffusion coefficient over time for the four models. From this figure, we can see a rapid growth of $D_{\rm pp}$ in the early stage, followed by a quasi-steady state as a plateau (at $t\sim 10~{}000$ for the model M1, while other models are a little bit later $t\sim 20~{}000$). In general, the plateau characteristics in $D_{\rm pp}$ are usually regarded as a typical second-order Fermi process (Pezzi et al. 2022; Gao & Zhang 2023). However, the current turbulent reconnection situation is more complicated because the acceleration time of individual particles is inconsistent. The calculation we provide here is statistical averaging of all particles, just to observe the diffusion behavior of particles in the reconnection region. In particular, the momentum diffusion in M1 exhibits differences from other models. This suggests that the main factor influencing momentum diffusion is the strength of the magnetic field, and the second factor may be linked to small initial velocity perturbations, the width of the current sheet, and the additional guide field.

Furthermore, Figure 9 shows spectral energy distributions of one million particles injected for the four models at 36 snapshots from $t=0$ to $70~{}000$ with an interval of $2~{}000$. From this figure, we can find the following main characteristics. First of all, as time evolves, the particle spectra from these models are broadened and moved to higher energy from the initial thermal distributions to the non-thermal ones. Their distributions at the final snapshot almost show a hard spectrum with the power-law relationship of $dN/{{\rm d}E_{\rm k}}\propto E_{\rm k}^{-1}$, which is consistent with the results of PIC simulations (Guo et al., 2014, 2015). Secondly, the spectral energy distributions enter the high-energy cut-off, because only a few fractions of the particles can be accelerated to such high energy. In addition, an important difference between the four models is that the maximum acceleration energy of model M1 (upper-left panel) reaches above $10^{4}$, while the other three models are around $10^{3.5}$. This indicates that the strong mean magnetic field promotes the reconnection acceleration of particles, allowing the accelerated particles to reach a higher energy upper limit.

To verify the reliability of the numerical results, we run two high-resolution models (one is $1024\times 512\times 1024$, and another is $512^{3}$) and compare the numerical results, including the power spectra (left panel) and anisotropy scalings of velocity (middle), as well as particles energy spectra (right) in Figure 10. Except for the resolution, other parameters of these models are kept consistent with model M3. Because of the limited computing resources, we only simulate it up to $t=30~{}000$, at which we see that the power spectra of velocities all show a Kolmogorov-like spectra $E(k)\propto k^{-5/3}$, and the spectral energy distributions of particles present as the hard spectrum with an exponent of $-1$. As shown in the middle panel, the phenomena that the anisotropic scale does not satisfy $\ell_{\parallel}\propto l_{\perp}^{2/3}$ but more conforms to $\ell_{\parallel}\propto l_{\perp}$ is due to the insufficient run time. When the evolution time is sufficiently long, it should show more significant anisotropy. Nevertheless, the numerical resolution has little effect on our results, and the simulation settings in this paper are sufficient for our moderate aim.