The role of guide field in magnetic reconnection driven by island coalescence

A number of studies have considered how the rate of magnetic reconnection scales in large and weakly collisional systems by the modelling of long reconnecting current sheets. However, this set-up neglects both the formation of the current sheet and the coupling between the diffusion region and a larger system that supplies the magnetic flux. Recent studies of magnetic island merging, which naturally include these features, have found that ion kinetic physics is crucial to describe the reconnection rate and global evolution of such systems. In this paper, the effect of a guide field on reconnection during island merging is considered. In contrast to the earlier current sheet studies, we identify a limited range of guide fields for which the reconnection rate, outflow velocity, and pile-up magnetic field increase in magnitude as the guide field increases. The Hall-MHD fluid model is found to reproduce kinetic reconnection rates only for a sufficiently strong guide field, for which ion inertia breaks the frozen-in condition and the outflow becomes Alfvenic in the kinetic system. The merging of large islands occurs on a longer timescale in the zero guide field limit, which may in part be due to a mirror-like instability that occurs upstream of the reconnection region.


I. INTRODUCTION
Magnetic reconnection is the process of changing magnetic field-line connectivity in highly conducting plasmas, 1,2 and is usually associated with the conversion of a large amount of magnetic energy into plasma kinetic energy. The question of how fast magnetic flux can be reconnected in large and weakly collisional systems is important in space weather modelling, 3,4 magnetic confinement fusion devices, [5][6][7][8] and the understanding of many astrophysical phenomena. A significant number of nonlinear simulation studies have considered this question with both fully kinetic and a number of reduced fluid models. One key result from these studies was that the Hall-MHD fluid model, in which ion inertia is responsible for setting the ion diffusion region thickness, was able to reproduce the reconnection rates of fully kinetic simulations. 9 The reconnection rate did not depend on the detailed kinetic physics of electrons or ions, [9][10][11][12] since this physics is absent in the Hall-MHD model.
The majority of these studies used the Harris-sheet setup, 13 in which a pre-formed current sheet is often highly unstable due to its kinetic scale thickness and macroscale length. Further, such a system does not include important features that are present in real reconnecting systems, such as the formation of the current sheet, and its coupling to a large scale system that supplies the magnetic flux to the reconnection site. The magnetic island coalescence problem [14][15][16][17][18][19][20] naturally includes such features, since the flux is supplied by two macroscale magnetic islands and a current sheet forms selfconsistently between the islands as they collide. Such a) Electronic mail: stanier@lanl.gov magnetic islands are two-dimensional representations of magnetic flux ropes, a fundamental building block of magnetised plasmas. [21][22][23] The interaction of islands, or flux-ropes, is thought to be a common process in the solar corona, the solar wind and the Earth's magnetosphere. In particular, island coalescence has been inferred both by indirect observations of coronal mass ejections, 24,25 and by in-situ detection in the magnetotail. 26 Recent simulation studies 20,27,28 have revisited the conclusions from earlier current layer studies using the island coalescence problem set-up. It was found that the Hall-MHD model was unable to reproduce the geometry of the ion diffusion region, the upstream magnetic field 'pile-up' strength, the outflow velocity, or the reconnection rate including its dependence on ion temperature or system size. 27 Further, these differences lead to different global evolution of the system -for Hall-MHD the islands merge when they first approach, but for the fully kinetic model the islands bounce off each other to give a much longer timescale for merging. Ion kinetic physics, and in particular the anisotropic and agyrotropic nature of the ion pressure tensor, was identified as the crucial physics missing from the Hall-MHD model. A hybrid model with mass-less fluid electrons and kinetic ions was able to reproduce many of the full kinetic results.
In the present study, the sensitivity of these results to the addition of a finite guide field is considered. Such a guide field is not reconnected, but it can modify the orbits of ions and electrons around the reconnection site. It is found that the reconnection rate, pile-up magnetic field, and outflow velocity increase with guide field for a limited range of guide field strengths that are applicable to the Earth's magnetosphere. The Hall-MHD model is able to reproduce the fully kinetic reconnection rate only in the limit of strong guide fields, when the ion pressure tensor physics becomes sufficiently localised around the x-point and the ion inertia is responsible for setting the ion diffusion region thickness. For large systems with kinetic ions, it is found that the time for the islands to fully coalesce is significantly shorter when there is a guide field present, compared with the zero guide field case.
The paper is organised as follows. In Sec. II, a description of the different plasma models used in this study is presented, and Sec. III gives the detailed parameters and initial conditions for the island coalescence simulations. The detailed role of the guide field in small systems is discussed in Sec. IV, and then additional effects in large systems are given in Sec. V. Finally, the important results are summarised and conclusions are drawn in Sec. VI.

II. DESCRIPTION OF PLASMA MODELS EMPLOYED
A key aim of this study is to evaluate the fidelity of several commonly used reduced physics models for the island coalescence problem with a guide field, by comparing them against fully kinetic Particle-In-Cell (PIC) simulations. A description of these models follows, along with a very brief description of the main numerical methods used.

A. Fully kinetic particle-in-cell
The fully kinetic simulations are carried out using the high performance electromagnetic relativistic particlein-cell (PIC) code VPIC. 29,30 The electric and magnetic fields are advanced with the Maxwell-Ampere and Maxwell-Faraday equations Here E is the electric field, B is the magnetic field, j is the current density, c is the speed of light and µ 0 is the vacuum permittivity. The particles are advanced using where u s,m = γ s,m v s,m /c is the normalised momentum of the particle m from species s, x s,m is the particle position, q s and m s are the charge and rest-mass for particles of species s = i, e for ions and electrons, and γ s,m = 1 + u 2 s,m . VPIC advances fields on a Yee-mesh with a finitedifference time-domain method and particles are pushed with a Boris algorithm modified to use a sixth order rotation angle approximation. The electromagnetic fields are interpolated to the particles using an energy conserving scheme, and the current density is gathered using the charge conserving Villasenor-Buneman method. See Refs. [29], [30] and references therein for further details.

B. Hybrid model
The hybrid and Hall-MHD models solve for the magnetic field advance using Eq. (2), but calculate the current density as j = ∇×B/µ 0 . These models also enforce charge neutrality as given for a hydrogen plasma by equal ion and electron number densities, n = n i = n e , and assume mass-less electrons, m e = 0. The electric field is then found from Ohm's law where p e is the scalar electron pressure, η is the resistivity, η H is the hyper-resistivity and e = |q e | is the elementary charge. The hybrid model treats ions kinetically using the PIC method, while electrons are treated as a mass-less fluid. The ions are advanced with Eqs. (3) and (4) in the nonrelativistic limit (γ s,m → 1), and the density n and bulk ion velocity v are calculated from velocity moments of the particle distribution. The electrons are assumed isothermal p e = T e0 n, where T e0 is the initial temperature. Eq. (2) is solved with the Runge-Kutta method and particles are pushed using the standard Boris method. See Ref. [31] and references therein for further details.

C. Hall-MHD model
The Hall-MHD fluid model solves additional fluid conservation equations for plasma number density, total momentum, and the isotropic total pressure Here, a constant temperature ratio approximation α = T i /T e = const. is used to give a single equation for the total pressure p = p i + p e = n (1 + α) T e carried by the effective velocity v * = v − j/[ne(1 + α)]. Simple heat conduction and heating terms are used, q = −κ∇T e and Q = ηj 2 +η H (∇j) 2 − ↔ Π : ∇v, respectively, where the ion viscous tensor is given by ↔ Π = −µ∇v. µ is the viscosity, κ is the heat conductivity, and γ is the heat capacity ratio.
The system of equations (2), (5), (6)-(8) is advanced using a fully implicit, preconditioned Newton-Krylov method. 32,33 Spatial discretisation is done using a cell centred finite volume scheme, 34 and the second order BDF-2 method is used for time advance.

D. 10-moment model
The 10-moment model 28,35 solves the equations for s = i, e, in addition to Eqs. (1) and (2) using the current calculated with j = s n s q s v s . Here, ↔ P s and ↔ Q s are the pressure and heat flux tensors for the species s in its rest-frame. Currently, the set of moment equations is truncated by the simple closure 28,35,36 where p s = trace[P s ]/3 is the isotropic pressure, and v T s is the thermal velocity. This relaxes the pressure tensor to an isotropic pressure at a rate given by |k s |v T s , allowing deviations from isotropy in the ion and electron pressure tensors at length scales less than 1/|k i | and 1/|k e | respectively. These k s are currently treated as free parameters, but work is being done to close the set of equations in the manner of Ref. [37]. These equations are solved using a second order locally implicit operator splitting approach, where the hyperbolic part uses a dimensionally split finite-volume wave propagation scheme. More details can be found in Refs. [28] and [35].

III. PROBLEM SET-UP
The simulations described in this paper are initialised with the exact Vlasov magnetic island equilibrium. 4,38 The magnetic field is given by B = ∇ × (A zẑ ) + B gẑ , where B g is the guide field strength and A z is the outof-plane magnetic potential given by Here, = 0.4 and B 0 is the asymptotic magnetic field. For equilibrium, it is required that the sum of initial ion and electron temperatures is T i0 + T e0 = B 2 0 / (2µ 0 n 0 k B ), where n 0 is the reference density equal to the density enhancement at the centre of the current layer in the limit = 0, and k B is the Boltzmann constant. In this paper the ratio of initial temperatures used is T i0 /T e0 = 1. The density profile is given by where n b = 0.2n 0 is a background density. The size of the simulation domain is For all codes we use perfectly conducting boundaries that are reflecting for particles at x = ±πλ, and periodic boundaries in the y-direction. We use a sinusoidal perturbation of magnitude 0.1B 0 to start the merging. Additional parameters specific to each model are as follows. For fully kinetic PIC simulations we use a ratio of electron plasma to gyro-frequency ω pe /Ω ce = 2, and mass-ratio m i /m e = 25. For the hybrid model, we use a ratio of the ion frequencies ω pi /Ω ci = 2000, a resistivity η = 10 −5 µ 0 d i v A0 , and a hyper-resistivity This value of hyper-resistivity is an order of magnitude smaller than in Ref. [27] for the hybrid model. The value η H = 10 −3 µ 0 d 3 i v A0 was found to be sufficient in the zero guide field case, but with finite guide field the ion diffusion region is thinner (see below) and a smaller η H is necessary to ensure good separation between ion and electron scales. This change in η H gives only a 2% reduction of the peak reconnection rate for the zero guide field run with λ = 5d i reported in Ref. [27].
For Hall-MHD, we use a fixed temperature ratio α = 1, ion viscosity µ = 10 −2 m i n 0 d i v A0 , resistivity η = 10 −5 µ 0 d i v A0 , hyper-resistivity η H = 10 −4 µ 0 d 3 i v A0 , and heat conduction κ = 10 −4 n 0 d i v A0 . Finally, for 10moment simulations we use ω pe /Ω ce = 2, m i /m e = 25, k e = 1/d e and k i = 1/3d i , where d e = m e /m i d i is the electron-skin depth. As for Harris-sheet geometry, we find that the precise value of the dissipation and m i /m e do not influence the reconnection rate, provided there is a sufficient separation between ion and electron scales.
To compare with previous studies, 20,27,28 we normalise the magnetic fields by the maximum value of the initial in-plane field in a line joining the island O-points, B m = 0.353B 0 . Velocities are normalised by v Am = B m v A0 /B 0 , and lengths are normalised by d i . Times are given in terms of the global Alfvén crossing time τ A = 4πλ/v A0 to easily compare between different system sizes. Finally, reconnection rates are measured as In this study, we vary the strength of the guide field B g , and the half-thickness of the equilibrium current layer λ. While varying λ, we keep fixed and the domain size proportional to λ. Thus, variations in λ are equivalent to varying the system size with respect to d i .

IV. EFFECT OF GUIDE FIELD ON RECONNECTION RATE
A. Reconnection rate Figure 1 shows how the peak reconnection rate E R , the aspect ratio of the ion diffusion region δ i /w i , the normalised pile-up magnetic field B in,i /B m , and the normalised peak outflow velocity v out,i /v Am scale with normalised guide field B g /B m for runs with λ = 5d i for the different plasma models. Here, δ i and w i are measured as the full-width half-maxima of the non-ideal electric field E z =ẑ · (E + v × B) in cuts across the ion diffusion region with x = 0 and y = 0 respectively. B in,i is the maximum value of the upstream magnetic field, and v out,i is the maximum outflow velocity. For the Harris-sheet problem, B g is often normalised by the asymptotic field B 0 . For our normalisaton, the value of As discussed in Ref. [27], for B g = 0 the Hall-MHD model is unable to reproduce the equivalent fully kinetic PIC values of the quantities plotted in Fig. 1. The simulations which retain the full physics of kinetic ions (hybrid and PIC) have a broader ion diffusion regions and show significantly reduced rates, pile-up and outflow velocities with respect to Hall-MHD. In particular, the PIC simulation with B g = 0 has a peak outflow velocity 7 times smaller than Hall-MHD, an aspect-ratio 3.2 times larger, and the pile-up at 80% of the Hall-MHD value. To illustrate how these differences modify the rate E R , it is useful to consider a quasi-steady and rectangular ion diffusion region (iDR) with constant inflow and outflow velocities along the edges. The reconnection rate is then set by the inflowing flux to the iDR v in,i B in,i , where the inflow velocity is constrained by the continuity equation v in,i ≈ v out,i (δ i /w i )(n out /n in ). For the numbers above, and with little change in n out /n in between Hall-MHD and PIC, this simple estimate gives a factor of 2.7 reduction in the rate for PIC that is broadly consistent with the ≈ 2 times smaller E R in Fig. 1.
For the range of B g in Fig. 1, the plotted quantities for the Hall-MHD simulations have weak guide field dependence. In particular, the peak rate E R decreases by only 13% from B g = 0 to B g /B m = 2.83. In contrast, the peak rates increase by 51% (53%) for PIC (hybrid) across this range to give E R = 0.665 (0.681) for B g /B m = 2.83 that are comparable to Hall-MHD (E R = 0.697). As far as we are aware, such an increase in the reconnection rate with guide field has been noted previously only in asymmetric reconnection, 39,40 where there are different thermal and magnetic pressures on either side of the current sheet. Symmetric current layer studies typically find that the reconnection rate is independent of B g for B g ≤ B 0 , then begins to decrease 41,42 as B g increases until it flattens again once the sound Larmor radius ρ s = c s /Ω ci , defined with the sound speed c s and ion cyclotron frequency Ω ci , falls below the electron-skin depth. 43,44 For larger guide fields, such as B g /B m = 5.66 shown in Fig. 1, E R begins to gradually decrease in the same way for Hall-MHD, hybrid and PIC. However, the focus of the present study is on the range of weak guide fields for which there is qualitatively different behaviour from the extended current layer simulations (B g /B m ≤ 2.83). This regime has application in magnetospheric plasmas.
In addition to the increase in E R , hybrid and fully kinetic PIC simulations show strong dependences on B g for δ i /w i , B in,i /B m and v out,i /v Am , but all the quantities are close to the Hall-MHD results at B g /B m ≥ 2.83.
Finally, as seen in Fig. 1, the 10-moment model is able to reproduce the same trends as fully kinetic PIC in all of the plotted quantities for the range of guide fields considered. In particular, both the rates E R and outflow velocities v out,i /v Am are in excellent agreement between the 10-moment and fully kinetic PIC simulations.

B. Ion diffusion region physics
The decrease in the aspect ratio δ i /w i with the increase in B g , for all of the models except Hall-MHD, follows from a significant decrease in the thickness of the ion diffusion region δ i . The length of this region w i also decreases with guide-field but by a smaller amount (not shown). Figure 2 shows the contributions to the nonideal electric field E z in 1D cuts across the ion diffusion region in the inflow direction (x = 0) for fully kinetic PIC simulations (top four panels) with λ = 5d i and varying guide field. The contributions to E z are given by the ion momentum equation, which in a normalised form is given by where the collisionless PIC simulations only have contributions from the ion inertia (blue) and ion pressure tensor (green) terms.
As mentioned above, δ i is measured as the full-width at half-maximum (FWHM) of the E z curve (black) in each case. This thickness decreases from δ i = 2.8d i for B g = 0 to δ i = 0.85d i for B g /B m = 2.83. For B g = 0 the ions exhibit characteristic meandering-type orbits as they cross the weak-field region, 10,12,27,45 giving rise to strong gradients in the off-diagonal elements of the ion pressure tensor. The extent of the ion orbits is reduced by the addition of a finite guide field, and the orbits may become chaotic if the gyro-radius becomes comparable to the radius of magnetic curvature, R c . 46 For a strong enough guide field, the gyro-radius will fall below both the meandering length and R c , and the ions start to become magnetised. Fig. 2 shows a transition at B g = 2.83B m to ion inertia (blue) supporting E z at the edges of the iDR, and thus setting the thickness. The pressure tensor term (green) continues to balance E z close to the stagnation point (y = 0) where ion inertia is small provided the diffusion region is quasi-steady, ∂ t (nv iz ) ≈ 0. For the island coalescence problem this quasi-steady phase is short in duration and occurs at the time of the peak reconnection rate. 47 For comparison, the contributions to E z are also plotted from the Hall-MHD simulation with λ = 5d i and B g = 2.83B m in the bottom panel of Fig. 2. Here, as in Ref. [27], P iz = −µ∇v iz is a collisional viscosity (green curve) and F coll,z = ηj z − η H ∇ 2 j z is the collisional friction (red curve), which is almost entirely from the hyperresistive term. The thickness of the ion-diffusion region for the Hall-MHD run with B g /B m = 2.83 is δ i = 0.56d i (35% smaller than the PIC value), which is reduced only slightly from the zero-guide field case (δ i = 0.62d i , not shown). In the Hall-MHD simulations, ion inertia supports E z at the edges of the ion diffusion region and sets the thickness δ i for the full range of guide-fields of Fig. 1. In this sense, there is qualitative agreement in the physics breaking the frozen-in condition for ions between the Hall-MHD and PIC models only for the strongest guide field case, B g /B m = 2.83. For Hall-MHD, the viscous and hyper-resistive sub-layers in E z are significantly thinner than the pressure tensor sub-layer in the B g = 2.83B m PIC result. Their thickness depends on the values of the dissipation coefficients µ and η H ; however, the reconnection rate is insensitive to their precise thickness provided they are thin in comparison to δ i . In the absence of a guide field and for a small system size (λ = 5d i ), it was demonstrated in Ref. [28] that the 10-moment model was able to reproduce either the Hall-MHD or kinetic values for E R and δ i , depending on the value chosen for the parameter k i . For k i = 1/d e the pressure tensor is allowed to depart from isotropy only on scales smaller than the electron skin-depth. In this case, the ion frozen-in condition is violated at d i -scale due to the ion inertial term. However, setting k i = 1/δ i , where δ i = 3d i is the approximate thickness of the ion diffusion region for the zero guide field PIC simulations, allows the ion pressure tensor to break the frozen-in condition on this larger scale. For B g = 0 the thickness of the iDR for the 10-moment model is indeed at this larger scale, δ i = 2.6d i , but for the stronger guide field case B g /B m = 2.83 the thickness reduces to δ i = 0.952d i and the contributions to E z are in qualitative agreement with the equivalent PIC result (fourth panel of Fig. 2). This suggests that the k i parameter only determines the maximum allowable thickness of the ion diffusion region for the 10-moment model, δ i ≤ 1/k i .
The modification of the ion kinetic physics at the Xpoint with guide field can be seen in the ion distribution functions, shown in Figure 3 from the same fully kinetic PIC simulations. These distribution functions were calculated by collecting ions in a square box with 2d e edge length, centred on the X-point and time averaged over 1Ω −1 ci . For B g = 0 there is a clear double peaked structure along the v y axis (the inflow direction), which is characteristic of the ion meandering orbits. This structure persists for weak guide-field B g = 0.28B m , but for B g = 1.13B m and B g = 2.83B m the distribution starts to become more gyrotropic about the guide field direction.

C. Pile-up magnetic field strength and outflow velocity
It is important to note that the reduction in δ i with guide field alone can not explain the increase in the reconnection rate with B g /B m from the PIC simulations in Fig. 1. For the estimate of the rate based upon a quasisteady rectangular ion diffusion region above, the rate would be expected to decrease as E R ∝ δ i . However, it is clear from Fig. 1 that there are changes in the other quantities that influence the reconnection rate.
Firstly, the strength of the pile-up magnetic field B in,i /B m increases with guide field for PIC, hybrid and 10-moment models until it reaches the Hall-MHD value for B g /B m ≈ 1.132. The strength of this pile-up field is set by the cross-scale coupling between the global motion of the interacting islands and the micro-scale physics of the diffusion region.
To interpret this increasing trend with guide field, it is helpful to consider previous studies of island coalescence. In resistive MHD 15 it was found that the pile-up increases as the current sheet thickness decreases (with increasing Lundquist number) until it saturates at a critical value set by pressure balance considerations. 17,48 With the inclusion of the Hall effect, the magnetic field is frozen-in to the faster electron flow below d i -scale and the pile-up is reduced. 16,49 However, recent studies 27,28 using compressible Hall-MHD, 10-moment, hybrid, and fully kinetic models have found that pile-up can still occur in larger systems just upstream of the ion diffusion region, and the strength of pile-up increases as the system-size becomes large with respect to the iDR thickness δ i . In Fig. 1 the increase in pile-up with guide field is consistent with the increase in system-size with respect to δ i , as δ i decreases significantly over this range of guide fields for all models except Hall-MHD.
The increased pile-up gives a larger upstream Alfvén speed and so may increase the outflow speed, if the outflow is accelerated by the relaxing magnetic tension in newly reconnected field lines. For Hall-MHD the outflow speed is indeed comparable, see Fig. 1, to the effective upstream Alfvén speed v out,i ≈ v Am (B in,i /B m ) n 0 /n set by the pile-up field B in,i /B m and upstream density n/n 0 ≈ 1.4 on the edge of the iDR. However, as mentioned above, the kinetic ion codes have an outflow velocity ≈ 7 times smaller for B g = 0. Also, unexpectedly, the ion outflow speed is found to increase significantly with guide field.
To understand this effect we consider the momentum balance through a wedge shaped region embedded within the outflow, following the method used in Ref. [50]. The wedge region shown in Fig. 4 is found by integrating along the separatrix field-lines on one side of the X-point to y = ±y max , then the contour is closed by joining the ends with a straight line intersecting the ion outflow jet.
Provided the exhaust is quasi-steady, where l is the closed contour bounding the wedge region of area S,n is the outward normal unit vector at a point on the contour line, is the total momentum tensor, p T = p ⊥i +p ⊥e +B 2 /2µ 0 is the total of the perpendicular thermal and magnetic pressures, and The magnetic tension term BB has non-zero contribution only at the straight end piece of the contour, so to compare the tension more fairly between fully kinetic PIC and Hall-MHD runs we choose the length of the end piece to be the same in all calculations, y max = 1.6d i . This value is chosen so that the end contour intersects the outflow jet in both Hall-MHD and PIC. For the Harris-sheet set-up considered in Ref. [50] there is a significant quasisteady period in which time derivative terms in the outflow jet region can be assumed small but, as mentioned above, for island coalescence this assumption is valid only at the peak rate. Figure 5 shows the relative contributions from the different terms in T ↔ to the contour integral in Eq. (16) from fully kinetic runs with λ = 5d i , B g = 0 (top) and B g /B m = 2.83 (middle), and the Hall-MHD simulation with λ = 5d i and B g /B m = 2.83 (bottom) at the time of the peak reconnection rate.
For the fully kinetic PIC simulation with B g = 0 the inertial term is weak (blue) and the combined magnetic tension (green) and total pressure (gold) is almost entirely balanced by the stress tensor term (red). This results in an outflow velocity being much less than the upstream Alfvén speed in this case. Within the stress tensor contribution the majority (85%) of this value is from the non-gyrotropic terms, of which the ion contribution is largest. Thus, the mechanism for which the outflow jet is slowed appears distinct from the pressure anisotropy mechanism discussed in Refs. [50] and [51].
There is a significant increase in outflow inertia for the PIC simulation with B g /B m = 2.83. However, this can not be explained fully by the increase in magnetic tension (12% increase) and the reduction in the stress tensor term (28% decrease, main contribution from ion non-gyrotropic part), as the main contribution is due to a factor of three increase in p T between B g = 0 and B g /B m = 2.83. This total presssure term cancels out a significant proportion of the stress tensor term, so that the magnitude of the inertial term is comparable to that of the magnetic tension. This is consistent with an outflow velocity close to the upstream Alfvén speed. In contrast to the kinetic result, the total pressure and stress tensor terms are much smaller in the Hall-MHD run with B g /B m = 2.83. The Hall-MHD runs have isotropic pressure (p i,e = p ⊥i,e ) such that the stress tensor term is due to the collisional ion viscosity The viscous stress tensor term again mostly cancels with the total pressure so that the inertial term is balanced by magnetic tension and the outflow is Alfvénic.
It is clear that there are differences in the outflow physics between the fully kinetic PIC and Hall-MHD result for B g /B m = 2.83. Such differences may be related to the thickness of the ion pressure tensor sub-layer in the non-ideal electric field, shown in the bottom two panels of Fig. 2. However, these differences do not have a strong influence on the magnitude of the outflow velocity or the reconnection rate between the two runs.

A. Global motions and average rates
The results presented thus far have focussed on the smallest system-size λ = 5d i , for which the peak rate E R differs by a factor of two between Hall-MHD and fully kinetic PIC for B g = 0. However, as was demonstrated in Ref. [27], such order unity differences in E R can lead to different global evolution of the system. For zero guide field with λ ≥ 10d i , the separation of the island O-points in hybrid and fully kinetic PIC simulations does not monotonically decrease with time -the islands exhibit "sloshing" behaviour. 15,17,48 When this occurs, reconnection can temporarily shut off in the periods of time that the islands move apart, resulting in a slower average reconnection rate and a longer timescale for the islands to fully merge. 27 Figure 6 shows the separation between the island Opoints L sep , normalised by the initial separation L 0 , versus global time for λ = (5 − 25)d i and B g /B m = 2.83. The zero guide field result for the largest system size (λ = 25d i ) is also plotted for reference. 27 Although there is reversal in the O-point separation for λ ≥ 10d i , this is significantly less than in the zero-guide field case as shown for the λ = 25d i runs. Figure 7 shows the dependence of the average reconnection rates < E R > on system size, where <> denotes the time average over 1.5τ A . For B g = 0 (top panel) the hybrid and fully kinetic PIC simulations that include the full kinetic ion physics have an average rate that decreases steeply with system size, < E R >∝ (λ/d i ) −0.8 for the PIC simulation, and < E R >∝ (λ/d i ) −0.65 for the hybrid simulation. 27 This is in sharp contrast to the Hall-MHD run with B g = 0 which has an average rate of < E R >∝ (λ/d i ) −0.25 . The 10-moment simulations agree with the fully kinetic runs for the smallest system-size (λ = 5d i ) but overestimate < E R > in larger systems, where < E R >∝ (λ/d i ) −0.2 , see below and Ref. [28]. B. Secondary island formation Figure 8 (Multimedia view) shows the current density and magnetic flux from the fully kinetic PIC simulations with λ = 25d i for the case of B g = 0 (left panel) and B g = 2.83B m (right panel). When the islands collide there are significant qualitative differences in the structure of the current layer, the upstream plasma within the islands, and the outflow jets. In contrast to the zero guide-field run, the simulation with B g /B m = 2.83 has a thinner current sheet with larger current density that is unstable to the formation of secondary magnetic islands. The repeated formation and ejection of such islands is thought to be important in regulating both the length of the layer and the rate of collisionless reconnection, 52 and secondary islands/plasmoids have been well studied in the limit of resistive MHD. [53][54][55][56] A complete analysis of the dynamics of these secondary islands is beyond the scope of the present study, but we report several interesting observations from these simulations. Firstly, in the case of zero guide-field we see no secondary magnetic island formation in either the Hall-MHD, hybrid or kinetic simulations for λ = 5 − 25d i . However, this does not exclude secondary islands in larger systems, and we have indeed observed a solitary island at late time (t ≈ 1.8τ A ) in hybrid and fully kinetic PIC runs with λ = 50d i . Secondly, as discussed in Ref. [28], the current layer in the zero-guide field 10-moment simulations is unstable to the formation of secondary islands for simulations with λ ≥ 10d i . These secondary islands do appear to influence the reconnection rate, and only the smallest simulation (λ = 5d i ) without islands is found to reproduce kinetic reconnection rates -see Fig. 7 and Ref. [28]. For the simulations with B g = 2.83B m , the formation of secondary islands differs between the models. For fully kinetic PIC simulations such islands are present for λ ≥ 10d i , but for Hall-MHD no islands are observed in the range 5 − 25d i and the current layer has an open X-point configuration. We note that Ref. [57] has demonstrated such an open X-point configuration is more favourable to the efficient release of magnetic energy in a reduced two-fluid model with electron diamagnetic and inertial effects. For the hybrid simulations, we also observe the open X-point configuration, and the absence of secondary magnetic islands in runs with λ = 5 − 15d i . For λ = 25d i there are a small number of islands formed in periods where the x-point is not fully opened up. Finally, for the 10-moment simulations, secondary islands are formed in the simulations with λ ≥ 15d i . The differences in secondary island formation between these models do not strongly influence the value of < E R > for the guide field simulations, as similar values are found for all models (the bottom panel of Fig. 7).

C. Instability within pile-up region
It was reported in Ref. [27] that both hybrid and fully kinetic PIC simulations with zero guide-field have finite pressure anisotropy (p ⊥ /p = 1) within the upstream flux pile-up region at the time that the islands first collide. For the range of system sizes considered, both the maximum pile-up magnetic field strength, see Fig. 9, and the maximum pressure anisotropy increase with system size. Figure 10 shows the ion pressure anisotropy p ⊥i /p i from a large fully kinetic PIC simulation with λ = 50d i and B g = 0. There is a visible long wavelength modulation of p ⊥i /p i within a wedge shaped region of increased pressure anisotropy. A similar modulation is also visi- ble in the current density from the PIC simulation with λ = 25d i and B g = 0 shown in Fig. 8 (Multimedia view). Although the instability occurs in a region with finite gradients in both plasma density and magnetic field strength, we note that these gradients are larger for smaller systems in which no instability is observed. For instabilities driven by the free energy associated with the anisotropic pressure in this region (p ⊥ /p > 1), there are three main candidate modes. In the linear regime, the mirror instability is obliquely propagating, has zero real frequency, and both ion and electron anisotropies contribute to the drive. 58 The other candidates are the ion and electron cyclotron anisotropy instabilities (ICAI and ECAI), driven by the ion and electron pressure anisotropies respectively. These propagate parallel to the field in the linear regime, and have finite real frequencies. Fig. 10 shows the mode at around the time of nonlinear saturation. The mode structure at this time is generally oblique to the local magnetic field direction. We also observe that the mode appears stationary with respect to the frozen-in field of the moving islands.
To further isolate the mode, we solve numerically the electromagnetic Vlasov linear dispersion relation based upon plasma parameters in the wedge shaped region at the time before the instability is visible. For simplicity, we assume a straight and uniform magnetic field, and a uniform bi-Maxwellian plasma velocity distribution. The assumption of a uniform plasma and magnetic field is not strictly valid for the region of interest, and so the calculated values should only be considered approximate.
Typical spatially averaged values of the plasma parameters within the wedge shaped region during the time leading up to the instability are p i⊥ /p i = 1.4, p e⊥ /p e = 1.1, β i = 15, β e = 15, and v A /c = 0.0165 with β i,e = 2µ 0 p i,e /B 2 . For these parameters, the fastest growing mode of the three is the mirror mode with growth rate γ mirror /Ω ci = 0.06, wavenumber kd i,eff = 0.32 and oblique angle to the magnetic field θ = 60 • , where d i,eff = d i n 0 /n is the effective ion skin depth based upon the local density. The predicted wavelength of the linear mode λ mirror = 19.6d i,eff is only slightly larger than the typical measured wavelength (λ ≈ 12d i = 14.2d i,eff based upon the local density n ≈ 1.4n 0 ). It should be noted that these measured values are taken at non-linear saturation of the instability when the modulation is clearly visible above the background PIC noise, and the fastest growing mode from linear Vlasov theory is not necessarily the largest amplitude mode at saturation. For the stated plasma parameters, we find the ECAI is stable and the ICAI has smaller growthrate γ ICAI /Ω ci = 0.03 and longer wavelength kd i,eff = 0.16.
The instability in Fig. 10 gives rise to local depressions (troughs) in p i⊥ /p i . It also leads to strong peaks in the magnetic field strength, where the magnetic flux contours in Fig. 10 are close together, in similar locations to the p ⊥i /p i troughs as well as local density depressions. Such magnetic peaks are a common feature of mirror unstable plasmas in the magnetosphere, see, e.g., Ref. [59].
The peaks in the upstream field give rise to a significant increase in B max /B m for the zero guide-field PIC simulation with λ = 25d i in Fig. 9, since this is measured simply as the maximum upstream magnetic field. There is a much weaker increase in pile-up for the hybrid run with λ = 25d i , where signatures of this instability are very weak. The absence of the instability in this hybrid run may be due to the use of isotropic electron pressure, which means that only ion pressure anisotropy can contribute to the drive. For example, with the plasma parameters stated above but with isotropic electrons p e⊥ /p e = 1 the growth rate is reduced by almost a factor of 2 to γ mirror /Ω ci = 0.034. For a larger hybrid run with λ = 50d i we do see structures in p i⊥ /p i and the magnetic field that are very similar to those in Fig. 10. The difference in system-size threshold at which the instability is clearly observed between hybrid and PIC also coincides with a greater difference in the reconnection rate between the two codes. In particular, Fig. 7 shows a significantly lower average reconnection rate in the λ = 25d i and B g = 0 PIC run where the instability is present, compared with the hybrid result for the same system size. If the average reconnection rate for the B g = 0 PIC simulations is calculated neglecting the λ = 25d i run, the system-size scaling is < E R >∝ (λ/d i ) −0.69 which is much closer to the hybrid scaling. It is conceivable that the presence of this mirrorlike instability upstream of the reconnection site leads to further reductions in the reconnection rate in large systems, however, further study is needed to address this issue.

VI. CONCLUSIONS
Simulations of magnetic reconnection using the island coalescence problem set-up naturally incorporate several key features of real reconnecting systems that are absent in simpler elongated current layer models. In particular, the island coalescence problem includes the self-consistent formation of the current layer as the magnetic islands collide, and the two-way coupling between the micro-physics of the diffusion region and the macro-scale system that supplies the magnetic flux. For accurate modelling of near Earth space weather, where global fully kinetic simulations are currently impractical, it is important to understand the accuracy and limitations of reduced physics models for the simulation of such systems.
Refs. [20], [27], and [28] found significant differences between reconnection during island coalescence and earlier current layer studies in the zero guide-field limit. The Hall-MHD model, which had been thought sufficient to reproduce the reconnection rates of fully kinetic simulations, 9 was found inadequate to describe the kinetic results for the island coalescence problem. The Hall-MHD model overestimates the reconnection rate, is unable to reproduce its dependence on ion temperature or systemsize, and gives incorrect values for the pile-up magnetic field strength, outflow velocity, and ion diffusion region thickness. 27 The hybrid model, with massless fluid electrons and fully kinetic ions, was demonstrated as the minimum sufficient model to reproduce the above features of the island coalescence problem in the zero guide field limit. 27 The 10-moment fluid model was found in good agreement with kinetic results for the smallest system-size, but for larger systems is unstable to the formation of many secondary magnetic islands which are not present in the full kinetic simulations. 28 In this study we considered the sensitivity of these results to the addition of a finite guide field, since this can modify the kinetic ion physics that is crucial in setting the reconnection rate. We consider a range of guide fields B g /B m = 0 − 5.66 that is suitable for magnetospheric applications. Unexpectedly it is found that for a limited range of guide fields (B g /B m < 2.83) the reconnection rate, outflow velocity and pile-up field strength increase with increasing B g in the full kinetic, hybrid and 10-moment models. The same quantities for Hall-MHD vary little over this range of B g , and for B g /B m ≥ 2.83 the values of these quantities are in good agreement between the Hall-MHD and fully kinetic simulations. To understand this behaviour, the physics breaking the ion frozen-in condition and contributing to the outflow jet momentum balance was considered.
Only for the large guide fields B g /B m ≥ 2.83 is there qualitative agreement in the physics setting the ion diffusion region thickness, where ion inertia is the first term to break the frozen-in condition for ions in both the Hall-MHD and kinetic models. At the X-point, this guide field is sufficient to suppress meandering motion, so that the distribution function becomes to lowest order gyrotropic around the guide field direction. The outflow jet speed is reduced significantly with respect to the upstream ion Alfvén speed v A,in in the zero guide-field kinetic case due to a quasi-viscous effect associated with the non-gyrotropic part of the ion pressure tensor. For B g /B m = 2.83 there are still significant differences in the strength of the ion stress tensor and total pressure contributions to momentum balance between Hall-MHD and kinetic, but these terms mostly cancel so that they do not influence the outflow speed. The magnetic tension is balanced by ion inertia to give an outflow speed comparable to v A,in .
For the island coalescence problem, small changes of the peak reconnection rate can have larger consequences on the evolution of the system. If the reconnection rate is too slow, the islands can bounce off each other and reconnection can be temporarily switched off until the islands come together again. It is found that in large kinetic systems, islands with strong guide field have larger average reconnection rates and merge on a faster timescale than those with zero guide field. In addition to the changes in outflow velocity and ion diffusion region geometry with guide field that are present in small systems, large fully kinetic systems with guide field are also unstable to repeated formation of secondary magnetic islands which may play a role in limiting the length of the reconnection layer. Also, the presence of a mirror-like instability found in the largest zero guide field kinetic systems may further reduce the reconnection rate, but further work is needed to confirm this conclusion. This instability is not present in the strong guide field runs.
Reconnection events at the Earth magnetosphere can have a wide range of guide fields, but in the Earth's magnetotail it is often small 60 compared with the reconnecting field. Indeed, many of the early observations from the Magnetospheric Multiscale Mission have found reconnection events with weak guide field at the magnetopause, see e.g. Ref. [61]. This study, along with Refs. [27] and [28], demonstrate the limitations of the Hall-MHD model for reconnection events in such locations. The hybrid model, with massless fluid electrons and kinetic ions, and the 10-moment fluid model both give results in better agreement with the full kinetic system.