{'fraction_non_alphanumeric': 0.02813497885561405, 'fraction_numerical': 0.010183826702338828, 'mean_word_length': 4.3995806150978565, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Two-way coupled MHD-PIC simulations of magnetic reconnection in magnetic island coalescence

We present simulations of magnetic reconnection with a newly developed coupled MHD-PIC code. In this work a global magnetohydrodynamic (MHD) simulation receives kinetic feedback within an embedded region that is modeled by a kinetic particle-in-cell (PIC) code. The PIC code receives initial and boundary conditions from the MHD simulation, while the MHD solution is updated with the PIC state. We briefly describe this coupling mechanism. This method is suitable for simulating magnetic reconnection problems, as we show with the example of reconnection in the coalescence of magnetic islands. We compare the MHD, Hall-MHD, fully PIC and coupled MHD-PIC simulations of the magnetic island coalescence. We find that the kinetic simulations are very different from the MHD and Hall-MHD results, while the coupled MHD-PIC simulations can remedy this discrepancy while saving computing time. The diffusion region is well resolved in the kinetic simulations, which is also captured by the coupled MHD-PIC model. The coupled simulation also reproduces the kinetic Hall magnetic fields correctly. We calculate the reconnection rates and find differences between the MHD and kinetic results. We find that the coupled MHD-PIC code can reasonably reproduce the kinetic reconnection rate when a larger PIC feedback region is used, while still saving significant computing time.


Introduction
Simulating plasmas is a very challenging task because plasma processes are multi-scale as well as multi-physics in nature. For example, turbulence can be driven at large scales by fluid processes, like solar flares or some magnetohydrodynamic (MHD) instabilities. The nonlinear interactions at the high Reynolds and Lundquist numbers in space and astrophysical plasmas give rise to energy transfer to multiple scales, going right down to the ion and electron scales, which are several orders of magnitude smaller [1,2] than the driving scales. In solar corona MHD processes operate at 10 9 m while electron scales lie at 10 −1 m [3]. Similarly in Earth's magnetosphere the solar wind driver operates at 10 9 m while electron scales lie at 10 3 m [4]. Simulating this range of scales is not possible even with the computing facilities of foreseeable future. However, an even greater challenge is that while the plasma can be treated as a fluid at the largest scales, the kinetic, particulate nature of plasma becomes critical at the small scales. MHD models are fluid models which work at scales larger than the ion gyroradius or skin depth [5] and are computationally ideal for global simulations. However, they do not capture the kinetic physics that results from the particulate nature of plasma. The Maxwell-Vlasov equations are typically used to describe a plasma kinetically and many advanced particle-in-cell (PIC) techniques have been developed to simulate them. However, these are significantly more computationally demanding then MHD as they simulate a 6D phase space compared to the 3D physical space of fluid models. Moreover, they need to resolve the finer kinetic length and time scales. Hence, implicit PIC methods have been developed to allow larger time-steps and gridspacing [6]. Despite all this, it is still very difficult to simulate global space and astrophysical plasmas kinetically for physically relevant times.
To overcome some of these problems, we present simulations with a coupling of MHD and PIC in physical space. We simulate the well-known magnetic island coalescence problem with a coupled MHD-PIC code where the global simulation is done with MHD, while a domain around the reconnection region is simulated with PIC. We briefly describe the details of this two-way coupling in Sec. 2. Sec. 3 describes the problem setup of magnetic island coalescence which is a useful problem to study flux-rope merging in tokamaks as well as solar corona [7,8]. We also do a comparison of the results produced by MHD, Hall-MHD, coupled MHD-PIC, and fully PIC simulations. We find that the coupled simulation saves computing time compared to a full PIC simulation while capturing the relevant kinetic physics effects which are not reproduced in MHD and Hall-MHD. In Sec. 4 we show the reconnection rates obtained in these simulations and how we can get the correct reconnection rates from the coupled MHD-PIC approach. The results are summarized and discussed in Sec. 5.

Two-way coupling of MHD with PIC
The coupling is implemented between the codes MPI-AMRVAC [9] and iPIC3D [10,11]. MPI-AMRVAC solves a variety of hyperbolic advection problems, including MHD, with conservative, shock-capturing finite volume and finite difference schemes. It offers versatility in terms of geometry of the grid, dimensionality, the spatial and temporal discretization, and computer platform [12,13,14]. It has the resistive MHD and Hall-MHD modules which are critical in coupling with a kinetic code. It also has a grid-block based adaptive mesh refinement scheme. In this work we have used the TVDLF spatial discretization with a CADA3 limiter function [15], along with a three-step RK time stepping scheme. The ∇ · B = 0 condition is maintained by a generalized lagrangian multiplier (GLM) scheme [16]. iPIC3D is a semi-implicit moment method for PIC simulations of plasmas. It solves a linearized version of the implicit Ampère's law to advance the electric field. Then the magnetic field is obtained from the induction equation and the particles are also moved implicitly. As a result of this semi-implicit scheme, iPIC3D can take larger grid spacing (10 × −50×) and larger time steps (5 × −10×) compared to explicit PIC codes.
This coupling is feasible because of the complimentary characteristics of the two codes. MPI-AMRVAC has an adaptive mesh which we can utilize to zoom in with a high resolution around the area where kinetic effects are expected to be important. As a result the grid spacing of the two codes can be made similar. The time-step of MPI-AMRVAC is determined by the CFL condition where Δt < Δx/v f , where v f is the characteristic speed of the fastest wave, typically the fast magnetosonic wave. On the other hand, the grid spacing of an explicit PIC code is limited by the Debye length Δx < λ D and the time step is restricted by the speed of light, Δt < Δx/c. These scales are very different from the scales of an MHD simulation model. However, with the semi-implicit approach of iPIC3D these constraints are relaxed. In iPIC3D the grid spacing is limited by the electron skin depth Δx < d e and the time step is limited by the electron thermal velocity Δt < Δx/u th,e . As a result the grid spacing of the MHD and PIC are not very different and with mesh refinement the MPI-AMRVAC grid size can be brought within the iPIC3D grid spacing. Assuming the same grid spacing in MHD and PIC, the ratio of the time steps in the two schemes is Thus by cleverly adjusting the electron beta β e , electron temperature T e , and the ion-to-electron mass ratio, we can have similar time steps in MPI-AMRVAC and iPIC3D, making their coupling feasible.

Figure 1.
Cartoon showing the coupling scheme in 2D. The blue region at the upper level represents a global MHD simulation. The green region in the middle (reproduced on the lower level) is the domain within the global simulation which is simulated with PIC. Yellow is the interface region between MHD and PIC. The MHD quantities like density (ρ), pressure (p), velocity (v) and magnetic field (B) are exchanged between the two codes for coupling.
The concept of coupling is shown in Fig. 1. This is very similar to the coupling implemented between the MHD code BATS-R-US and iPIC3D [17,18]. The global domain is simulated with MPI-AMRVAC. At some time step, a region within the global domain is selected by the user to be simulated with iPIC3D. In MHD to PIC coupling, information from MPI-AMRVAC is used to provide initial conditions to iPIC3D to start the simulation. The primitive variables of MHD namely density (ρ), pressure (p), velocity (v), and magnetic field (B) are supplied to iPIC3D. Using this the electromagnetic fields are initialized. The particles are initialized with a Maxwellian distribution using the thermal velocity derived from MHD pressure and drift velocity derived from MHD velocity and current density. After this MPI-AMRVAC is advanced by one time step and the time averaged state at the interface between MHD and PIC is supplied to iPIC3D as the boundary conditions for its time step. After iPIC3D advances its state to the next time step, moments of the distribution function are gathered to obtain the MHD primitive variables. In PIC to MHD coupling this information from iPIC3D is used to update the MPI-AMRVAC solution. The MHD solution is overwritten by a weighted average of PIC solution and the MHD solution in the PIC domain. The weight function varies smoothly across the interface region such that the updated MHD solution smoothly changes from purely MHD at the edge to purely PIC in the interior. This is the two-way coupling method. A more detailed explanation with several tests is provided in Ref. [19].
This coupling optimizes the parallelization in an efficient way. In Fig. 2 we show the parallelization strategy of this coupling. MPI-AMRVAC is an adaptive mesh code which utilizes the distribution of MPI-AMRVAC grid blocks on different processors according to the Morton curve, shown as an example in Fig. 2. There is a region of refined grid blocks in the center out of which the dashed rectangle is selected to be simulated with PIC. The grid blocks are arranged in a complex manner on different processors, shown by 4 different colors for 4 different processors. On the right we zoom into the iPIC3D region. The PIC simulation is distributed and run on the same set of processors as MPI-AMRVAC (represented by the same colors). This is done by running the iPIC3D as a library of MPI-AMRVAC. As a result of this each processor handles some part of the MHD domain as well as some part of the PIC domain. This means that each processor is continuously advancing the MHD and PIC simulations, alternately, maximizing the use of all processors. However, the physical domains of MHD and PIC simulations processed by a processor can be totally different. As can be seen, for ex., the red processor processes very different physical regions in the MHD and PIC parts. Therefore there are more communication costs in this approach. The communication is currently handled via master processors of MHD and PIC which gather/distribute information and communicate it to the other program via MPI. We now apply this method to simulate the coalescence of magnetic islands.

Numerical setup of magnetic island coalescence
We setup the problem of magnetic island coalescence given by the Fadeev equilibrium [20]. This has the magnetic vector potential given as [21,22] where = 0.4. The density is given by The asymptotic magnetic field is kept at B 0 = 0.1. The background density is taken as n b = 0.3n 0 where the ion mass is normalized to unity and n 0 is also unity. The velocity is zero and the pressure is given by to maintain the equilibrium. The initial condition is perturbed by a sinusoidal perturbation of the form given in Ref. [23] to initiate the coalescence of the two islands. Fig. 3 shows the state of the MHD solution with the coupled PIC feedback after 0.5τ A , where τ A is the Alfvén time. The black box in Fig. 3 represents the domain where PIC feedback is supplied. The PIC domain is inside −16 ≤ x ≤ 16 and −16 ≤ y ≤ 16. This coupling is done with MHD. Initially the magnetic islands start out with a separation of 2πλ = 30d i which is at the edge of the PIC box. After t = 0.5τ A they have smoothly moved inside the PIC box. Across all the 8 primitive variables of MHD we find that the coalescing islands smoothly merge across the MHD-PIC interface. In Fig. 3(e) we see that the inflow velocity v x varies continuously across the MHD-PIC interface. We use a smoothing function of width 2d i to average the MHD and PIC quantities at the interface [19].
The simulations are carried out in a box of size L x = L y = 60d i , where d i is the ion skin depth. The parameter λ which is the island width is set by the relation L x = 4πλ. The grid is taken of (600, 600) cells giving a resolution of Δx = Δy = 0.1d i . We simulate this setup with the MHD module of MPI-AMRVAC, with the Hall-MHD module, with the coupled MHD-PIC approach and with a fully iPIC3D simulation. For the MHD simulation the resistivity is taken to be η = 0.003 in code units which is required at this resolution in order to have a resolved current sheet forming between the two coalescing islands. For lower resistivity values we find that the thickness of the current sheet decreases as the resolution is increased, implying that it is not resolved. With η = 0.003 the Lundquist number is S λ = λv A /η ≈ 160. In the coupled MHD-PIC simulations, the resistivity in the MHD is kept to zero. The Hall-MHD module adds the extra Hall-term (η h /ρ)J × B in the Ohm's law equation [9]. In units where the length is normalized to a reference ion skin depth and velocity is normalized to speed of light, the Hall parameter η h is equal to the reference Alfvén velocity normalized to c which is simply unity. In the coupled MHD-PIC simulations, the resistivity (η) and the Hall parameter (η h ) in the MHD part are kept zero. For the MHD and Hall-MHD simulations the boundary conditions are symmetric for the density, pressure, tangential velocity and normal magnetic field component, while antisymmetric for the normal velocity and tangential magnetic field along the x-boundaries. It is continuous for the density and pressure, symmetric for tangential velocity and tangential magnetic field, while antisymmetric for the normal velocity and magnetic field in the y-boundaries. PIC simulations the MHD domain has boundary conditions as described above for MHD and Hall-MHD while the PIC domain receives boundary conditions from the MHD state [19]. The time step for the MHD simulations with these parameters is kept at 0.2ω −1 pi which satisfies the CFL condition, and the same time step is used in the PIC simulations as well as in the coupled MHD-PIC simulations. The time step for Hall-MHD is significantly reduced to around 0.008ω −1 pi due to the dispersive whistler waves.   Fig. 3(f)). Therefore the x axis represents the coordinate along the x axis in all figures, except in (f) where it shows the y co-ordinate.The green line is taken from the MHD simulation, the blue line is from the coupled MHD-PIC simulation, while the red line is taken from the fully PIC simulation. This is the state of the simulations after 1.0 Alfvén time.
The solutions along the X and Y axes are shown in Fig. 4 at 1 Alfvén time. This is the time around which the reconnection rate peaks. Three simulations are shown in this figure, the green line is the MHD simulation, the blue line is from the coupled MHD-PIC simulation, while the red line is from the fully PIC simulation. The magnetic field B x is close to zero along the x axis because that is where it flips sign, therefore it is noisy in all the three simulations. The vertical field B y agrees well between the coupled and PIC simulations, but is stronger and more concentrated in the MHD simulation. The out-of-plane magnetic field is zero in MHD, but non-zero in the coupled and PIC simulations, which agree with each other. The density and pressure match well across the three simulations, with a closer match between the coupled and PIC simulations. There is a spike in the pressure at the central current sheet for the MHD which is not present in the coupled and PIC simulations. is taken along the y-axis unlike the other panels, because there is no interesting behavior along the x axis. Here the PIC and coupled solution match well except at the boundaries where the PIC solution develops some numerical effects. We see that the outflow is much stronger and very different in the MHD simulation than in the coupled and PIC simulation. If we look at the 2D profile of the outflow for MHD, we find that this is due to the fact that the outflow in MHD is much more collimated while it is much more spread out in PIC (Fig. 3(f)). The out-of-plane velocity v z is again zero in the MHD description but shows a strong negative value near the central region in both coupled and PIC simulations that agree quite well with each other. Outside the PIC feedback, the coupled solution diverges from the PIC solution.
In Fig. 5 we show the out-of-plane Hall magnetic field around the diffusion region for the 4 different simulations at 1 Alfvén time. MHD does not develop a Hall field whereas Hall-MHD produces a quadrupolar structure around the central current sheet. However, this structure is different from that produced in the coupled MHD-PIC and fully PIC simulations. The coupled and PIC simulations match quite well showing a thicker region with prominent B z field. The strength of B z is also much stronger in the PIC and coupled simulations compared to Hall-MHD. The current density for the four different simulations is showed in Fig. 6. The current sheet that forms in MHD looks very different than the current sheet in PIC and coupled simulations. The Hall-MHD current sheet has a lower aspect ratio compared to MHD and is closer to a point-like geometry. We have found that to obtain a correct X-point like Hall-MHD structure, we need to go to significantly lower resistivity, higher resolution, and also larger system size. However, both MHD and Hall-MHD show sheets that are thinner than the ion skin-depth [22]. On the other hand, the coupled and PIC simulations show well-resolved current sheets. Their thickness is close to 1 ion skin depth d i and does not change with resolution. This is similar to results from simulations of turbulence as well [24,25]. Also their structure is closer to an X-point geometry.

Reconnection rates
The reconnection rate of the merging islands is an important property of the underlying physics of coalescence. To calculate the reconnection rate, the magnetic flux along the x axis is defined as ψ(x) = x 0 B(y)dx. This is calculated by integrating along the x-axis starting from the X-point which lies at the origin. The reconnection rate (E R ) is calculated as [21] where ψ X and ψ O are the fluxes at X and O points respectively. The O-point is identified as the center of the island where the flux is maximum. For comparison with Ref. [21] it is normalized to the Alfvén velocity v Am and magnetic field B m at the point between the islands where the magnetic field is maximum, which is around 0.035 in our setup.  Table. 1. The coupled MHD-PIC simulations match with the fully PIC simulations for larger embedded PIC domains.
We calculate the reconnection rates for several different simulations and they are plotted in Fig. 7 as a function of time. The 4 simulations described above are labeled as MHD 0, Hall 0, Coupled 0, and PIC 0. In addition, there are several other simulations performed with the coupled MHD-PIC code where different sizes of the embedded PIC region were used. For ex., in Coupled 0 the embedded PIC domain was within −16 ≤ x ≤ 16, −16 ≤ y ≤ 16. In Coupled -1, a slightly smaller PIC region is used, lying within −12 ≤ x ≤ 12, −12 ≤ y ≤ 12. These sizes are specified in Table. 1 for the different simulations. We find that the MHD 0 and Hall 0 simulations show similar time evolution of the reconnection rate, with the Hall simulation giving a higher peak reconnection rate than MHD. They peak later than the PIC simulation, closer to 1 Alfvén time and after t = 1.2τ A their reconnection rate drops significantly compared to PIC 0. The Coupled and PIC simulations agree in their reconnection rates up to 0.7τ A . After that the coupled simulations start to show a higher reconnection rate. We see that as the domain of PIC feedback is reduced, the reconnection rate shoots higher, overestimating the actual PIC rate. As  In Table. 1 we summarize the reconnection characteristics of the different simulations. The peak reconnection rate is reported and also the time at which it occurs τ P R . The reconnection rate is averaged from time 0.5τ A to 1.5τ A and is also reported. The peak reconnection rate of 0.438 in PIC 0 simulation matches well with the result of Ref. [21]. The MHD, Hall MHD, as well as coupled MHD-PIC simulations overestimate the peak and average reconnection rates compared to the PIC result. However, the match between coupled and PIC simulations gets better as the size of PIC feedback region increases. The result of the Coupled -1 simulation is reasonably close enough to the PIC result with an error of less than 20% in the average reconnection rate. However, the Coupled -1 simulation took only 26% of the computing time of the PIC 0 simulation, leading to significant computational savings.

Conclusions
We have simulated the magnetic island coalescence problem with different simulation methods. This is a subtle problem where neither MHD nor Hall-MHD are able to resolve the reconnection dynamics properly, compared to a fully kinetic code. The MHD and Hall-MHD simulations produce very thin current sheets at Lundquist numbers greater than 10 3 which need a very high resolution to resolve. On the other hand, the kinetic simulations show a well resolved current sheet which does not become thinner than 1 ion skin depth. The kinetic simulations show a X-point like geometry, while the MHD produces a Sweet-Parker geometry and Hall-MHD shows a lower aspect ratio current sheet. We also simulate this with a new coupled code where a PIC simulation is embedded in a smaller region within a global MHD simulation. This code coupling works reasonably well for this problem. It reproduces the correct current sheet structure and magnetic field structure as a fully PIC simulation, while taking significantly less computing time than a fully PIC simulation. The out-of-plane Hall magnetic field is actually much broader and stronger in the kinetic simulation than in Hall-MHD. For the out-of-plane velocity and the inflow velocity, the coupled method reproduces the kinetic solution well inside the PIC feedback region, but there is some divergence outside of it.
The reconnection rates are also calculated. We find that the MHD and Hall-MHD show a different reconnection rate time evolution compared to a fully PIC solution. The peak and average reconnection rates are overestimated in the MHD and Hall-MHD simulations. Coupled MHD-PIC simulations also overestimate the reconnection rate when the size of PIC feedback region is small. As the size of PIC feedback increases, we recover the PIC reconnection rate, but we still save significant computing time. This problem is a test case for flux tube merging in solar corona and tokamaks. It will be useful to apply this for larger system sizes to obtain scaling results of the reconnection process with system size. This will also help us in studying this process at large, realistic Lundquist numbers.