Particle-in-Cell Simulations of Collisionless Magnetic Reconnection with a Non-Uniform Guide Field

Results are presented of a first study of collisionless magnetic reconnection starting from a recently found exact nonlinear force-free Vlasov-Maxwell equilibrium. The initial state has a Harris sheet magnetic field profile in one direction and a non-uniform guide field in a second direction, resulting in a spatially constant magnetic field strength as well as a constant initial plasma density and plasma pressure. It is found that the reconnection process initially resembles guide field reconnection, but that a gradual transition to anti-parallel reconnection happens as the system evolves. The time evolution of a number of plasma parameters is investigated, and the results are compared with simulations starting from a Harris sheet equilibrium and a Harris sheet plus constant guide field equilibrium.


I. INTRODUCTION
Magnetic reconnection is one of the most fundamental plasma processes, and plays an important role in the magnetic activity of many astrophysical and laboratory plasmas [1,2]. It allows the conversion of stored magnetic energy into bulk flow, thermal and non-thermal energy, through changes in magnetic connectivity. In many astrophysical plasmas, the effects of particle collisions are negligible, and various aspects of collisionless reconnection have previously been studied in great detail . One particular aspect which has been investigated by a number of authors (e.g. Refs. 8, [10][11][12][13][14][15][16][17]19, and 21) is the influence of a guide field on the reconnection process. Most of these studies have used a Harris sheet [25] with a constant guide field as an initial current sheet configuration.
The addition of a constant guide field to the Harris sheet affects the evolution in a number of ways (see e.g. Ref. 2 for a more comprehensive overview than what we describe here). Some important points to note are as follows: (a) A constant guide field (of sufficient magnitude) has been shown to reduce the reconnection rate [8,12,13].
(b) The structure of the diffusion region is changed with the addition of a constant guide field [2]. In the anti-parallel (Harris sheet) case, the different outflow trajectories of the ions and electrons generate in-plane current loops (Hall currents), which in turn generate a quadrupolar out-of-plane magnetic field [4,9]. The addition of a constant (outof-plane) guide field results in a distortion of this quadrupolar field [8,10]. Furthermore, there is a * Electronic Mail: fw237@st-andrews.ac.uk † Electronic Mail: tn3@st-andrews.ac.uk strong parallel component to the out-of-plane electric field, which generates strong out-of-plane currents, and in-plane components of the parallel electron flows produce a density asymmetry along the separatrices.
(c) As a result of the density asymmetry described in point (b), in guide field reconnection there is a rotation of the reconnecting current sheet(s) [10,12,14,15,17,21].
(d) A guide field affects the particle orbits in the electron diffusion region [21] -it can destroy the bounce motion which occurs across the field reversal in the anti-parallel case, and so the length scales characterising the orbits in each case are different. In the guide field case, the relevant scale is the electron Larmor radius in the guide field, whereas in the anti-parallel case it is the electron bounce width in the reconnecting field component.
(e) A consequence of point (d) is that the addition of a guide field leads to thinner current sheets than in the anti-parallel case [17].
In this paper, we wish to address the following question: does the reconnection process change (and, if so, how?) if we use an initial current sheet configuration with a non-uniform guide field? We present results of a 2.5D particle-in-cell (PIC) simulation, in which we use an exact self-consistent equilibrium for the force-free Harris sheet as an initial condition [26,27].
Since the equilibrium guide field of the force-free Harris sheet (here B y = B 0 / cosh(z/L)) decreases with distance from the centre of the current sheet, we expect that the system will initially show features of guide field reconnection, but that a gradual transition to anti-parallel reconnection should take place, because plasma with smaller guide field strength should be transported towards the reconnection region as the system evolves in time. We will investigate whether and how this transition takes place, and also how it is reflected in the time evolution of arXiv:1512.07844v1 [astro-ph.SR] 24 Dec 2015 plasma quantities relevant for collisionless reconnection, such as the off-diagonal, non-gyrotropic elements of the electron pressure tensor.
Three-dimensional PIC simulations have previously been carried out for a magnetic field profile similar to that of the force-free Harris sheet [15], but with an additional constant guide field added in the same direction as the non-uniform guide field. The initial particle distribution functions were taken to be drifting Maxwellian distributions, which do not represent an exact initial equilibrium for this configuration. This leads us to discuss another motivation for our work -we are not aware of any previous study of collisionless reconnection for which exactly force-free initial conditions have been used for a nonlinear force-free field. The only known studies to use exactly force-free initial conditions have started from a linear force-free configuration [28][29][30][31][32]. Exact collisionless equilibria for such 1D linear force-free fields were first found approximately five decades ago [33,34], but the first exact equilibria of this type for nonlinear forcefree fields were found only very recently [26,27,[35][36][37][38]. Hence, only preliminary investigations have been carried out into the linear and nonlinear collisionless stability and dynamics of these configurations [39,40].
The structure of the paper is as follows. In Sect. II we discuss the simulation setup, followed by a detailed description of the results in Sect. III. We conclude with a summary and conclusions in Sect. IV.

A. Overview of Initial Configuration
For the main simulation run to be discussed, the initial magnetic field configuration is a force-free Harris sheet with added perturbation B p = B xpx + B zpẑ : where L is the current sheet half-width. The perturbation components have the form where L z is the half-width of the numerical box in the z-direction, a 0 = 0.1 and x m = L z /2. This gives an X-point reconnection site at the centre of the numerical box, and allows the nonlinear phase of the evolution to be studied without considering the cause of the reconnection onset. The x-component of the force-free Harris sheet magnetic field (when B p = 0) has the same spatial structure as that of the Harris sheet [25], and there is a nonuniform guide field in the y-direction, which is chosen in such a way that the total magnetic field strength is spatially uniform, and is given by B 2 0 = B 2 x + B 2 y . The resulting current density is parallel to the magnetic field, and hence the equilibrium is force-free [26,27]. A further consequence is that both the plasma density and P zz , the component of the pressure tensor that keeps the equilibrium in force balance, are spatially uniform. The equilibrium also has non-zero current density component in both the x-and y-directions, given by To initialise the particle positions and velocities in our main simulation run, we use the distribution function [26,27] where H s = (m s /2)(v 2 x + v 2 y + v 2 z ) is the particle energy, and p xs = m s v s + q s A x and p ys = m s v y + q s A y are the xand y-components of the canonical momentum (for mass m s , charge q s and vector potential components A x = 2B 0 L arctan(e z/L ) and A y = −B 0 L ln [cosh(z/L)]). The parameter β s is defined as β s = (k B T s ) −1 , where T s is the constant temperature of species s. Additionally, f 0s , a s , b s , u xs and u ys are constant parameters.
The distribution function (4) consists of a part which is equal to the Harris sheet distribution function [25], and an extra part which arises from the non-uniform guide field of the force-free Harris sheet. It should be noted that it can have a non-Maxwellian structure in velocity space. For further details of the properties of this function, see Refs. 26 and 27. To analyse the expected transition from guide field to anti-parallel reconnection in the force-free Harris sheet case, we will also present results from two other simulation runs: one which starts from a Harris sheet, and the other from a Harris sheet plus uniform guide field of B y = B 0 .

B. Normalisation and Parameters
To study the reconnection process, we use a 2.5D fully electromagnetic particle-in-cell code, which has been frequently used by Hesse and co-authors (see, for example, Refs. 5 and 9). The normalisation is as follows: the magnetic field is normalised to B 0 ; the number density to a free parameter, n 0 ; times to Ω −1 i = (eB 0 /m i ) −1 (the inverse of the ion cyclotron frequency in the equilibrium magnetic field); and lengths to the ion inertial length, c/ω i , where ω i = (e 2 n 0 / 0 m i ) 1/2 is the ion plasma frequency. Furthermore, velocities are normalised to the ion Alfvén velocity, v A = √ µ 0 m i n 0 , and so current densities and electric fields are normalised to B 0 /(µ 0 c/ω i ) and v A B 0 , respectively. In all simulation runs, we use an ion-electron mass ratio of m i /m e = 25. The total number of particles is 1.44 × 10 9 .The grid spacing in x and z is n x = 1200, n z = 600, and hence there are 2000 particles per cell. The numerical box has length L x = 64.0 and width L z = 32.0, which gives a grid spacing of ∆x = ∆z = 0.053. The boundary conditions are periodic at the x-boundaries, and specularly reflecting at the z-boundaries. The time step chosen is dt = 0.5/ω e (where ω e is the electron plasma frequency), with smaller time steps used occasionally. The ratio ω e /Ω e is set to equal 5. The ion-electron temperature ratio is equal to unity, with T i + T e = 0.5, so that T i = T e = 0.25. The current sheet half-thickness is equal to one ion inertial length: The various parameters from the force-free Harris sheet distribution function (4) have the following values:  Figure 1 shows the time evolution of the reconnected flux for the three simulation runs, and reconnection rates are shown in Figure 2. The maximum reconnection rate is highest in the Harris sheet case, occurring at t = 18. It has been observed in previous work that the effect of a constant guide field (of significant magnitude) is to reduce the maximum reconnection rate [8,12,13]. We see here that in the force-free run the maximum reconnection rate is further reduced from that of the constant guide field. It should be noted, however, that we used a parameter combination such that the initial electron number density is 25% higher in the force-free case than in the other two cases, which will have an effect on the reconnection rate. Figure 3 shows the y-component of the current density (in colour) and the projection of the magnetic field lines on to the x-z-plane, at various times for the force-free run. The figures show how reconnection leads to global changes in the structure of both quantities. At t = 0, it can be seen that the perturbation (2) to the magnetic field gives an initial X-point in the centre of the box. As time proceeds initially, a strong current sheet develops in the central region, and is slightly inclined, which is a typical feature of guide field reconnection [21]. As time proceeds beyond t = 20, the current sheet becomes more aligned with the x-axis, which could be a sign of a transition from guide field to anti-parallel reconnection.
Looking closely at Figure 3 for t = 20, it can be seen that a small magnetic island has started to form, which is a result of the bifurcation of the original X-point reconnection site into two new reconnection regions [9] one to the left of the island, and one to the right, which can be seen more clearly at later times. This is a feature commonly seen in reconnection simulations. Beyond t = 20, the island proceeds to move to the left, and eventually disappears, as the right-hand X-point begins to dominate over the left-hand one. By the end of the simulation, at t = 80, the island is no longer visible, and there is only one remaining reconnection region, which has shifted back towards the centre of the box. There is still a relatively strong current in this region though, which is higher than the original j y . Figure 4 shows the evolution of the non-uniform guide field in the force-free case. It can again be seen how the magnetic island starts to form around t = 20, and eventually disappears. At t = 40 and at subsequent times, a modified quadrupolar structure of B y can be seen around the X-point. This structure is qualitatively similar to that seen in Harris plus constant guide field simulations [8,10], and so we do not see a transition to the quadrupolar structure seen in Harris sheet simula- tions [4,8,9]. Figure 5 shows the variation of B y at the dominant X-point with time. It can be seen that, on the whole, there is a downward trend as time proceeds, representing a gradual transition from guide field to antiparallel reconnection (where B y would be close to zero). From around t = 35 onwards, B y fluctuates around a value of approximately 0.15. Of course, we do not have totally anti-parallel reconnection by the end of the simulation, but B y has clearly been significantly decreased from its initial value of 1.0 at the X-point. Figure 6 shows the x-component of the current density in the force-free case. The equilibrium j x is antisymmetric (see Eq (3)). As time proceeds, there is a build up of j x in the magnetic islands. These regions of strong j x correspond to regions where there is a strong gradient in the y-component of the magnetic field (see Figure 4). Similar behaviour has been seen in linear force-free simulations, and also in preliminary force-free Harris sheet simulations [39]. We have not included similar plots for the Harris and Harris plus constant guide field runs, but comment that j x is more prominent in the force-free case, which is to be expected since the other two cases have zero equilibrium j x .

B. Electron Larmor Radius and Bounce Width
In order to further investigate the expected transition from guide field to anti-parallel reconnection, we now consider the relevant length scales for the reconnection electric field. In the case of a guide field of significant magnitude, the electrons are strongly magnetised in the electron diffusion region, and r Ley = v th,e /(eB y /m e ), the thermal electron Larmor radius in the guide field B y , is the characteristic length scale [21]. As the guide field gets weaker, however, the important scale length is the electron bounce width in the reconnecting field component on the electron orbits is significant if When the condition (6) is satisfied at the reconnection site, therefore, we would expect to see mainly signatures of guide field reconnection and, when it is no longer satisfied, we would expect that this has coincided with a gradual transition towards anti-parallel reconnection, and would expect to see some signatures of this. In Figure 7, the ratio r Ley /λ z is plotted as a function of time. It first goes above unity between t = 24 and t = 25. We will consider t = 25 to be the 'transition time' towards anti-parallel reconnection, since after this time the ratio ceases to fluctuate around unity. Figure  8 shows the y-component of the current density and the magnetic field lines at this time (for the force-free run), together with plots at t = 18 and t = 16 for the Harris plus constant guide field and Harris runs, respectively (these are the times at which the reconnected flux in both cases matches that of the force-free case at t = 25). On the macroscopic level, the field-line structure looks more like that from the Harris sheet case, with an island separating two X-points. The central current sheet in the force-free case is still slightly inclined, but not as much as seen in Figure 3 at t = 20, and this inclination is also not as strong as seen in the constant guide field case.

C. The Reconnection Electric Field
In a 2D setup, the reconnection electric field is given by where the first bracket represents convection, the second represents the effect of the off-diagonal pressure tensor components, and the last bracket represents the effect of bulk inertia. Figures 9 and 10 show, for the force-free case, the contributions from each of the terms on the right-hand-side of Eq. (7) to the reconnection electric field, along x and z, through the average position of the dominant X-point, for data averaged between t = 24 and t = 26. The time we chose to average around is the 'transition time' discussed in Section III B, where the dominant scale for the evolution switches from the Larmor radius in the guide field B y , to the electron bounce width λ z . The pressure gradient terms are graphed as green lines, the convection term as purple lines, and the inertial term as black lines. The sum of these three terms, referred to as 'eytest', is plotted as a blue line in both plots. Although this fluctuates due to random noise, it can be seen that in both plots it matches reasonably well with the E y that is calculated on the numerical grid in the code (indicated by red lines).
From Figure 9, it can be seen that the pressure gradient term in x is significantly enhanced around the dominant X-point (x = 1.39). This increase in pressure coincides with a decrease (towards zero) of the convection term. Such behaviour can also be seen at x ≈ −3.75, which corresponds roughly to the position of the second, less dominant X-point (see Figure 8). In comparison with the other terms, the inertial term is small. The convection term should of course vanish at any X-points, since they are stagnation points where v s = 0, and so the pressure gradient term acts to support the reconnection electric field. This is in agreement with what has been found previously for Harris sheet and Harris sheet plus constant guide field simulations [3, 5-8, 12, 13, 16].
From Figure 10, it can be seen that, at z = 0 (the position of the dominant X-point), the convection term can be seen to drop to zero, and again the main contribution to E y comes from the pressure term. The inertial term is close virtually zero everywhere, apart from in the small region surrounding the X-point.

D. Pressure Tensor Components
We now focus on the structure of the off-diagonal components of the electron pressure tensor in the diffusion region, restricting attention to the electron quantities, since they are the dominant current carriers. Of particular importance are the non-gyrotropic components, which are given by [21] P e,ng = P e − P e,g , where is the gyrotropic component. The term (∇ · P e,g ) y vanishes at any X-points, since B x and B z vanish, and so non-gyrotropies of the pressure are required to give a contribution to the reconnection electric field [12]. Figures 11 to 13 show plots of the xy-and yzcomponents of the electron pressure tensor, together with the corresponding non-gyrotropic parts, at an early stage of the evolution, at which the total reconnected flux is the same in each case. The data has been averaged between t = 20 and t = 22 for the force-free run, t = 12 and t = 14 for the Harris run, and t = 13.87 and t = 15.87 for the constant guide field run. Note that we only show the non-gyrotropic components in the Harris case (Figure 12) because they are almost identical to the plots of the total P xye and P yze . The average location of the Xpoint under consideration is indicated by a green square. From Figure 11 for the force-free case, it can be seen that both P xye and P xye,ng have a gradient primarily in z, which is comparable to that from the constant guide field case in Figure 13. The structure of P yze and P yze,ng are also comparable to that from the constant guide field case in the vicinity of the X-point -these components all have gradients in x. Note, however, that P yze,ng in the constant guide field case also has a significant gradient in z, and so there is a significant difference in this component between the force-free and constant guide field cases. The structure of all pressure components in the vicinity of the X-point for the force-free case differ considerably from that of the Harris sheet case, which clearly has horizontal gradients in P xye and vertical gradients in P yze . It can be said, therefore, that in the early stages of the evolution the pressure in the force-free case exhibits (qualitatively) more features of guide field reconnection than anti-parallel reconnection.
As we discussed in Section III B, there is a change in the important scale length for the evolution around t = 25, from the Larmor radius in the guide field B y , to the electron bounce width λ z . Figure 14 shows nongyrotropic pressure plots for the force-free case, for data averaged around this transition time (between t = 24 and t = 26). On the whole, in the vicinity of the Xpoint, the structures remain qualitatively more similar to those from the constant guide field case (Figure 13), than the Harris case ( Figure 12).
To further investigate the transition, therefore, we now focus on a later time in the evolution. Figures 15 to 17 show the pressure components for data averaged between t = 45 and t = 47 for the force-free case, t = 25 and t = 27 for the constant guide field case, and t = 20 and t = 22 for the Harris case. As with the earlier Harris plot (Figure 12), we only show the non-gyrotropic components in Figure 16 because again they are almost identical to the plots of the total P xye and P yze . The structure of both P xye and P xye,ng in the force-free case is now significantly different than at earlier times. Focussing on the non-gyrotropic component, P xye,ng , the gradient is now primarily in the horizontal direction, and looks comparable (qualitatively) to P xye,ng for the Harris sheet. The other non-gyrotropic component, P yze,ng , now has significant gradients in both x and z, and still looks more similar to P yze in the guide field case than in the Harris sheet case. From Figures 15 to 17, we can conclude that some sort of transition has taken place in the structure of the pressure, since we see some signatures of anti-parallel reconnection. We can also conclude from this that the transition is not as simple as being from purely guide field reconnection to purely anti-parallel reconnection, but in- stead we see initially primarily signatures of guide field reconnection, and signatures of both guide field and antiparallel reconnection as the system evolves. This may be due to the fact that while B y at the dominant reconnection site (see Figure 5) decreases over time, it does not actually vanish completely, and Figure 4 shows that there is a modified quadrupolar structure of B y at later times -so not a transition to the quadrupolar structure seen in Harris sheet simulations. We speculate that this could cause some features of guide field reconnection to persist. This is clearly a point which should be investigated in future studies.

IV. SUMMARY AND CONCLUSIONS
In this paper, we have investigated how the reconnection process differs when adding a non-uniform guide field to the Harris sheet, instead of a constant guide field. We have presented results from a 2.5D fully electromagnetic particle-in-cell simulation of collisions magnetic reconnection, starting from a force-free Harris sheet with added perturbation and using the exact collisionless distribution function solution from Ref. 26 to initialise the particle velocities. For comparison, we have also presented results from a Harris sheet simulation, and a Harris sheet plus uniform guide field simulation.
We have found, as expected, that as time evolves in the force-free Harris sheet simulation, there are signs of a transition from guide field to anti-parallel reconnection. Firstly, on the macroscopic level, the initially rotated current sheet (similar to the constant guide field case) becomes more horizontally oriented (more like the Harris sheet case) as time progresses. Secondly, there is a gradual decrease in the guide field B y at the dominant X-point, indicating that it becomes less important as time proceeds. Thirdly, the transition can also be seen by looking at the ratio of the electron Larmor radius in the guide field B y , and the electron bounce width in the reconnecting field component, B x . The effect of the guide field on the electron orbits is significant if the ratio is less than unity [21]. At the beginning of the simulation, the ratio is well below unity, and begins to increase, eventually becoming greater than unity at a time of around t = 25. Finally, there are signs of a transition in the structure of the off-diagonal components of the electron pressure tensor. Initially in the force-free case, the structure and direction of the gradient in the vicinity of the X-point is more similar (qualitatively) to the constant guide field case, but at a later time in the evolution the structure looks more similar to the Harris case. It should be noted, however, that the transition we see is not as clear as going from purely guide field reconnection to purely anti-parallel reconnection, but instead we see initially primarily signatures of guide field reconnection, and signatures of both guide field and anti-parallel reconnection as the system evolves. This may be due to the fact that while B y at the dominant reconnection site decreases over time, it does not vanish completely, and there is a modified quadrupolar structure of B y at later times -not a transition to the quadrupolar structure seen in Harris sheet simulations. This could cause some features of guide field reconnection to persist, and is certainly a point open to further investigation.
The dominant contribution to the reconnection electric field, E y , was found to come from gradients of the offdiagonal components of the electron pressure tensor, in agreement wth previous findings for Harris and Harris plus constant guide field setups [3, 5-8, 12, 13, 16].
In this investigation, we have used only one set of parameters for the force-free run, which corresponds to a case where the ion distribution function is single-peaked in v y , and has a double maximum in the v x -direction, for small values of z around zero. The electron distribution function is single-peaked in both v x and v y . The distribution functions can of course both be single-peaked in v x for other sets of parameters, and can also have more pronounced double maxima in v x , as well as a double maximum in v y [27]. A future study could investigate how the evolution of the system depends on the initial velocity space profile for this equilibrium. The dependence of the evolution on other parameters could be investigated, such as mass ratio, temperature ratio or initial current sheet thickness.