3.1. Reconnection rate and structure

The reconnected flux and reconnection rate of simulation A are shown in figure 1(a). The reconnected flux $\Psi$ is computed as the difference between the maximum and minimum of the out-of-plane component of vector potential $A_z$ along the upper current sheet. The reconnection rate $E_r$ is calculated as $E_r=\partial \Psi /\partial t$ and is normalized by the in-plane asymptotic magnetic field $B_0$ and the inflow Alfvén speed $V_A=B_0/\sqrt {\mu _0 \rho }$ . We find that reconnection goes through a slow ‘build-up’ phase at first and this lasts until approximately $t=80\omega _{ce}^{-1}$ . Then the reconnection rate increases steadily and reaches the first peak at $t=265\omega _{ce}^{-1}$ . Afterwards, it decreases first and later returns to increasing from $t=380\omega _{ce}^{-1}$ , which is caused by the coalescence of magnetic islands. A second peak is then reached, followed by a new decrease. It is noted that the first peak rate is larger. In our simulations, due to the relatively small size of the box in the $y$ direction, the lack of ’fresh’ plasma to process in the reconnection region prevents reaching and maintaining the steady-state (Wan & Lapenta Reference Wan and Lapenta2008).

Figure 1. Reconnected flux and reconnection rate in the three cases: (a) case A ( $B_g=5B_{x0}$ ); (b) case B ( $B_g=10B_{x0}$ ); (c) case C ( $B_g=20B_{x0}$ ). The orange dashed lines represent the time snapshots used to compare the three simulations. These three snapshots have the same reconnected flux and roughly the same island size. The reconnection rate $E_r$ is normalized by the inflow Alfvén speed $V_A$ and in-plane asymptotic magnetic field $B_0$ .

The reconnected flux and reconnection rate of simulations B and C are illustrated in figure 1(b,c), respectively. They have similar patterns as case A. The reconnection is very slow initially, and then after the ‘build-up’ phase, both the reconnected flux and the reconnection rate increase gradually. Then the rate reaches a peak value and decreases later while the reconnected flux still increases. However, cases B and C do not exhibit a second peak, unlike case A. This is because simulations B and C are not run to the point where magnetic islands coalesce. The coalescence does not affect our results since we carry out the analysis at an earlier time (as shown below). We find that the reconnection is delayed with large guide fields due to reducing electron gyroradius. The reconnection rate starts to increase from $t\approx 100\omega _{ce}^{-1}$ in case A, $t\approx 300\omega _{ce}^{-1}$ in case B and $t\approx 1000\omega _{ce}^{-1}$ in case C, respectively. This is because when the guide field increases, the ratio between the initial current sheet thickness (which is kept fixed in our simulations) and the electron gyroradius becomes larger and it takes longer for the sheet to become narrow enough and for the field lines to reconnect, as pointed out by Shay et al. (Reference Shay, Drake, Swisdak and Rogers2004). Furthermore, it is noted that the maximum reconnection rate decreases when the guide field increases (to compare the reconnection rate, as done in Ng et al. (Reference Ng, Huang, Hakim, Bhattacharjee, Stanier, Daughton, Wang and Germaschewski2015), Stanier et al. (Reference Stanier, Daughton, Simakov, Chacon, Le, Karimabadi, Ng and Bhattacharjee2017) and Guan et al. (Reference Guan, Lu, Lu, Huang and Wang2023), we use the maximum reconnection rate here since our simulations do not reach a steady state). These results suggest that large guide fields delay the onset of reconnection and decrease the maximum reconnection rate, which is analogous to previous ion-coupled reconnection works (Hesse et al. Reference Hesse, Kuznetsova and Hoshino2002; Ricci et al. Reference Ricci, Brackbill, Daughton and Lapenta2004; Muñoz & Büchner, Reference Muñoz and Büchner2016). In the following sections, we will analyse and compare the three simulations when the plasmoid instability is well-developed. The times chosen for the comparison are indicated by the orange dashed lines in figure 1, and we select those times because they correspond to the same reconnected flux. Since the reconnected flux is a measure of the islands’ size, the islands at these three times have roughly the same size.

In figure 2(a–c), we compare the out-of-plane magnetic field variations $\delta B_z=B_z(t)-B_z(t=0)$ of the three runs at the times indicated in figure 1. Here, $B_z(t)$ denotes the out-of-plane magnetic field at time $t$ (so $B_z(t=0)$ represents the initial out-of-plane magnetic field). Solid lines represent the in-plane component of magnetic field lines and we notice the presence of a well-developed reconnection region in all three runs, with magnetic islands of comparable size. We find that $B_z$ increases along the separatrix in the upper left-hand and lower right-hand quadrants (where $\delta B_z\gt 0$ ), whereas it decreases along the other separatrix in the upper right-hand and lower left-hand quadrants (where $\delta B_z\lt 0$ ). This result is the same as the quadrupolar structure in the Hall reconnection model (Zweibel & Yamada Reference Zweibel and Yamada2009; Muñoz Sepúlveda Reference Sepúlveda and Alejandro2015) and is consistent with Sharma Pyakurel et al. (Reference Sharma Pyakurel, Shay, Phan, Matthaeus, Drake, TenBarge, Haggerty, Klein and Cassak2019). It is noted that the perturbation is relatively small and would not change the sign of the total out-of-plane magnetic field. In addition, we also find the structure of $\delta B_z$ varies as the guide field increases. In case A, $\delta B_z$ is less localized, exhibiting wide spatial variations, as shown in figure 2(a). Furthermore, $\delta B_z$ is negative at the X-point and its neighbouring areas along the neutral sheet. However, in cases B and C, the structure of $\delta B_z$ along the separatrices becomes more localized and appears to be nearly antisymmetric with respect to the neutral sheet, as shown in figure 2(b,c), respectively. This localization effect of the quadrupolar structure is sensitive to the guide fields, as it becomes more and more prominent with the increase of the guide field strength, as illustrated in figure 2(a–c).

Figure 2. (a,d,g) Here $t=330 \omega _{ce}^{-1}$ in case A ( $B_g=5B_{x0}$ ); (b,e,h) $t=780 \omega _{ce}^{-1}$ in case B ( $B_g=10B_{x0}$ ); (c,f,i) $t=2440 \omega _{ce}^{-1}$ in case C ( $B_g=20B_{x0}$ ). (a–c) Structure of the out-of-plane magnetic field variations $\delta B_z$ ; (d–f) electron density $n_e$ and electron velocity streamlines; (g–i) electron velocity $V_{e,x}$ and ion velocity $V_{i,x}$ along $x$ at $y=15.0$ . In each of (a–c) and (d–f), the three panels share the same colourbar, which is located on the far right-hand side of the row. In (g–i), the red line represents the electron velocity $ V_{e,x}/V_A$ , plotted against the left-hand vertical axis. The blue line represents the ion velocity $ V_{i,x}/V_A$ , plotted against the right-hand vertical axis. Note that the velocity scales for electrons and ions are different, with $ V_{e,x}/V_A$ on the left-hand side and $ V_{i,x}/V_A$ on the right-hand side.

Asymmetric density structures have been observed around the separatrices in standard ion-scale guide field reconnection (Kleva et al. Reference Kleva, Drake and Waelbroeck1995; Pritchett & Coroniti Reference Pritchett and Coroniti2004; Pritchett Reference Pritchett2005; Birn & Priest Reference Birn and Priest2007; Lapenta et al. Reference Lapenta, Markidis, Divin, Goldman and Newman2011; Markidis et al. Reference Markidis, Lapenta, Divin, Goldman, Newman and Andersson2012). It is pointed out that electron parallel motion accelerated by the parallel electric field (along newly reconnected magnetic field lines) contributes to density cavities (low density regions) along one separatrix and density enhancements (high density regions) along the other separatrix. Similar density structures are also obtained in our simulations. Figure 2(d–f) illustrate the electron density (represented by the coloured contour plots) in the three simulations. Electron density enhancements develop along the separatrix in the lower left-hand and upper-hand right quadrants. In contrast, electron density cavities reside on the other separatrix in the upper left-hand and lower right-hand quadrants. For clarity, we hereafter refer to them as the ‘high-density separatrix’ and the ‘low-density separatrix’, respectively. The in-plane motion of electrons in all cases is shown in figure 2(d–f), where streamlines with arrows illustrate the in-plane component of the fluid velocities of electrons. The electron fluid motion exhibits patterns correlated to the out-of-plane magnetic field variations and density asymmetry structures as expected, since these structures are due to the still magnetized electrons carrying the magnetic field in the out-of-plane direction. It is remarked that the spatial extent of both quadrupolar structures and density enhancements/cavities is confined to a few electron inertial lengths in E-REC, typically less than one ion inertial length ( $d_i$ ) (along the outflow direction). This contrasts with the standard ion-coupled reconnection regime, where such structures generally span several ion inertial lengths (Pritchett & Coroniti Reference Pritchett and Coroniti2004). It is noted that the system size in the $x$ direction for all three simulations corresponds to $40\rho _e$ , $80\rho _e$ and $160\rho _e$ , respectively. These values confirm that the box size is sufficiently large, ensuring that the boundary conditions are not a significant factor in influencing the spatial structure of the reconnection region. As expected, electron and ion velocities also differ in magnitude. Figure 2(g–i) illustrate the velocities of electrons and ions (normalized by the Alfvén speed $V_A$ ) in the $x$ direction along $y=15.0$ (marked by the yellow dashed lines in figure 2 a–c) in all three cases. The electron outflow jets are observed in the results. The ion velocity in the cases A and B is very small as shown in figure 2(g,h) (the ion velocities are still nearly zero quite far away from the X-point along the current sheet). In standard ion-coupled reconnection, the ion outflow velocity increases linearly from the X-point until reaching the upstream Alfvén speed at the boundary of the ion diffusion region. As a result, the ion flow velocity at the boundary of the electron diffusion region (EDR) (where the electron outflow velocity reaches its peak) should be $(d_e/d_i)V_A=0.1V_A$ (since the mass ratio is 100 in our simulations), which could be used as a criterion to determine the extent of ion-coupling to the reconnection process. In both case A and case B, the ion velocity $V_{i,x}$ does not reach the $0.1V_A$ . The results are consistent with Sharma Pyakurel et al. (Reference Sharma Pyakurel, Shay, Phan, Matthaeus, Drake, TenBarge, Haggerty, Klein and Cassak2019). Furthermore, this also implies that the in-plane currents are predominantly carried by electrons (Kleva et al. Reference Kleva, Drake and Waelbroeck1995). In case C, $V_{i,x}$ reaches approximately $0.067V_A$ at the boundary of the EDR. This value is below the commonly cited $0.1V_A$ threshold for ion-coupled reconnection, but it is noticeably higher than in cases A and B, indicating an enhanced ion response.

The guide field has noticeable effects on density structures and particle motion patterns (Kleva et al. Reference Kleva, Drake and Waelbroeck1995). With increasing guide field strength, the regions of density enhancement and cavity become narrower, as illustrated in figure 2(d–f). Additionally, the maximum (minimum) values in the density enhancement (cavity) regions increase (decrease) with growing guide fields. The $\delta B_z$ (perturbed magnetic field in the out-of-plane direction), which either adds to or subtracts from the guide field, changes the magnetic pressure. Given the low $\beta$ , the gas pressure must adjust to approximately balance the total (magnetic plus gas) pressure. Consequently, we obtain $\delta nk_B T \approx \delta (B^2/(2\mu _0)) = B_{z0}\delta B_z/\mu _0$ . Therefore, the density variation $\delta n$ becomes more pronounced when the guide field strength is stronger, as previously discussed in Kleva et al. (Reference Kleva, Drake and Waelbroeck1995). The narrowing effect also appears in the electron flow results. The in-plane electron velocity exhibits a rotational pattern along local density cavities and enhancement regions. It turns more localized with increasing guide fields, as depicted in figure 2(d–f). The difference of $\delta B_z$ structure is correlated to the density enhancement/cavity region change. The electrons moving along separatrices undergo gyromotion around the magnetic field lines, and their gyroradius becomes smaller as the guide field increases, which leads to the stronger localization of density variation regions. Since the structure of $\delta B_z$ is caused by the in-plane current, especially the motion of electrons, its structure also narrows with $B_g$ increasing (Pritchett & Coroniti Reference Pritchett and Coroniti2004; Swisdak et al. Reference Swisdak, Drake, Shay and McIlhargey2005). We also observe that the electron velocity increases slightly with larger $B_g$ . We also plot the value of $\sqrt {V_{e,x}^2+V_{e,y}^2}$ as shown in figure 3, which confirms this result further and exclude the effect of jet rotation (here $V_{e,x}$ and $V_{e,y}$ are the x and y components of electron velocities, respectively). Starting from the electron momentum equation, the Lorentz force term has been identified as contributing to the in-plane electron outflow velocity in E-REC simulations (Guan et al. Reference Guan, Lu, Lu, Shu and Wang2024). Consequently, with an increasing guide field, its influence becomes more significant, leading to an enhanced outflow velocity. We find the density variations along the enhancement/cavity regions are relatively modest compared with the background density, a feature notably different from the dramatic variations observed in previous standard ion-coupled reconnection studies (Pritchett & Coroniti Reference Pritchett and Coroniti2004; Pritchett Reference Pritchett2005). Furthermore, we observe that the density does not show substantial variations inside the islands in the outflow. This is a remarkable difference with respect to ion-coupled reconnection, where density typically increases inside magnetic islands, as reported in Pritchett & Coroniti (Reference Pritchett and Coroniti2004) and Pritchett (Reference Pritchett2005). However, our simulation set-up, including the initial equilibrium, boundary conditions, guide field strength and simulation domain, differs significantly from those used in Pritchett’s studies. These differences could play a role in the differences we observed, and other in-depth investigations are required to pinpoint the reason for the mismatch.

Figure 3. The in-plane velocity of electrons in three cases are shown. Here $V_{ex}$ and $V_{ey}$ are the x and y components of electron velocities, respectively.

3.2. Anisotropic electron heating

Previous studies on guide field reconnection have revealed anisotropic heating observed at different locations in both experiments and simulations (Guo et al. Reference Guo, Horiuchi, Cheng, Kaminou and Ono2017; Pucci et al. Reference Pucci, Usami, Ji, Guo, Horiuchi, Okamura, Fox, Jara-Almonte, Yamada and Yoo2018; Shi et al. Reference Shi, Srivastav, Barbhuiya, Cassak, Scime and Swisdak2022a ,Reference Shi b , Reference Shi, Scime, Barbhuiya, Cassak, Adhikari, Swisdak and Stawarz2023). In this study, we investigate this problem with a broader range of intense guide fields. Figure 4(a–c) depict the electron temperature anisotropy $A_e=T_{e,\perp }/T_{e,\parallel }-1$ at the same time as figure 2 for the three simulations, with $T_{e,\perp }$ and $T_{e,\parallel }$ representing the electron temperature perpendicular and parallel to the local magnetic field, respectively. Positive values of $A_e$ indicate that perpendicular temperature exceeds parallel temperature locally (caused by perpendicular heating and/or parallel cooling), whereas negative values of $A_e$ imply that parallel temperature is larger than perpendicular temperature locally (caused by parallel heating and/or perpendicular cooling). Our analysis reveals a consistent presence of significant anisotropy along separatrices in all cases. Specifically, negative anisotropy (depicted by blue regions) is observed along the high-density separatrix and at the X-point, which means that parallel temperature is dominant in these regions. Conversely, positive anisotropy (illustrated by red regions) is situated along the low-density separatrix, indicating relatively larger perpendicular temperature in these areas. However, notable differences are observed with varying guide field strengths. As the guide field increases, the highest value of anisotropy increases, the anisotropic regions narrow, and the zones of positive anisotropy shift closer to the X-point (although $A_e$ remains negative at the X-point). This shift is caused by the reduction of the thickness of the diffusion region. The EDR becomes thinner with the decrease of the electron gyroradius stemming from the larger guide fields (Ricci et al. Reference Ricci, Brackbill, Daughton and Lapenta2004). Since the reconnection rate does not change much, the length of the reconnection region also decreases in the outflow direction. As a result, the positive anisotropy along the separatrix moves closer to the X-point. Apart from the separatrices, pronounced anisotropy also emerges in the outflow regions. In case A (figure 4 a), negative anisotropy dominates within the outflow regions. Nevertheless, negative anisotropy decreases in magnitude and becomes closer to zero at a certain location (the faint white region near $x\approx 11, y\approx 16$ shown in figure 4 b) in the outflow regions when the guide field increases. In case C, the anisotropy at this location is observed to become positive, as illustrated in figure 4(c).

Figure 4. Electron temperature anisotropy $A_e$ , normalized deviation of perpendicular and parallel temperatures: $T_{e, \perp }/T_{e,0}-1$ and $T_{e, \parallel }/T_{e,0} -1$ in the three cases. Here $T_{e,0}$ is the initial electron temperature. The coloured boxes represent the positions at which EVDFs are calculated in these three cases. The rectangle regions in (a–c) represent the areas selected for plotting scatter plots of $T_{e, \parallel }$ and $(n/B)^2$ which are shown in details below. In each row, the three panels share the same colourbar, which is located on the far right-hand side of the row.

To gain insight into the underlying mechanisms driving the observed temperature anisotropy, we examine the parallel and perpendicular electron temperature (where ‘parallel’ and ‘perpendicular’ refer to directions relative to the total local magnetic field, including the guide field), which are defined as $T_{e,\parallel } = \boldsymbol{b} \cdot \boldsymbol{P}_e \cdot \boldsymbol{b} / n_e$ and $T_{e,\perp } = (3\, T_e - T_{e,\parallel })/2$ , respectively. Here $\boldsymbol{P}_e$ is the electron pressure tensor, $\boldsymbol{b}=\boldsymbol{B}/|\boldsymbol{B}|$ is the magnetic field unit vector and $T_e=\frac {1}{3}Tr(\boldsymbol{P}_e)/n_e$ is the isotropic electron temperature (with $Tr(\cdot )$ being the trace operator). In figure 4(d–f), we show the normalized deviation of electron perpendicular temperature ratio $T_{e, \perp }/T_{e, 0}-1$ ( $T_{e, 0}$ represents the initial isotropic electron temperature) in all three cases. We find perpendicular temperature increase ( $T_{e, \perp }/T_{e, 0} -1 \gt 0$ ) along the low-density separatrix. In contrast, perpendicular temperature decrease ( $T_{e \perp }/T_{e, 0}-1\lt 0$ ) is observed along the high-density separatrix. Additionally, we observe that the regions of perpendicular temperature increase and decrease narrow as the guide field strength increases.

Recent experimental studies found that E-REC with strong guide fields is accompanied by parallel heating along one separatrix (Shi et al. Reference Shi, Scime, Barbhuiya, Cassak, Adhikari, Swisdak and Stawarz2023), which is attributed to the parallel electric field. Here we show the normalized deviation of electron parallel temperature ratio $T_{e, \parallel }/T_{e, 0}-1$ from our three simulations in figure 4(g–i). Our results reveal parallel temperature increases along the high-density separatrix, consistent with experiments (Shi et al. Reference Shi, Scime, Barbhuiya, Cassak, Adhikari, Swisdak and Stawarz2023). However, we also observe parallel temperature decreases along the low-density separatrix. Furthermore, the parallel temperature increase and decrease are significantly influenced by the magnitude of the guide field. With increasing guide field strength, the parallel temperature increase/decrease regions become thinner and more localized, and the magnitude of the parallel temperature increase/decrease changes notably. It is also evident from the colourbar in figures 4(c–f) and 4(g–i) that the variations in the parallel temperature ratio $T_{e, \parallel }/T_{e, 0}$ exceed those in the perpendicular temperature ratio $T_{e, \perp }/T_{e, 0}$ , particularly in case C where the guide field is stronger. As a consequence, the parallel temperature ratio $T_{e,\parallel }/T_{e,0}$ shows more correlation with the spatial distribution of $A_e$ than perpendicular temperature ratio $T_{e,\perp }/T_{e,0}$ . Both the spatial distribution of $T_{e,\parallel }/T_{e,0}$ and $A_e$ illustrate a transition from a quadrupolar to a hexapolar to an octopolar structure as the guide field increases. Specifically speaking, in simulation A, $A_e$ shows a quadrupolar pattern, with two positive bands located on the low-density separatrix and the other two negative bands along the other separatrix and within the outflow regions. In simulation B, $A_e$ exhibits a hexapolar structure. In addition to the quadrupolar structure similar to that observed in figure 4(a), there is a faint band (where $A_e$ is almost zero near $x\approx 11, y\approx 16$ downstream of the separatrix, which is absent in simulation A. These two faint bands (there is also a faint band in the left-half of the domain) segregate the negative bands along the separatrix from those within the outflow regions, culminating in a pronounced hexapolar structure in both $A_e$ and $T_{e,\parallel }$ . In simulation C, an octopolar structure appears, which features two positive bands on one separatrix paired with two negative bands along the other separatrix. Additionally, there are two bipolar alternations of positive and negative temperature anisotropy bands inside the outflow regions and downstream of the separatrix, which is a unique characteristic in this case. Essentially, our analysis shows that E-REC produces different electron temperature anisotropy patterns under varied guide fields, with more parallel temperature increase/decrease features developing in magnetic islands as the guide field strength increases.

In order to exclude the possibility of the results being influenced by plasma recirculation, We have calculated the recirculation times for two conditions of the guide field, $ B_g = 10 B_{x0}$ and $ B_g = 20 B_{x0}$ (since the new features appear in these two cases).

1. (i) For $ B_g = 10 B_{x0}$ , where the box size $ L = 2 d_i$ and the maximum electron outflow velocity in the x direction $ V_{ex,\max } = 0.006c$ (c is the speed of light), the recirculation time $\Delta t$ is computed as

   (3.1) \begin{equation} \Delta t = \frac {L}{V_{ex,\max }} \approx 333 \omega _{pi}^{-1} = 666 \omega _{ce}^{-1}. \end{equation}
   Given that the reconnection evolution time from figure 1(b) is approximately $ 530 \omega _{ce}^{-1}$ (starting from when the reconnection rate is above zero (which is $\approx 250 \omega _{ce}^{-1}$ ) and ending at the time when we do analyses), it is evident that the observed structures are not a result of recirculation.

2. (ii) For $ B_g = 20 B_{x0}$ , with $ V_{ex,\max } = 0.007c$ , the calculated recirculation time is

   (3.2) \begin{equation} \Delta t \approx 286 \omega _{pi}^{-1} = 1143 \omega _{ce}^{-1}. \end{equation}
   Considering the evolution time of $ 1740 \omega _{ce}^{-1}$ from figure 1(c) (starting from when the reconnection rate is above zero (which is $\approx 700 \omega _{ce}^{-1}$ ) and ending at the time when we do analyses), and recognizing that the actual velocity is much lower than the maximum value used here, the real recirculation time would be considerably longer.

To provide more compelling evidence, we have conducted an additional simulation using a larger simulation box, specifically $L_x=6d_i, L_y=3d_i$ (case D), with a guide field of $B_g=20B_{x0}$ . This simulation yielded similar anisotropy results. In this case, the analysis was undertaken at a time corresponding to that depicted in figure 4(c). The electron temperature anisotropy results for case D are shown in figure 5. We also calculated the recirculation time for case D. We found the electron outflow jet velocities are comparable to that of case C. Since the box side lengths are three times larger, the corresponding recirculation times are also proportionally longer: $\Delta t \approx 3429 \omega _{ce}^{-1}$ for case D. This indicates that the recirculation time in this larger simulation box is significantly greater than the time at which positive anisotropy in the exhaust is observed (i.e. $ 1740 \omega _{ce}^{-1}$ as mentioned above). Notably, the recirculation time also exceeds the total analysis time (the time we do analyses, i.e. $t=2440\omega _{ce}^{-1}$ ), providing robust evidence supporting our argument. It should be noted that the choice of vertical length $L_y=3d_i$ in case D does not result in interactions between the two current sheets, as the reconnection along the lower sheet remains in an early phase and no interactions are observed at the time of analysis. Furthermore, the reconnection in case D remains within the electron-only regime. Consequently, we affirm the results are not caused by recirculation.

Figure 5. In case D, with a significantly larger simulation box ( $L_x = 6d_i$ , $L_y = 3d_i$ ), the electron anisotropy results align with those shown at the same time snapshot as figure 4(c). Similar patterns of anisotropy are observed, with positive $A_e$ (red regions) appearing along one separatrix and within the outflow regions.

3.3. Electron velocity distribution functions

In the previous section, we analysed the electron temperature anisotropy observed in the three cases, highlighting the notable variations that take place as the guide field strength increases. In this section, we further examine the EVDFs at different spatial locations. First, we analyse the EVDFs at the X-point in the three cases. We also focus on regions exhibiting pronounced anisotropy, where we observe a transition from a quadrupolar to a hexapolar to an octopolar structure as the guide field intensity increases. Figure 4 illustrates the selected positions where EVDFs are calculated, denoted by coloured square boxes whose side is approximately 0.3 $d_e$ , each containing approximately 500000 electron macroparticles. We only sample macroparticles within the bands in the right-half of the domain since the structure of the reconnection region is symmetrical, and we find analogous results in the left-half of the domain (not shown here).

Figure 6 presents the 2-D projections of the EVDFs at different locations in case A. Each row has four panels, with the initial three depicting the normalized microscopic electron velocity components centred to their mean values, i.e. $w_{\parallel }/v_{th,e}=(v_{\parallel }-v_{\parallel ,{\rm mean}})/v_{th,e}$ , $w_{\perp 1}/v_{th,e}=(v_{\perp 1}-v_{\perp 1,{\rm mean}})/v_{th,e}$ and $w_{\perp 2}/v_{th,e}=(v_{\perp 2}-v_{\perp 2,{\rm mean}})/v_{th,e}$ , respectively. The last panel in each row represents the one-dimensional (1-D) histogram of these three components with 256 bins. The velocity components are defined as follows (the same as Goldman, Newman & Lapenta (Reference Goldman, Newman and Lapenta2016) and Arrò et al. (Reference Arrò, Pucci, Califano, Innocenti and Lapenta2023)): $v_{\parallel }$ aligns with the local magnetic field $\boldsymbol{B}$ , $v_{\perp 1}$ points towards $\boldsymbol{B} \times \hat {z}$ and $v_{\perp 2}$ is oriented in the direction of $\boldsymbol{B} \times \boldsymbol{v}_{\perp 1}$ . Here $v_{\parallel ,{\rm mean}}, v_{\perp 1,{\rm mean}}$ and $v_{\perp 2,{\rm mean}}$ denote the fluid drift velocity in each direction and have been subtracted. Here $v_{th,e}$ represents the initial electron thermal velocity. The colourbar of 2-D plots and the y axis of the 1-D histogram indicate the number of macroparticles.

Figure 6. The EVDFs in case A are depicted as follows: (a–d) correspond to the magenta box (the X-point with negative $A_e$ ); (e–h) correspond to the yellow box (the outflow region with negative $A_e$ ); and (i–l) correspond to the cyan box (the separatrix in the lower right quadrant with positive $A_e$ ). The colourbar indicates the particle count. Panels (a,e,i), (b,f,j) and (c,g,k) display the projections of EVDFs in distinct planes, while (d,h,l) illustrate the 1-D EVDFs, with the y axis denoting the number of particles.

Figure 6(a–d) corresponds to the magenta box (X-point). Figures 6(a) and 6(b) indicate evident skewness in $v_{\parallel }$ , implying a net parallel heat flux density. This observation is supported by figure 6(d). We also find the agyrotropy is very small from figure 6(c) where the projection of the EVDFs on the $v_{\perp 1}-v_{\perp 2}$ plane is almost rotationally symmetric. Agyrotropy represents the asymmetry in the velocity distribution of particles within the plane perpendicular to the magnetic field. It is used to measure the degree to which the velocity distribution deviates from gyrotropy. The near rotational symmetry indicates that the EVDFs are close to gyrotropic in our simulations, though small finite agyrotropy also exists. Additionally, in figure 4(g), parallel heating is evident in this region, corroborated by the slightly irregular oval-shaped EVDFs in figure 6(a,b). Figure 6(e–h) depicts the EVDFs of the yellow box (outflow region). Figure 6(e,f) also display irregular oval-shaped EVDFs, consistent with the parallel heating observed in the outflow region in figure 4(g). Skewness is also apparent in figure 6(h), while figure 6(g) illustrates small agyrotropy in this region. Similar agyrotropic characteristics are observed in the cyan box (the separatrix in the lower right-hand quadrant), as depicted in figure 6(k). This region contrasts with the former two, exhibiting perpendicular heating and parallel cooling, as evidenced in figure 4(d,g). As a consequence, the EVDFs in figure 6(i,j) exhibit slightly irregular vertical oval shapes. Figure 6(l) indicates that the perpendicular heating is relatively small in this region since the broadening of the perpendicular velocities is not very pronounced with respect to the parallel component.

Figure 7 illustrates the 2-D projections of the EVDFs at four locations in case B. The first row corresponds to the magenta box (X-point), where the skewness similar to case A is observed in figures 7(a) and 7(b). The skewness is also evident in the 1-D histogram (figure 7 d). Figure 7(c) displays a distribution that is predominantly gyrotropic, with small but finite agyrotropy, which is also observed for EVDFs calculated inside the other three boxes (figures 7 g, 7 k and 7 o) in case B. Figure 7(e–h) shows the EVDFs of the yellow box. Due to the parallel heating in this area, the velocity distribution becomes broader in the parallel direction. Consequently, figures 7(e) and 7(f) display oval-shaped distributions, which are more pronounced than those observed in the yellow box of case A. The parallel heating is also highlighted in figure 7(h). On the other hand, the cyan box experiences parallel cooling and perpendicular heating. Therefore, the EVDFs appear elongated in the directions perpendicular to the local magnetic field, as demonstrated in figures 7(i) and 7(j). Figure 7(l) displays the perpendicular heating more distinctly as the two perpendicular components become wider. Figure 7(m–p) demonstrates the EVDFs of the lime box, which also presents prominent parallel heating. Correspondingly, figures 7(m) and 7(n) show that 2-D EVDFs are elongated in the parallel direction. However, their shapes are irregular, as also shown in figure 7(p), which is significantly different from the yellow box in case A (also located in the outflow region).

Figure 7. The EVDFs in case B are depicted as follows: (a–d) correspond to the magenta box (the X-point with negative $A_e$ ); (e–h) correspond to the yellow box (the separatrix in the upper right-hand quadrant with negative $A_e$ ); (i–l) correspond to the cyan box (the separatrix in the lower right-hand quadrant with positive $A_e$ ); and (m–p) correspond to the lime box (the outflow region with negative $A_e$ ). The colourbar indicates the particle count. Panels (a,e,i,m), (b,f,j,n) and (c,g,k,o) display the projections of EVDFs in different planes, while the (d,h,l,p) represent the 1-D EVDFs, with the y axis denoting the number of particles.

In figure 8, we show the 2-D projections of the EVDFs at five locations in case C. Figure 8(a–d) illustrate the EVDFs of the magenta box (X-point). Figures 8(a) and 8(b) suggest a skewed parallel velocity component, similar to the magenta boxes in cases A and B. Figure 8(d) represents the 1-D EVDFs that confirm this skewness, indicating a significant heat flux density near the X-point in case C. Figure 8(c) reveals that the EVDFs in the magenta box are close to gyrotropic, a feature that becomes more notable as the guide field increases as the core region becomes more rotationally symmetric in figure 8(c) compared with figures 6(c) and 7(c). In fact, the EVDFs of all five boxes show predominantly gyrotropic features in case C, as shown in figures 8(g), 8(k), 8(o) and 8(s). Remarkable anisotropy and parallel heating along the separatrix in the upper right-hand quadrant (the yellow box) is observed in figure 4(c,i). Figure 8(e–h) shows the EVDFs of the yellow box. Due to the parallel heating, oval-shaped non-Maxwellian EVDFs appear, as shown in figures 8(e) and 8(f). Figure 8(h) also clearly displays that the parallel velocity component becomes wider. It is also noted that the oval-shaped EVDFs become more regular and elongated in case C, compared with cases A and B. Correspondingly, the $T_{e, \parallel }$ of the yellow box in case C also becomes larger than in the other two cases, as shown in figure 4(g–i). Figure 8(i–l) depicts the EVDFs of the cyan box. Both parallel cooling and perpendicular heating exist along this separatrix in the lower right-hand quadrant, leading to the formation of vertically elongated EVDFs, as illustrated in figures 8(i) and 8(j). The parallel cooling is illustrated more clearly in the 1-D histograms as $w_{\parallel }/v_{th,e}$ becomes narrower (figure 8 l). The EVDFs in the cyan box, which are close to gyrotropic and exhibit small finite agyrotropy, suggest that the heating is close to isotropic within the plane perpendicular to the local magnetic field. Figure 8(m–p) illustrates the EVDFs of the lime box. Figures 8(m) and 8(n) show EVDFs that are stretched along the parallel direction, corresponding to the parallel heating observed in this region in figure 4(i). In addition, the skewness in figure 8(p) denotes relatively large heat flux density in this area. However, the EVDFs projected in the parallel direction become irregular. Figure 8(q–t) demonstrates the EVDFs of the purple box, which exhibit different features from the other boxes. A flat-top parallel velocity distribution is observed for particles with small random velocity (the velocity of particles relative to the bulk velocity) in the parallel direction, while parallel cooling emerges for particles with large random velocity in the parallel direction, as shown in figures 8(q), 8(r) and 8(t), where 2-D EVDFs appear stretched in perpendicular directions and compressed in the parallel direction.

Figure 8. The EVDFs in case C: (a–d) correspond to the magenta box (the X-point with negative $A_e$ ); (e–h) correspond to the yellow box (the separatrix in the upper right-hand quadrant with negative $A_e$ ); (i–l) correspond to the cyan box (the separatrix in the lower right-hand quadrant with positive $A_e$ ); (m–p) correspond to the lime box (the outflow region with negative $A_e$ ); (q–t) correspond to the purple box (the outflow region with positive $A_e$ ). The colourbar indicates the particle count. Panels (a,e,i,m,q), (b,f,j,n,r) and (c,g,k,o,s) are the projections of EVDFs in different planes. Panels (d,h,l,p,t) represent the 1-D EVDFs, with the y axis denoting the number of particles.

To explore the causes of the temperature anisotropy and the EVDFs observed along separatrices, we analyse the relationship among electron density $n$ , magnetic field strength $B$ and electron parallel temperature $T_{e \parallel }$ . According to the CGL theory, the following scaling relation holds (Bittencourt Reference Bittencourt2013):

(3.3) \begin{equation} \frac {T_{e \parallel } B^2}{n^2} = \text{constant}. \end{equation}

We plot $T_{e \parallel }$ against ${n^2}/{B^2}$ for each of the three cases to examine their alignment with CGL theory prediction or the Boltzmann condition (Chew, Goldberger & Low Reference Chew, Goldberger and Low1956; Snyder et al. Reference Snyder, Hammett and Dorland1997; Wetherton et al. Reference Wetherton, Egedal, Lê and Daughton2019). From the CGL relation, we can derive

(3.4) \begin{equation} T_{e \parallel } \propto \left (\frac {n^2}{B^2}\right ), \end{equation}

which means $T_{e \parallel }$ and ${n^2}/{B^2}$ have a positive linear relationship. In comparison, the Boltzmann response is a horizontal line, as $T_{e \parallel }$ remains constant regardless of $n$ and $B$ . The scatter plots of $T_{e \parallel }$ versus ${n^2}/{B^2}$ provide insight into whether the simulations adhere to the CGL predictions (Wetherton et al. Reference Wetherton, Egedal, Le and Daughton2021; Oka et al. Reference Oka, Birn, Egedal, Guo, Ergun, Turner, Khotyaintsev, Hwang, Cohen and Drake2023).

For each of the three simulations with varying guide field strengths, we select data points residing on the separatrix regions, which encompass both high-density and low-density areas. In order to also show their relationship under moderate density, we also cover the regions neighbouring to the separatrices. The regions we selected are shown in figure 4(a–c), represented by the cyan and yellow rectangles. From these regions, we extract values of electron density $n$ , magnetic field strength $B$ and electron parallel temperature $T_{e \parallel }$ . We then create scatter plots of $T_{e \parallel }$ versus ${n^2}/{B^2}$ for each simulation as shown in figure 9 (the three quantities ( $T_{e \parallel }$ , $n$ and $B$ ) are all normalized by their respective mean values).

Figure 9. Scatter plot of the relationship between electron parallel temperature $T_{e \parallel }$ and $(n/B)^2$ in three cases, where $n$ is electron number density and $B$ is magnetic field strength. The three quantities ( $T_{e \parallel }$ , $n$ and $B$ ) are all normalized by their mean value.

1. (i) Case A – ( $B_g = 5 B_{x0}$ ). At lower values of ${n^2}/{B^2}$ (corresponding to the low-density regions since the variations of $B^2$ are very slight in the domain), $T_{e \parallel }$ remains approximately constant, indicating Boltzmann-like behaviour. As ${n^2}/{B^2}$ increases, $T_{e \parallel }$ also increases modestly, indicating the onset of CGL effects. This suggests that the system predominantly exhibits Boltzmann-like behaviour at lower ${n^2}/{B^2}$ values, where electron heat conduction effectively maintains temperature isotropy. The slight increase in $T_{e \parallel }$ at higher ${n^2}/{B^2}$ indicates the onset of CGL relation, showing that adiabatic effects start to influence the electron temperature.

2. (ii) Case B – ( $B_g = 10 B_{x0}$ ). The relationship shows a more pronounced increase in $T_{e \parallel }$ with ${n^2}/{B^2}$ , exhibiting a transitional behaviour between Boltzmann and CGL effects. This indicates that the system is in an intermediate regime where both heat conduction and adiabatic effects significantly influence the electron temperature. The increased guide field strength enhances the influence of adiabatic processes, resulting in more noticeable CGL characteristics.

3. (iii) Case C – ( $B_g = 20 B_{x0}$ ). The scatter plot shows an almost linear relationship between $T_{e \parallel }$ and ${n^2}/{B^2}$ , suggesting a strong alignment with CGL theory and adiabatic processes dominating the electron temperature.

However, the presence of finite parallel heat flux in all simulations (as shown in figure 11) may be considered deviations from the ideal CGL behaviour, where zero heat flux is assumed. We note, nevertheless, that within the CGL model the heat flux is assumed to be zero as part of the mathematical closure of the fluid moment hierarchy. As a result, the heat flux density is not necessarily zero in fully kinetic simulations of a system approaching ideal CGL behaviour (see discussion in Allmann-Rahn, Trost & Grauer (Reference Allmann-Rahn, Trost and Grauer2018)). In this context, we do not interpret the presence of non-zero heat flux density in our simulations as a contradiction to CGL-like behaviour, and we consider that our three cases lie in a transitional regime where both thermal conduction and adiabatic effects coexist. It is also noted that the parallel heat flux intensity increases as the guide field strength enhances, which might be caused by the presence of multiple populations (Hwang et al. Reference Hwang2019; Lavraud et al. Reference Lavraud2021) and warrants further investigation. These findings provide an explanation for the electron parallel temperature variations observed along the separatrices. Since the parallel velocity component in the EVDFs must be consistent with the parallel temperature variations, they also elucidate the causes of the EVDFs’ shape along the separatrices we observed in all three cases.

The electric field causes a skewed velocity distribution function along its direction by accelerating or decelerating particles (Whealton & Woo Reference Whealton and Woo1972). In our simulations, the parallel electric field mainly resides near the X-point (here ‘parallel’ means the direction parallel to the local magnetic field), which accelerates electrons in the negative direction and causes the negative-skewed EVDFs features at the X-point in all three cases (here ‘negative-skewed’ means the parallel velocity components of the EVDFs at the X-point is skewed towards the negative direction, as shown in figures 6 d, 7 d and 8 d).

It is also noted that the selected yellow box for run A in figure 4(a) is not exactly adjacent to the separatrix, unlike those in runs B and C. This is because in run A, the region near the separatrix exhibited very limited variation in local density, while showing a significant density difference between the two boxes. As shown in figure 2(d), the density enhancement and cavity regions in run A occupy a relatively larger area compared with runs B and C. If we had chosen a box very close to the separatrix, the individual boxes would contain nearly uniform density values, and the span of $(n/B)^2$ would be narrow. Consequently, the resulting comparison between $T_{e\parallel }$ and $(n/B)^2$ would not sufficiently reveal the scaling trend expected from the CGL relation. To ensure a meaningful spread in $(n/B)^2$ , which is crucial for visualizing the CGL-like trend, we deliberately selected the current box placement. As a result, the deviation from CGL predictions in run A, as shown in figure 9(a), might be related to the box location. Nevertheless, we would like to emphasize that even in run A, a linear trend between $T_{e\parallel }$ and $(n/B)^2$ is still visible in the high-density portion (i.e. large $(n/B)^2$ values) of figure 9(a), indicating that CGL-like behaviour is still present to some extent.

In summary, by investigating the EVDFs in our three runs, we find that with the increase of guide field strength and the resultant reduction of electron gyroradius, the electron parallel velocity component changes remarkably (especially the electrons in the outflow regions). Different types of EVDFs are observed as the guide field strength increases, and the most various features are obtained in the case of the strongest guide field studied in this work. Notably, the EVDFs of all the examined locations in the three cases are close to gyrotropic. Figure 10 illustrates the agyrotropy in three cases (calculated using the formula in Swisdak (Reference Swisdak2016)) and its magnitude is very small. This indicates that the EVDFs are close to gyrotropic in such intense guide field conditions, while the agyrotropy could be large when the guide field is relatively low (Gao et al. Reference Gao, Tang, Guo, Li, Khotyaintsev, Graham, Turner, Yang and Wang2023).

Figure 10. Electron agyrotropy $G_e$ in all three cases. The three panels share the same colourbar, which is located on the far right-hand side of the figure.

Figure 11 illustrates the electron heat flux density in the parallel direction, normalized to the reference electron heat flux density $q_0=n_{0}\,m_e\,v_{th,e}^{3}$ (where $n_0$ is the initial density), for all three runs. These plots show the presence of finite positive parallel electron heat flux density at the X-point, extending to the separatrices as the guide field increases, which is consistent with the EVDFs observed above. Intense heat flux density is also evident in the outflow regions, albeit with a direction opposite to that observed in the aforementioned regions. As the guide field intensifies, the region exhibiting intense heat flux density becomes narrower, and the magnitude of the normalized heat flux density along the separatrices and within the outflow regions increases. Positive heat flux density is also present at the centre of the magnetic islands, but its amplitude decreases with increasing guide field strength.

Figure 11. Electron heat flux density in the parallel direction (with respect to the local magnetic field) in all three cases, normalized to the reference electron heat flux density $q_0 =n_{0}\,m_e\,v_{th,e}^{3}$ . The three figures share the same colourbar, which is located on the far right-hand side of the figure.