Gyrokinetic and kinetic particle-in-cell simulations of guide-field reconnection. I. Macroscopic effects of the electron flows

With our independent set of codes, the PIC code ACRONYM and the GK code GENE, we obtained similar values for the reconnection rates (see Fig. 1) as the ones reported by the previous comparison work, Ref. TenBarge et al., 2014. First, the normalized reconnection rates have the same value for each high/low $\beta$ case for different guide fields (although they decrease with $b_{g}$ when measured without normalization, consequence of keeping constant the total plasma $\beta$). The constancy of normalized reconnection rates in the strong guide field limit agrees with some previous studies, even when the total plasma $\beta$ changesLiu et al. (2014) (and correspondingly it was expected to have different reconnection regimes). Note that is not valid for Harris CS in the very low guide field regime, where it is known that the normalized reconnection rates are reduced with increasing $b_{g}$ (see, e.g., Refs. Pritchett, 2001; Pritchett and Coroniti, 2004; Ricci et al., 2004), as a result of the enhanced incompressibility of the plasma.

Second, we also found smaller reconnection rates in the high beta case compared to the low beta one, in agreement with previous studies of gyrokinetic simulations of magnetic reconnection.Numata and Loureiro (2014, 2015) However, we can notice in Fig. 1 that the reconnection values for both low and high plasma $\beta$ cases are still a substantial fraction of $(d\Psi/dt)/\dot{\Psi}_{N}\sim 0.1$. Therefore, fast magnetic reconnection takes place in both cases, even though the low beta runs are in a regime where there should be not dispersive waves that may explain these rates, according to the two fluid analysis in Ref. Rogers et al., 2001 (see also Ref. Ricci et al., 2004). This contradiction was already addressed in Ref. TenBarge et al., 2014 and also other recent works.Liu et al. (2014); Stanier et al. (2015); Cassak et al. (2015) Although there is no a full definitive answer to this controversy yet, all these evidences seem to indicate that fast dispersive waves, such as whistler and kinetic Alfvén waves, seem not to play the essential role to explain fast magnetic reconnection in collisionless plasmas. Since the purpose of this paper is different, we are not going to discuss further that issue in this work.

It is worth to mention that some details in the evolution of the system are not exactly the same as in Ref. TenBarge et al., 2014. For instance, in the low beta case (Fig. 1(a)), although the peak in the reconnection rate is reached more or less at the same time, we observed that the CS is more prone to the development of secondary islands than reported there. They start to appear just after the reconnection peak ($t\gtrsim 40\tau_{A}$) instead of $t\gtrsim 75\tau_{A}$. Probably, this is due to the different initial noise level in the PIC runs, since it is random and dependent on the specifics of the algorithms implementation. Therefore, the locations and timing of the secondary magnetic islands are expected to be different between the results given by the VPIC (as used in Ref. TenBarge et al., 2014) and ACRONYM codes (and also with the GK code GENE).

All the results and plots to be shown from now on are based in the low beta case. The slightly different conclusions for the high beta case will be briefly analyzed only in Sec. VI.

As indicated in Ref. TenBarge et al., 2014, based in the two fluid analysis of Ref. Rogers, Denton, and Drake, 2003, the thermal pressure and magnetic field fluctuations should display antisymmetric (odd-parity) structures in the separatrices in the strong guide field regime $b_{g}\gg 1$, assuming small enough fluctuations and asymptotic plasma beta $\beta_{y}=2\mu_{0}n_{0}k_{B}T_{i}/B_{\infty y}^{2}\gtrsim 1$ (although a later studyHosseinpour and Mohammadi (2013) showed that this assumption is not necessary for the appearance of an odd-parity magnetic field structure). In addition, both quantities should scale inversely with the guide field $b_{g}$ in the following way (see Eqs. (17) and (19) of Ref. Rogers, Denton, and Drake, 2003):

where $l_{x}$ is the typical length scale of variation of these quantities across the CS, away enough from the X points. The fluctuations $\delta$ are calculated by subtracting the (initial) equilibrium quantities $B_{z}(t=0)\approx B_{g}$ (for $b_{g}\gg 1$) and $P_{th}(t=0):=n_{e}k_{B}T_{e,\perp}+n_{i}k_{B}T_{i,\perp}$ (for brevity, we omit the subscript $\perp$ in the thermal pressure). Note that in deriving the previous expressions, the pressure equilibrium condition in the limit of strong guide field has been usedRogers, Denton, and Drake (2003) (also obtained in the framework of reduced MHD), formally equivalent (by assuming a scalar pressure) to the perpendicular force balance of the more general GK equations as given in Eq. (1) (satisfied exactly to order $\epsilon$). This can be written in the following way by combining the previous expressions for $\delta P_{th}$ and $\delta B_{z}$:

Note that Eqs. (5) and (6) predict that the fluctuations $\delta P_{th}$ and $\delta B_{z}$ are proportional to $1/b_{g}$ providing $b_{g}\gg 1$ and $l_{x}/d_{i}$ constant for different guide fields, which is valid recurring to the estimate of order of magnitude $l_{x}/d_{i}\sim\sqrt{\beta_{i}}$ given in Ref. Rogers, Denton, and Drake, 2003. The last assumption requires additional evidence in order to be applied to our simulations, and that is why we proved this claim via the method explained in detail in the Appendix B. Using the suitable normalizations explained in Appendix A, we show the estimates for $\delta B_{z}$ and $\delta P_{th}$ in Fig. 2 and Fig. 3, respectively, for a time shortly after the reconnection peak. These and all the other figures from the PIC runs shown in this paper, unless stated otherwise, have been averaged over $t=0.5\tau_{A}$ to reduce the effects of the numerical noise. In addition, the color scheme is scaled between $\pm$ the mean plus 3.5 standard deviations of the plotted quantity, a representative maximum as explained in the discussion of Fig. 5.

From Figs. 2 and  3 we confirm the result already found in Ref. TenBarge et al., 2014 and predicted by the previously sketched two fluid model:Rogers, Denton, and Drake (2003) the convergence of the PIC results towards the GK ones in the limit of strong guide field, in both odd symmetry and scaling with the guide field. As discussed by these authors, this is valid only in the region close to the separatrices. However, we can immediately notice a fact not discussed in the previous comparison work: not only the symmetry between both separatrices is broken in the low guide field regime, but also the appearance of an additional magnetic field in the secondary magnetic islands not visible in GK or in the limit of strong PIC guide field. This is the main difference and new result to be the focus of the remainder of this work. For brevity, other differences between PIC/GK will be addressed in a follow-up paper.

The core magnetic field in the secondary magnetic islands and $y$ boundaries is unbalanced, in the sense that the thermal pressure does not decrease at these locations (compare Figs. 2 with 3). This is a violation of the pressure equilibrium condition in the strong guide field limit as given in Eq. (1) or Eq. (7). Since it is satisfied exactly to order $\epsilon$ in GK, it is built-in in the equations that the GENE code solves, while PIC codes allow large deviations from it for finite guide fields, since they solve the full collisionless Vlasov-Maxwell system of equations. Note, in particular, that Eq. (7) represents a straight line with slope $-\beta_{i}$ in the plane $\delta B_{z}$ - $\delta P_{th}$. Therefore, a convenient way to display this relation and deviations is by means of a 2D (frequency) histogram of these fluctuating quantities, with results shown in Fig. 4(middle row). These plots were generated by selecting the interesting region close enough to the center of the CS with $J_{z}$ above $10\%$ of its initial value, indicated in Fig. 4(top row). In order to show the locations of deviations in the pressure equilibrium condition, we also plot in Fig. 4(bottom row) the corresponding fluctuations in the total pressure

with $P_{mag}=B^{2}/(2\mu_{0})$.

Fig. 4(middle row) shows that most of the locations in the CS are along the line $-\beta_{i}$ for the case PIC $b_{g}=20$, with an increasing spread in the lower $b_{g}$ regime. The corresponding GK results follow very accurately that line. Note a distinctive, pressure unbalanced, “bump” in the region $\delta P_{th}\sim 0$ with $\delta B_{z}\gtrsim 0$, being very noticeably for the case PIC $b_{g}=5$. Fig. 4(bottom row) shows that these regions are mostly in the secondary magnetic islands and in the $y$ boundaries. This excess of total pressure generates a net force towards the exterior of the magnetic islands, leading to an expansion in reconnection time scales ($\sim\tau_{A}$). $\delta P_{\rm total}$ increases going to the lower guide field regime, even though this quantity has been linearly scaled to $b_{g}$. That is to be expected, since the high fluctuation level in these cases breaks down the small $\epsilon$ assumption in which Eqs. (5),(6) or (7) (and also Eq. (1)) are based.

Note that the histograms for low $b_{g}$ show a strongly asymmetric distribution in comparison with GK or PIC high $b_{g}$ regime (see Fig. 4(middle row)): there are more points located in the right bottom quadrant ($\delta P_{th}>0$ and $\delta B_{z}<0$) than in the left upper quadrant ($\delta P_{th}<0$ and $\delta B_{z}>0$). This is an indication of (and an efficient way to quantify) the asymmetry in the separatrices for the PIC low guide field regime (see the respective contour plots Figs. 2 and 3). We will explain the reasons in Sec. V.1.

We can summarize these findings by saying that, although the magnetic field fluctuations $\delta B_{z}$ predicted by the GK simulations can be comparable to the PIC ones in the low guide field regime ($b_{g}=5$ and $10$) close to the separatrices, this is not true inside of the secondary magnetic islands or the periodic $y$ boundaries, due to the deviation in the pressure equilibrium Eq. (1).

The purpose of this subsection is to analyze the dynamical evolution of the thermal and magnetic fluctuations, complementing the work by TenBarge et al. (2014). This allow us to determine the physical mechanisms producing the deviations in the pressure equilibrium condition, as well as the asymmetric separatrices in the PIC low guide field regime. We quantify this by means of tracking the “maximum” of $\delta P_{th}$ and $\delta B_{z}$ in the entire region shown in Figs. 2 and 3. For the PIC simulations, in order to decrease the effects of the numerical noise, this “maximum” is chosen as equal to the mean plus 3.5 standard deviations. The absolute maximum is not a good choice since is very prone to outlier values. In addition, the initial value of the respective fluctuating quantity is subtracted (note that it is not only the initial value $P_{th,0}$ and $B_{z,0}$), because: (1) a zero offset is required by GK since the fluctuating quantities are initially zero and (2) it mostly measures the numerical PIC noise, being enhanced for higher $b_{g}$. The results are shown in Fig. 5.

First of all, Fig. 5(a) shows that the maximum of $\delta P_{th}$ is reached around $t\sim 50\tau_{A}$ for all the cases, after the reconnection peak time $t\sim 40\tau_{A}$ when also the secondary magnetic islands start to form. $\delta B_{z}$ reaches maximum values even later. This in an indication of an additional mechanism generating $\delta B_{z}$, different from the one produced by the Hall currents due to the reconnection process itself (close to the X point, in the separatrices), and deeply in the non-linear phase. This is also the justification for showing the contour plots at a time $t=50\tau_{A}$ in most of the figures of this paper.

By means of the results shown in Figs. 5(a) and  5(b), we also confirm (complementing the work of Ref. TenBarge et al., 2014) the convergence of the time evolution of both PIC $\Gamma\delta P_{th}$ and $\Gamma\delta B_{z}$ toward the GK curves in the strong guide field limit $b_{g}\gtrsim 20$. As can be expected, the GK values follow a similar trend as the strongest PIC guide field ($b_{g}=50$). However, there are important deviations in some curves, noticeable earlier for smaller $b_{g}$. For example, the PIC run $b_{g}=5$ shows deviations from the GK curve of $\Gamma\delta B_{z}$ already in $t\gtrsim 20\tau_{A}$, while $b_{g}=30$ only after $t\gtrsim 45\tau_{A}$. Note that the deviations in $\Gamma\delta P_{th}$ are smaller than those of $\Gamma\delta B_{z}$. The (interesting) physical reasons of the differences for $\Gamma\delta P_{th}$ will be addressed in a follow-up paper. Here we explain and analyze with more detail the differences in $\Gamma\delta B_{z}$ shown in Fig. 5(b).

Although useful, Fig. 5 cannot provide us with information about the relation of these maxima with the pressure equilibrium condition Eq. (7). For this purpose, in Fig. 6 we show the 2D histograms relating these fluctuations in the plane $\delta P_{th}-\delta B_{z}$ for three characteristic times during the evolution of the case PIC $b_{g}=5$.

The asymmetry along the pressure equilibrium line located in the right bottom part ($\delta P_{th}>0$ with $\delta B_{z}<0$) starts well before the reconnection peak time, indicating a process that is active since the start. Because the points tracked in Fig. 6(a) (for $t=20\tau_{A}$) are mostly located in the separatrices, we will see in Sec. V that it is associated with the initial shear flow. On the other hand, the “bump” $\delta P_{th}\sim 0$ with $\delta B_{z}\gtrsim 0$ (indicating violation of the pressure equilibrium condition and core magnetic field) starts to develop just before of the reconnection peak time ($t\sim 30\tau_{A}$, Fig. 6(b)). These points are located inside of the secondary magnetic islands or the $y$ boundaries. Therefore, and different from the mechanism leading to the asymmetric separatrices, this is an indication of being caused via a process that needs to be build up during the evolution of reconnection. Later, in Fig. 6(c) for $t\sim 40\tau_{A}$, the points spread even further in the “bump region”, indicating the shift of the largest deviations of the pressure equilibrium condition from the separatrices to the “bump”. In Sec. V we will see that the mechanism is a combined effect of both shear flow and magnetic islands generated via magnetic reconnection.

We can summarize this section by saying that PIC and GK simulations predict similar evolution for both magnetic and thermal pressures especially during the linear phase of reconnection. The formation of secondary magnetic islands breaks this similarity, producing deviations especially in the magnetic fluctuations $\delta B_{z}$, that reach maxima for times much later than the reconnection peak time.

The core magnetic field seen in the PIC low guide field runs (see, e.g., Figs. 3(a) and 3(b)) is due to the development of strong in-plane currents growing from a shear flow present in the PIC initialization. Indeed, in addition to the component $J_{z}$ sustaining the CS (associated with $B_{y}(x)$), the force free initialization with $B_{z}(x)$ implies, via Ampère’s law, the presence of an in-plane component $J_{y}$, chosen to be carried only by the electrons, and given by

This represents a counterstreaming (shear) flow of electrons along the center of each CS (see Fig. 7(b)), with maximum value $V_{e,y}=-J_{e,y}/(en_{0})\propto B_{\infty y}$. Thus, due to the normalization used (see Appendix A), its magnitude will decrease as $1/b_{g}$ (see Fig. 7(a)), while the associated kinetic energy of the shear flow has an even stronger dependence ($\propto 1/b_{g}^{4}$). Therefore, the shear flow strength in the PIC cases converges very quickly to the GK initialization, where it is completely absent. The zero initial shear flow in this plasma model is because of the perpendicular GK Ampère’s law, $(\nabla\times\delta\vec{B})_{\perp}=\nabla_{\perp}\delta B_{z}=\mu_{0}\delta\vec{J}_{\perp}$ (leading to the pressure equilibrium condition Eq. (1)). Indeed, when $J_{\perp}$ is calculated from the force free particle distribution function, it turns out to be of second order in $\epsilon$. If taken into account, this would imply a $\delta B_{z}$ of order $\epsilon^{2}$, which is ruled out by construction from the GK equations. Therefore, any shear flow (in-plane $\delta J_{\perp}$) does not enter into the GK equations since they are of order $\epsilon^{2}$.

The initial shear flow is responsible for the asymmetric separatrices (see discussion of Figs. 4 and 6), especially regarding the preferential $\delta P_{th}>0$ over $\delta P_{th}<0$. This behavior was already pointed out by Cassak (2011), where it was attributed to the dynamic pressure of the shear flow, piling up electrons preferentially in one pair of the separatrices over the other one if strong enough, contributing to the increase of density, temperature and thermal pressure (see Fig. 2). An extended discussion about this topic will be given in a follow-up paper.

The effects of the shear flow need to be carefully considered when normalized to $V_{A}$ (dependent on $b_{g}$) or $v_{th,i}$ (independent on $b_{g}$). The ratio between these two speeds is shown in Fig. 7(a), decreasing with $b_{g}$ according to the normalization used (see Appendix A). This speed ratio plays an essential role in the GK model, because of the perpendicular drift approximation. This is the assumption that the perpendicular bulk velocities $\vec{V}_{\perp}$ are caused only by the $\vec{E}\times\vec{B}$, diamagnetic $\nabla P_{e}$ and other similar drifts. The perpendicular approximation breaks down when the in-plane speeds are close to the ion thermal speeds, such as the typically obtained in magnetic reconnection outflows when $V_{A}\sim v_{th,i}$. It can also be violated by the presence of waves with high enough frequency (possibly only captured by the fully kinetic plasma description), leading to a cross-field diffusion of particles across the magnetic field.Chen, Le, and Fok (2015) Therefore, deviations from the real physical behavior of a Vlasov plasma modeled via PIC simulations are expected in the GK approach when the ratio $V_{A}/v_{th,i}\sim 1$. This is precisely the situation when comparing the GK to PIC simulations with $b_{g}=5$ or $b_{g}=10$ (the latter is the critical guide field for which $V_{A}/v_{th,i}\approx 1$).

As a result of the initial shear flow, PIC runs with sufficiently low guide field build up a net vortical current inside of the secondary magnetic islands and the periodic $y$ boundaries, as can be seen in Fig. 8. The formation of secondary magnetic islands wraps up magnetic field lines around them, deflecting the electron shear flow in the same direction. Therefore, a net out-of-plane magnetic field is generated (see Fig. 3). Note that the direction of the curl of $\vec{J}$ coincides with the direction of the out-of-plane magnetic field ($+z$ direction points into the page). As expected, there is a working dynamo process in these places, i.e.: $\vec{J}\cdot\vec{E}<0$, involving a transfer of energy from the bulk electron motion to the magnetic field (plots not shown here). High guide field PIC or GK runs do not show this effect. (Negative values of $\vec{J}\cdot\vec{E}$ are seen only near the outflows, due to the bulk motion of the plasma.) A comparison of this process to the dissipation $\vec{J}\cdot\vec{E}>0$ close to the X points will be given in a follow-up paper.

Two other consequences of the shear flow in reconnection found in previous works are worth to mention in this point. First, the tilting of magnetic islands and generation of concentric vortical flows inside of them as result of sub-Alfvénic flows (in agreement with our parameter regime) has been seen in both 2D MHD and Hall-MHD simulations (see Ref. Shi et al., 2005 and references therein). Second, the increasing symmetry of the core magnetic field for low PIC $b_{g}$ can also be explained due to the shear flow, as investigated via a two fluid analysis of tearing mode with shear flow of Ref. Hosseinpour and Mohammadi, 2013. But there is also opposite evidence to this claim: Ref. Karimabadi et al., 1999 found, by means of hybrid simulations of Harris sheets, that the “S” shape of $\delta B_{z}$ in the secondary magnetic islands is only consequence of guide field reconnection (nothing to do with shear flow). This might be due to the small guide field regime analyzed: $b_{g}<1$, not being correct to be applied for our case.

The current inside of the magnetic islands, significant only in the low guide field regime, is produced due to the Hall effect: a decoupling of motion between electrons and ions, as shown in Fig. 9. For the lowest PIC guide field run $b_{g}=5$, the ions follow the outflow due to reconnection from the X point (Fig. 9(b1)), but the electrons keep their initial shear flow and are only weakly deflected (Fig. 9(a1)) by that reconnection outflow, following a vortical flow pattern inside of the magnetic islands. This characteristic (antisymmetric) flow pattern is barely visible for a higher guide field of $b_{g}=20$ (see Figs. 9(a2)-(b2)) and totally absent for the GK run (see Figs. 9(a3)-(b3)), where the internal flow has a symmetric structure following the features of the reconnection outflow.

It is important to mention that because the generation of magnetic field is due to a Hall effect, their effects will be stronger when the CS is on the order of/thinner than $\rho_{s}=\sqrt{k_{B}(T_{e}+T_{i})/m_{i}}/\Omega_{ci}=\rho_{i}=0.1d_{i}$, the sound Larmor radius (for $T_{i}=T_{e}$), which applies very well to our case ($L=2\rho_{i}$). Therefore, in thicker CS these effects should be significantly reduced and an agreement between PIC and GK simulations of magnetic reconnection will be more easily reached. Note that the halfwidth is also on the order of electron skin depth, since $\rho_{i}=d_{e}$, and so the electron inertial effects are important:Nakamura, Fujimoto, and Otto (2008) electrons are also unmagnetized (not fulfilling the frozen-in condition) for distances $l\lesssim L/2$.

In Fig. 9, we can estimate that the speed of the electron outflows from the main X point is a few times the asymptotic Alfvén speed $V_{e,y}\sim 2.2V_{A}$. It varies weakly with the guide field in the PIC runs. On the other hand, the ion outflow speeds only reach sub-Alfvénic values $V_{i,y}\sim 0.8V_{A}$, and are much more constant among different PIC guide fields and the GK runs. These values agree, to a first order, with a two fluid theory of magnetic reconnection by Shay et al. (2001): the outflow speeds from the X point should be of the order of the in-plane Alfvén speed $V_{A}$ for ions and in-plane electron Alfvén speed $V_{Ae}=\sqrt{m_{i}/m_{e}}V_{A}$ for electrons. But the initial shear flow is of course strongly dependent on the guide field (see Fig. 7). The combination of both effects determines the critical PIC guide field for which the shear flow can generate the current that builds up the core magnetic field ($b_{g}\lesssim 20$).

Since in Sec. V.2 we showed that the electron/ion outflows do not depend strongly on $b_{g}$ when normalized to $V_{A}$, the in-plane current $\vec{J}_{\perp}\propto\vec{V}_{i,\perp}-\vec{V}_{e,\perp}$ will display similar values among different PIC $b_{g}$ when the same normalization is used, as can be seen in Fig. 8. But according to the Ampère’s law, the generation of $\delta B_{z}$ depends on the unnormalized $J_{\perp}$, changing with $b_{g}$ according to Eq. (16). It can be estimated by approximating the curl $\nabla\times\vec{B}$ by the gradient scale length $1/\Delta L$:

Since $B_{g}$ practically does not change for different PIC guide fields $b_{g}$ (result of choosing an invariant $B_{T}$, recall Sec. II), by defining the constant $\Lambda=\mu_{0}\rho_{i}/B_{g}$ we infer that (in order of magnitude), $\delta B_{z}/B_{g}\sim\Lambda J_{\perp}$ when $\Delta L\sim\rho_{i}$, not dependent on the guide field. Note that this estimate relies on the fact that $\rho_{i}$ is approximately constant during the evolution, which in turn depends on the fact that $T_{i}$ does not have to change too much. We checked that the ion temperature does not increase more than $30\%$ of its initial value throughout the evolution of the system, and therefore we can safely take $\rho_{i}$ and so $\Lambda$ constant for our purposes. This is equivalent to a proportionality between $\delta B_{z}$ and $J_{\perp}$, a relation always satisfied in a force free equilibrium. This is a valid first order approximation since the magnetic islands grow at ion time scales ($\tau_{A}$), and therefore, at electron time-scales that relation can be satisfied quasi-adiabatically. On the other hand, in our runs $\Delta L$ can be estimated as the size across the $x$ direction of either the secondary magnetic islands close to the $X$ point ($\sim 10\rho_{i}$), or the magnetic island at the $y$ boundaries ($\sim 18\rho_{i}$). Thus, the time history of the maximum value of both components of $\Lambda J_{\perp}$ is shown in Fig. 10.

The presence of the initial $J_{y}$ due to the shear flow can be seen in Fig. 10(b), increasingly important for lower PIC guide fields. But the in-plane component $J_{y}$, and to a lesser extent $J_{x}$, start to grow just after the formation of secondary magnetic islands. There is a good agreement (in order of magnitude) of the maximum values that $\Lambda J_{y}$ reach with the corresponding maximum values of $\delta B_{z}$ (see Fig. 5(b)) for different guide fields, justifying empirically our estimation Eq. (10). Note that the factor $\Gamma$ should be removed in the latter for a proper comparison (e.g., for $b_{g}=5$, the maximum is $\delta B_{z}/B_{g}\sim 0.034$). It is also important to remark that the maximum values of $J_{y}$ are reached in the boundaries of the locations where $\delta B_{z}$ peaks: around the secondary magnetic islands and the $y$ boundaries (neglecting the region in the separatrices with almost zero curl). The most important conclusion that can be inferred from Fig. 10 is that the maximum values of $J_{y}$, and so $\delta B_{z}$, become negligible in the PIC higher guide field limit when measured in absolute units. Therefore, magnetic field generation is only effective in the PIC low guide field regime.

In principle, another possibility for the core-magnetic field generation that it is interesting to analyze is due to a pinch effect and the associated magnetic flux compression inside of the magnetic islands (via conservation of magnetic flux/frozen-in condition). We checked that the divergence $\nabla\cdot\vec{v}_{\perp}$ for both electrons and ions has a complex spatial structure and mostly odd symmetry inside and around the islands, not correlating well with the properties of $\delta B_{z}$ (plots not shown here). Moreover, we checked that the absolute value of $\nabla\cdot\vec{v}_{\perp}$ (normalized to the natural time scale $\tau_{A}^{-1}$), is reduced for higher guide fields, implying a more incompressible plasma. This is consequence of the fact that the lowest order perpendicular drifts are incompressible, a condition satisfied better in this regime. However, $\nabla\cdot\vec{v}_{\perp}$ decays much faster with the guide field as the core-field: already for $b_{g}=20$, this quantity is comparable with the noise level but, on the other hand, there is still significant $\delta B_{z}$ in the islands for this guide field case. Therefore, a pinch effect $\nabla\cdot\vec{v}_{\perp}<0$ is not likely to be responsible for the generation of $\delta B_{z}$ in the magnetic islands.

All these are indications that the generation of the core magnetic field $\delta B_{z}$ is due to a combination of the effects of the shear flow and the formation of magnetic islands, becoming increasingly important for lower PIC guide fields, and that is why it should not be taken into account when comparing PIC with GK simulations of magnetic reconnection.

We already mentioned that in the PIC low guide regime, besides of the secondary magnetic islands, a core magnetic field is also generated in the $y$ boundaries. This is due to a combination of two mechanisms. The main one is because of the compression by the colliding outflows generated from reconnection in the main X point and the periodic boundary conditions (equivalent to a configuration with multiple X points). This effect is stronger for PIC low guide field regime since the electron outflows are faster (in absolute units, since they have same values in units of $V_{A}$). Note that this process could be avoided by choosing longer boxes.Karimabadi et al. (1999) The second mechanism contributing to the core magnetic field generation is the same as for the secondary magnetic islands, since the O point of the primary magnetic island is located precisely at the $y$ boundaries.

As we can see in Fig. 1, the reconnection rates are smaller for the lower guide field cases $b_{g}=5$ and $10$ compared to the PIC high guide field regime or the GK runs. This can be understood in terms of three different reasons. First, the outflows from the X points should be less efficient (slower) because of the additional magnetic pressure due to $\delta B_{z}$ in the secondary magnetic islands (and $y$ boundaries), after being normalized to the Alfvén speed. This process tends to inhibit the CS thinning. Note that this additional out-of-plane guide field has higher relative importance with respect to the asymptotic magnetic field precisely in cases with low PIC $b_{g}$ (see, e.g., Fig. 5). Second, as indicated in Ref. Cassak, 2011 with the Hall-MHD model without guide field, the magnetic tension in the reconnected magnetic field lines should be released by the shear flow, reducing the outflow speed produced by them and thus decreasing reconnection rates by a factor of $\left(1-V_{e,y}^{2}/V_{A}^{2}\right)$. Similarly, numerical solutions of a kinetic dispersion relation and 2D PIC simulations for thin CSRoytershteyn and Daughton (2008) demonstrated a reduction in the growth rate of the collisionless tearing mode for increasing shear flows, the instability associated with the onset of magnetic reconnection. And third, part of the available magnetic energy for reconnection that should be converted into particle energy is given back when the Hall currents are formed in the secondary magnetic islands. Overall, these reasons suggest that reconnection rates might be reduced in the PIC low guide field regime, if the total plasma $\beta$ is kept constant.