2.1 CS formation

Consider a reversing, time-dependent magnetic field with a Harris (Reference Harris1962) profile without a guide field,

(2.1)\begin{equation} \boldsymbol{B}_{\rm r}(t,x) = {B}_{\rm r}(t)\tanh\left[\frac{x}{a(t)}\right]\hat{\boldsymbol{y}}, \end{equation}

which is frozen into an incompressible flow given by

(2.2)\begin{equation} \boldsymbol{u}(t,x,y) = \frac{1}{\varGamma(t) \tau_{\rm cs}} ({-}x\hat{\boldsymbol{x}} + y\hat{\boldsymbol{y}} ) \quad {\rm with} \ \varGamma(t) \doteq 1 + \frac{t}{\tau_{\rm cs}} . \end{equation}

This flow compresses the CS in the $x$ direction and stretches it in the $y$ direction over a characteristic CS formation timescale $\tau _{\rm cs}$, such that the CS maintains its profile but with a decreasing thickness $a(t)=a_0/\varGamma (t)$ and increasing field strength ${B}_{\rm r}(t)=B_{{\rm r}0}\varGamma (t)$ (Tolman, Loureiro & Uzdensky Reference Tolman, Loureiro and Uzdensky2018, § 2). (Here, and for the remainder of this paper, the ‘0’ subscript denotes an initial value.) The characteristic Alfvén Mach number of the flow (2.2) is

(2.3)\begin{equation} M_{{\rm A}0} \doteq \frac{a_0}{v_{{\rm A}0}\tau_{\rm cs}}, \end{equation}

where $v_{{\rm A}0}$ is the initial Alfvén speed of the reconnecting field. Although the magnetic-field profile $\boldsymbol {B}_{\rm r}(t,x)$ has no intrinsic lengthscale in the $y$ direction, any two points initially separated by a distance $L_0\hat {\boldsymbol {y}}$ will move apart and eventually be separated by $L_0\varGamma (t)\hat {\boldsymbol {y}}$, consistent with the incompressibility of the flow. Accordingly, we define a lengthscale $L(t)=L_0\varGamma (t)$ and refer to it as the ‘length’ of the CS; although the constant $L_0$ is somewhat arbitrary at this point, it will find a practical definition in our numerical simulations as the initial length of the computational domain in the $y$ direction (see § 3).

The thinning of this CS brings with it two important consequences. First, because the aspect ratio $L(t)/a(t)=(L_0/a_0)\varGamma ^{2}(t)$, the CS will become increasingly susceptible to tearing. Namely, the Harris-sheet tearing instability parameter

(2.4)\begin{equation} \varDelta'(t,N) = \frac{2N}{L(t)} \left[ \frac{1}{N^{2}} \frac{L^{2}(t)}{a^{2}(t)}-1\right] \end{equation}

will grow positively for all mode numbers $N\doteq k_{\rm t}(t)L(t)={{\rm const.}}$, where $\boldsymbol {k}_{\rm t}(t) = [k_{{\rm t}0}/\varGamma (t)]\hat {\boldsymbol {y}}$ is the tearing wavenumber made time-dependent due to its advection by the flow (2.2). Even for an initially tearing-stable CS satisfying $2{\rm \pi} a_0/L_0\geqslant 1$ (for which $\varDelta '(0,N)\leqslant 0$ for all $N$), the CS will eventually become unstable as it continually thins. Uzdensky & Loureiro (Reference Uzdensky and Loureiro2016) considered this time-dependent problem for the onset and evolution of the tearing instability of a thinning CS in the limit of resistive MHD. In § 2.2, we adapt their arguments for the case in which the resistive diffusion of the magnetic field is fourth order in space, viz. $\eta \nabla ^{2}\boldsymbol {B}\rightarrow \eta _4\nabla ^{4}\boldsymbol {B}$, relevant to the simulations described in § 3.

The second consequence of the thinning of the CS, at least for a sufficiently collisionless plasma, is the adiabatic production of pressure anisotropy in the inflowing fluid elements and, at sufficiently large $\beta$, the eventual triggering of the mirror instability. These events will of course affect any subsequent tearing of the CS, indeed, this is the entire point of this paper, but it is instructive to ask first how the CS would evolve without such interference.

2.2 Hyper-resistive tearing of a thinning CS without mirrors

Our first goal is to establish for how long a given CS must thin until tearing modes are able to onset and disrupt its formation. For that, we adapt the arguments given in Uzdensky & Loureiro (Reference Uzdensky and Loureiro2016) to the case with hyper-resistive tearing, relevant to the simulations described in § 3 that use fourth-order magnetic dissipation. The idea is as follows: because the instability parameter $\varDelta '(t,N)$ (see (2.4)) increases in time for each unstable tearing mode $N$, the growth rate $\gamma _{\rm t}(t,N)$ of each tearing mode increases as well, with the onset of tearing occurring only once $\gamma _{\rm t}(t,N)\tau _{\rm cs}\gtrsim 1$. A further complication in the case of tearing is that, as the CS lengthens, more and more modes (viz., larger values of $N$) become successively unstable, i.e. their $\varDelta '(t,N)$ becomes positive. One must then know how $\gamma _{\rm t}$ depends on $N$ to assess whether the tearing is ultimately dominated by a single island or by multiple islands: a question that depends upon whether $\varDelta '$ is small (${\ll }\delta ^{-1}_{\rm in}$, where $\delta _{\rm in}$ is the thickness of the resistive inner layer, corresponding to the ‘constant-$\psi$’ or ‘FKR’ regime; Furth, Killeen & Rosenbluth Reference Furth, Killeen and Rosenbluth1963) or large (${\sim }\delta ^{-1}_{\rm in}$, corresponding to the ‘Coppi’ regime; Coppi et al. Reference Coppi, Galvăo, Pellat, Rosenbluth and Rutherford1976).

For a Harris sheet subject to Ohmic resistivity, the tearing growth rate in the FKR regime satisfies $\gamma _{\rm t} \propto k^{-2/5}_{{\rm t}}$ for tearing wavenumber $k_{\rm t}$ satisfying $k_{\rm t}a\ll 1$. The fastest-growing FKR mode is then the longest one that fits into the CS (i.e. $N\sim 1$). With hyper-resistivity, however, the growth rate in the FKR regime is roughly independent of $k_{\rm t}$ for $k_{\rm t}a\ll 1$ (Huang, Bhattacharjee & Forbes Reference Huang, Bhattacharjee and Forbes2013):

(2.5a)\begin{align} \gamma _{\rm t} &\sim \frac{v_{\rm A}}{a} S^{{-}1/3}_a (\varDelta' a)^{2/3} (k_{\rm t}a)^{2/3} \end{align}

(2.5b)\begin{align} * \mbox{} &\sim \frac{v_{\rm A}}{a} S^{{-}1/3}_a \left( 1 - k^{2}_{\rm t} a^{2}\right)^{2/3} , \quad\text{where } S_a \doteq \frac{a^{3}v_{\rm A}}{\eta_4} \end{align}

is the hyper-resistive Lundquist number and $v_{\rm A}$ is the Alfvén speed associated with the reconnecting field ${B}_{\rm r}$. In the corresponding Coppi regime, Huang et al. (Reference Huang, Bhattacharjee and Forbes2013) find slower growth with $\gamma _{\rm t}\propto k^{4/5}_{\rm t}$ for $k_{\rm t}a\lesssim S^{-1/6}_a$ (see their figure 1). As a result, the fastest-growing unstable mode will always be an FKR-like mode with $\gamma _{\rm t} \sim (v_{\rm A}/a)S^{-1/3}_a$.

With all such modes growing at approximately the same rate, we argue that it will be the mode that first becomes unstable as the CS thins, viz. $N\sim 1$, that will ultimately win out (it having a head start over the others). Demanding that its growth rate satisfy $\gamma _{\rm t}\tau _{\rm cs}\gtrsim 1$ requires that $a \lesssim v_{\rm A} \tau _{\rm cs} S^{-1/3}_a \ll L$. Using our CS model with $L(t)/L_0=a_0/a(t)=v_{\rm A}(t)/v_{{\rm A}0}=1+t/\tau _{\rm cs}$ (see § 2.1) then specifies constraints on the critical time $t_{\rm cr}$ at which $\gamma _{\rm t}\tau _{\rm cs}\gtrsim 1$:Footnote ¹

(2.6)\begin{equation} \left(1+\frac{t_{\rm cr}}{\tau_{\rm cs}}\right)^{2/3} \frac{a_0}{L_0} \ll S_{a0}^{1/3} M_{{\rm A}0} \lesssim \left( 1 + \frac{t_{\rm cr}}{\tau_{\rm cs}} \right)^{8/3} . \end{equation}

With a view towards the results presented in § 4, we note that our numerical simulations generally have $S_{a0} \sim 10^{6}$ and $M^{-1}_{{\rm A}0} \in [1,16]$; the inequality (2.6) then implies values of $t_{\rm cr}/\tau _{\rm cs}$ that exceed unity by a factor of a few. In other words, for a thinning Harris CS without the production of pressure anisotropy and consequent excitation of mirror modes, we expect linear hyper-resistive tearing of our forming CS to onset at $N\sim 1$ after a few $\tau _{\rm cs}$. It is important to bear this point in mind for what follows.

2.3 Production and destabilisation of pressure anisotropy

We now demonstrate that the arguments given in § 2.2 for the onset of tearing in a thinning CS are likely to be irrelevant in a collisionless, high-beta plasma.

Consider a flux-frozen fluid element initially located at $x=\xi _0$ and moving towards $x=0$ in the velocity field (2.2). It is straightforward to show that, as the CS thins, the magnetic-field strength as seen by this element steadily increases as $B_{{\rm r}0}\tanh (\xi _0/a_0)\varGamma (t)$. If the first and second adiabatic invariants of the plasma are approximately conserved during this thinning, then the components of the pressure tensor parallel ($p_\parallel$) and perpendicular ($p_\perp$) to the local magnetic field become unequal, with the latter outpacing the former. Namely, if the plasma is initially pressure-isotropic, then the pressure anisotropy $\varDelta _p \doteq p_\perp /p_\parallel - 1$ in the fluid element should evolve according to

(2.7)\begin{equation} \varDelta_p(t,\xi(t)) = [\varGamma(t)]^{3} - 1 \approx \frac{3t}{\tau_{\rm cs}} , \end{equation}

the approximation being accurate for $t/\tau _{\rm cs}\ll 1$. Thus, pressure anisotropy increases in all fluid elements at the same rate.

In a plasma whose collision timescale is much longer than the CS formation timescale, the pressure anisotropy will continue to grow in time following (2.7) until one of two things occurs: either (i) the tearing instability onsets and disrupts the steady thinning of the CS or (ii) the pressure anisotropy grows large enough to surpass the mirror instability threshold,

(2.8)\begin{equation} \varLambda_{\rm m}\doteq\varDelta_p-\frac{1}{\beta_\perp}>0 , \end{equation}

where $\beta _\perp \doteq 8{\rm \pi} p_\perp / B^{2}$. With $\beta _\perp (t,\xi (t)) = \beta _\perp (0,\xi _0)/\varGamma (t)$, the latter scenario is possible in any inflowing fluid element so long as $\beta _\perp (0,\xi _0) > 4^{1/3}/3 \simeq 0.53$. Let us assume for the moment that case (ii) happens first, likely a good assumption when ${\beta _\perp (0,\xi _0)\gg 1}$, because in this situation the plasma becomes mirror-unstable (viz. $\varLambda _{\rm m}>0$) after a small fraction of the CS formation time. Namely, setting (2.7) equal to $1/\beta _\perp (t,\xi (t))$ and Taylor-expanding the result in $t/\tau _{\rm cs}\sim 1/\beta _\perp (0,\xi _0)\ll 1$, we find that the plasma becomes mirror-unstable at a time $t_{\rm m}$ that satisfies

(2.9)\begin{equation} \frac{t_{\rm m}}{\tau_{\rm cs}} \approx \frac{1}{3\beta_\perp(0,\xi_0)} \ll 1. \end{equation}

In the next two subsections, we use this estimate to predict the evolution of mirror fluctuations, including when they should regulate the pressure anisotropy and what the consequences are for tearing in a mirror-infested sheet.

2.4 Mirror instability in a forming CS

For a uniform plasma with $\varLambda _{\rm m}>0$ threaded by a static, uniform magnetic field, the maximum growth rate of the mirror instability is given by $\gamma _{\rm m}\sim \varOmega _i \varLambda ^{2}_{\rm m}$, where $\varOmega _i$ is the ion-Larmor frequency (Hellinger Reference Hellinger2007). In the asymptotic limit $\beta _\perp \varLambda _{\rm m}\ll 1$, the field-parallel and field-perpendicular wavenumbers at which this growth occurs satisfy $k_{\parallel,{\rm m}}\rho _i \sim (k_{\perp,{\rm m}}\rho _i)^{2}\sim \varLambda _{\rm m}$, where $\rho _i$ is the ion-Larmor radius. Two things complicate the application of these formulae to our forming CS. First, the reconnecting field (2.1) is non-uniform, with a null line occurring at the centre of the CS. This means that both $\varOmega ^{-1}_i$ and $\rho _i$ are smaller away from the neutral line, with ions unable to execute Larmor motion for $|x|\lesssim (\rho _i a)^{1/2}$ (Parker Reference Parker1957; Dobrowolny Reference Dobrowolny1968; Chen & Palmadesso Reference Chen and Palmadesso1984). As a result, mirror modes will grow faster and on smaller physical scales away from the CS ($|x|\gtrsim a$), where the plasma is better magnetised.

Secondly, the mirror growth rate and wavenumbers are made time-dependent by the increasingly positive pressure anisotropy, and mirror modes can only emerge in the magnetised regions of the CS if their growth rate is much larger than the rate of CS formation, i.e. $\gamma _{\rm m}(t) \tau _{\rm cs}\gg 1$. With $\gamma _{\rm m}(t) \sim \varOmega _i \varLambda ^{2}_{\rm m}(t)$, this means that the instability parameter $\varLambda _{\rm m}$ must reach a value ${\gg }(\varOmega _i \tau _{\rm cs})^{-1/2}$ before the mirrors can grow fast enough to outpace the more leisurely production of pressure anisotropy. Thereafter, $\varLambda _{\rm m}$ will decrease as the plasma converts its free energy into magnetic fluctuations. We therefore anticipate a maximal value of $\varLambda _{\rm m}$ given by

(2.10)\begin{equation} \varLambda_{{\rm m},{\rm max}} = C_{\rm m} (\varOmega_i\tau_{\rm cs})^{{-}1/2} , \end{equation}

where $C_{\rm m}\gg 1$ is some constant (determined by the numerical simulations in § 4 to be ${\approx }14$).Footnote ² By Taylor-expanding $\varLambda _{\rm m}(t)$ about $t=t_{\rm m}$, we find that (2.10) is attained after a time $t_{{\rm m},{\rm reg}}$ that satisfies

(2.11)\begin{equation} \frac{t_{{\rm m},{\rm reg}}-t_{\rm m}}{\tau_{\rm cs}} \approx \frac{C_{\rm m}}{3} (\varOmega_i\tau_{\rm cs})^{{-}1/2}, \end{equation}

beyond which the pressure anisotropy is regulated towards the mirror-instability threshold, $\varLambda _{\rm m}=0$. Due to the mirrors’ super-exponential growth, this regulation will occur very rapidly after $t_{{\rm m},{\rm reg}}$.

As $\varLambda _{\rm m}$ tends towards zero, the transfer of free energy from pressure anisotropy to mirror fluctuations drives the fluctuation amplitude $\delta B/B$ to a value of approximately $\varLambda ^{1/2}_{{\rm m},{\rm max}}$ (Kunz et al. Reference Kunz, Schekochihin and Stone2014a). However, as the CS continues to thin, $\varDelta _p$ continues to be driven positively, and the mirror fluctuations must grow in amplitude to maintain a marginally unstable plasma. For a linear-in-time drive, they do so secularly with $\delta B^{2} \propto t^{4/3}$ (Schekochihin et al. Reference Schekochihin, Cowley, Kulsrud, Rosin and Heinemann2008; Kunz et al. Reference Kunz, Schekochihin and Stone2014a; Rincon, Schekochihin & Cowley Reference Rincon, Schekochihin and Cowley2015), as an increasing fraction of large-pitch-angle particles become trapped in the deepening magnetic wells. If tearing does not intercede and disrupt their evolution, these mirrors would grow in amplitude all the way to $\delta B/B \sim 0.3$, independent of $\varLambda _{{\rm m},{\rm max}}$, after which the ions pitch-angle scatter off sharp ion-Larmor-scale bends in the mirroring field and saturate the instability. This scattering breaks adiabatic invariance and thereby maintains marginal stability by severing the link between $\varDelta _p$ and changes in $B$ (Kunz et al. Reference Kunz, Schekochihin and Stone2014a; Riquelme et al. Reference Riquelme, Quataert and Verscharen2015). In the context of our forming CS, reaching this saturated state would take a time of approximately $\tau _{\rm cs}$.

2.5 Hyper-resistive tearing of a mirror-infested CS

The emergence of mirrors in the forming CS has two effects on any subsequent tearing. First, by appreciably wrinkling the reconnecting field on kinetic scales, nonlinear mirrors cause $\varDelta '$ to depart from that corresponding to an undisturbed Harris sheet (2.4). Proposing a simple model for the magnetic profile of a mirror-infested CS, Alt & Kunz (Reference Alt and Kunz2019) calculated the resulting $\varDelta '(k_{\rm t})$ and showed that it becomes positive (and, thus, the CS becomes tearing-unstable) at wavenumbers $k_{\rm t}$ up to that corresponding to the inverse location of the innermost magnetic mirror (which effectively replaces $a$ as the characteristic thickness of the CS). Those authors then argued that the location of the innermost magnetic mirror is set by the perpendicular wavenumber of that mirror, using $k_{\perp,{\rm m}}\rho _i \sim \varLambda ^{1/2}_{{\rm m},{\rm max}}$ with $\varLambda _{{\rm m},{\rm max}}$ given by (2.10) and $\rho _i$ and $\varOmega _i$ taking on their local values. In the case of a very weak guide field, and accounting for the Lagrangian compression during CS formation, this corresponds to a distance $x_{\rm m}$ that satisfies

(2.12)\begin{equation} \frac{x_{\rm m}}{a(t)} \sim \left(\frac{\rho_{i0}}{a_0}\right)^{4/7} (\varOmega_{i0}\tau_{\rm cs})^{1/7} \end{equation}

(see their (4.11)). For $a_0/\rho _{i0}\gg (\varOmega _{i0}\tau _{\rm cs})^{1/4}$, this dramatically extends the range of tearing-unstable wavenumbers (from $k_{\rm t}\sim 1/a$ up to ${\sim }1/x_{\rm m}$). It also expedites the onset of reconnection by boosting tearing growth rates, since the effective $\varDelta ' \sim 1/x_{\rm m}$ is independent of $k_{\rm t}$ and so (2.5a) implies $\gamma _{\rm t} \propto k^{2/3}_{\rm t}$. Secondly, the perturbations to the reconnecting field caused by the mirror instability will nonlinearly seed tearing modes with $k_{\rm t}\gtrsim k_{\parallel,{\rm m}}$, giving them a head start over any traditional tearing modes of the thinning Harris sheet (if they are indeed unstable). As a result, mirror-stimulated tearing should onset at a wavenumber $k_{\rm t}$ satisfying $x^{-1}_{\rm m} \gtrsim k_{\rm t} \gtrsim k_{\parallel,{\rm m}}$.

Taken together, these two effects suggest that the onset of reconnection in an evolving CS, driven by mirror-stimulated tearing modes, likely occurs earlier and at smaller scales than it would have without the mirrors. More quantitatively, if pretearing mirrors were to wrinkle the CS and effectively reduce the CS thickness by a factor of approximately $(x_{\rm m}/a)$ (see (2.12)), then the right-hand side of (2.6) would acquire a multiplicative factor of approximately $(a/x_{\rm m})^{4/3} \sim (a_0/\rho _{i0})^{16/21}(\varOmega _{i0}\tau _{\rm cs})^{-4/21}\gg 1$, thereby reducing the critical time significantly.

The remainder of the paper is dedicated to testing these ideas using numerical simulations of CS formation and reconnection onset.

3.1 Hybrid-kinetic treatment of a thinning CS

For our numerical simulations, we adopt the hybrid-kinetic approximation, in which a non-relativistic, quasi-neutral, collisionless plasma is modelled by kinetic ions (mass $m_i$, charge $e$) that interact with a massless electron fluid via an electric field $\boldsymbol {E}$ given by the following generalised Ohm's law:

(3.1)\begin{equation} \boldsymbol{E} ={-}\frac{\boldsymbol{u}}{c}\times\boldsymbol{B} - \frac{T_e}{en} \boldsymbol{\nabla} n + (\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{B})\boldsymbol{\times}\frac{\boldsymbol{B}}{4{\rm \pi} en} . \end{equation}

Our notation is standard: $n$ and $\boldsymbol {u}$ are the ion number density and bulk-flow velocity, respectively; $T_e$ is the (assumed constant and isotropic) temperature of the electron fluid; $c$ is the speed of light; and the magnetic field $\boldsymbol {B}$ satisfies Faraday's law of induction,

(3.2)\begin{equation} \frac{\partial {\boldsymbol{B}}}{\partial {t}} ={-}c\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{E} . \end{equation}

The ions are treated using the PIC method, in which they are represented by finite-sized macro-particles whose positions $\boldsymbol {r}_p$ and velocities $\boldsymbol {v}_p$ are governed by the characteristic equations

(3.3)$$\begin{gather} \frac{{\rm d} {\boldsymbol{r}_p}}{{\rm d} {t}} = \boldsymbol{v}_p , \end{gather}$$

(3.4)$$\begin{gather}* \frac{{\rm d} {\boldsymbol{v}_p}}{{\rm d} {t}} = \frac{e}{m_i} \left[ \boldsymbol{E}(t,\boldsymbol{r}_p) + \frac{\boldsymbol{v}_p}{c}\boldsymbol{\times}\boldsymbol{B}(t,\boldsymbol{r}_p)\right] . \end{gather}$$

These equations are solved using the second-order-accurate code Pegasus++ (Arzamasskiy et al., in preparation); we refer the reader to Kunz, Stone & Bai (Reference Kunz, Stone and Bai2014b) for algorithmic details.

To model the thinning of the CS, we impose an externally driven, immutable flow that incompressibly expands the plasma along the sheet direction while contracting it in the perpendicular direction. To avoid having the box itself deform, we perform the simulation in the frame of the flow by applying a continuous coordinate transformation to the equations during each time step of the integration. To do so, we introduce the time-dependent Jacobian transformation matrix ${\boldsymbol{\mathsf{\Lambda}}}(t)\doteq \partial \boldsymbol {r}/\partial \boldsymbol {r}'$ tying the lab frame $(t,\boldsymbol {r},\boldsymbol {v})$ to the comoving frame $(t'=t,\boldsymbol {r}'={\boldsymbol{\mathsf{\Lambda}}}^{-1}\boldsymbol {r},\boldsymbol {v}'={\boldsymbol{\mathsf{\Lambda}}}^{-1}\boldsymbol {v})$, and specify its form via

(3.5)\begin{equation} {\boldsymbol{\mathsf{\Lambda}}}(t) = \varGamma^{{-}1}(t) \hat{\boldsymbol{x}}\hat{\boldsymbol{x}} + \varGamma(t) \hat{\boldsymbol{y}}\hat{\boldsymbol{y}},\quad{\rm with}\ \varGamma(t) \doteq 1 + \frac{t}{\tau_{\rm cs}} . \end{equation}

This transformation has the consequence that the strength of the reconnecting field increases linearly in time as the CS thickness $a(t)\propto 1/\varGamma (t)$, as in § 2.1. Note that $\lambda \doteq \det {\boldsymbol{\mathsf{\Lambda}}}=1$ at all times, thereby preserving the area of each cell and thus their density. In this comoving (primed) frame, (3.1)–(3.4) become, respectively,

(3.6)$$\begin{gather} \boldsymbol{E}' ={-} \frac{\boldsymbol{u}'}{c}\boldsymbol{\times}\boldsymbol{B}' - \frac{T_e}{en'}\boldsymbol{\nabla}' n' + (\boldsymbol{\nabla}'\boldsymbol{\times}\boldsymbol{B}')\boldsymbol{\times}\frac{{\boldsymbol{\mathsf{\Lambda}}}^{2}\boldsymbol{B}'}{4{\rm \pi} e n'\lambda}, \end{gather}$$

(3.7)$$\begin{gather}\frac{\partial {\boldsymbol{B}'}}{\partial {t'}} ={-}c\boldsymbol{\nabla}'\boldsymbol{\times}\boldsymbol{E}' - \eta_4 \nabla'^{4} \boldsymbol{B}', \end{gather}$$

(3.8)$$\begin{gather}{}\frac{{\rm d} {\boldsymbol{r}'_p}}{{\rm d} {t'}} = \boldsymbol{v}'_p , \end{gather}$$

(3.9)$$\begin{gather}{}\frac{{\rm d} {\boldsymbol{v}'_p}}{{\rm d} {t'}} = \frac{e}{m_i} {\boldsymbol{\mathsf{\Lambda}}}^{{-}2} \left[ \boldsymbol{E}'(t',\boldsymbol{r}'_p) + \frac{\boldsymbol{v}'_p}{c}\boldsymbol{\times}\boldsymbol{B}'(t',\boldsymbol{r}'_p) \right] - 2{\boldsymbol{\mathsf{\Lambda}}}^{{-}1} \frac{{\rm d} {{\boldsymbol{\mathsf{\Lambda}}}}}{{\rm d} {t'}}\, \boldsymbol{v}'_p , \end{gather}$$

where $\boldsymbol {E}'={\boldsymbol{\mathsf{\Lambda}}}\boldsymbol {E}$, $\boldsymbol {B}'=\lambda {\boldsymbol{\mathsf{\Lambda}}}^{-1}\boldsymbol {B}$, $n'\doteq \lambda n$, and $\boldsymbol {u}'\doteq {\boldsymbol{\mathsf{\Lambda}}}^{-1}\boldsymbol {u}$ (see Hellinger & Trávnıček Reference Hellinger and Trávnıček2005, Appendix A). To facilitate the reconnection of magnetic-field lines, we have appended to (3.7) a fourth-order hyper-resistive diffusion with constant coefficient $\eta _4$ (discussed further in the following). This is the set of equations solved by Pegasus++; the final (velocity-dependent) term in (3.9) is straightforwardly incorporated into the semi-implicit Boris algorithm for solving particle trajectories alongside the $\boldsymbol {v}'_p\boldsymbol {\times} \boldsymbol {B}'$ rotation. Quantities in the lab frame are easily obtained ex post facto.

Before proceeding any further, we pause here for a moment to explain our adoption of a hybrid-kinetic model with a fourth-order hyper-resistivity, rather than an alternative approach in which the electron kinetics play a role and the reconnection is facilitated by electron inertia and pressure-tensor effects. In some ways, our use of hybrid kinetics is of a purely pragmatic nature: our focus is on the macroscale evolution of a forming CS, the generation of pressure anisotropy in the inflowing fluid elements, and the effect of emergent mirror modes on the stability of the thinning sheet. All of this physics is captured within hybrid kinetics, which has the added bonus of being significantly cheaper numerically than a fully kinetic approach (especially at high $\beta$ with large-scale separations). Moreover, hyper-resistivity is often used to mimic the role of anomalous electron viscosity in facilitating reconnection, and yields tearing growth rates in the FKR regime that are approximately independent of $k_{\rm t}$ for a Harris sheet (recall (2.5b)), just as in the fully collisionless case (e.g. Drake & Lee Reference Drake and Lee1977; Karimabadi, Daughton & Quest Reference Karimabadi, Daughton and Quest2005; Fitzpatrick & Porcelli Reference Fitzpatrick and Porcelli2004, Reference Fitzpatrick and Porcelli2007).

Practical considerations aside, however, there is some astrophysical justification for resolving the ion-Larmor scale while breaking the frozen-in condition through a resistive term at sub-$\rho _i$ scales. In the hot and dilute ICM, $\rho _i \sim 1\ {\rm npc}$ is far below any meso- or macro-scale (including the Coulomb mean free path) but is only a factor of a few larger than the Ohmic-resistive scale (assuming Coulomb collisions; Schekochihin & Cowley Reference Schekochihin and Cowley2006). Such a scale hierarchy places the ICM (and our simulations) in the ‘semi-collisional’ regime specified by the ordering $d_e \ll \delta _{\rm in}\ll \rho _i \ll a$ where $d_e$ denotes the electron skin depth, in which MHD is no longer a sufficient description but the frozen-flux constraint is broken by resistivity rather than electron inertia (Cowley, Kulsrud & Hahm Reference Cowley, Kulsrud and Hahm1986; Bhat & Loureiro Reference Bhat and Loureiro2018, although their focus was on the low-beta limit). All this is to say that our hybrid-kinetic approach is not without physical utility or direct application to actual systems.

Of greater potential consequence is our neglect of electron pressure anisotropy, which could affect the reconnection dynamics in a number of ways. First, both in observations of reconnection in the Earth's magnetotail (e.g. Øieroset et al. Reference Øieroset, Lin, Phan, Larson and Bale2002) and in previous theoretical work on collisionless reconnection (Egedal et al. Reference Egedal, Le and Daughton2013), pressure anisotropy in the electron species (with $p_\parallel > p_\perp$) has been found to influence the inner diffusive layer of the CS. The idea is that adiabatic trapping of the electrons by magnetic mirrors and parallel electric fields within the reconnection region cause strong parallel heating of the electrons, producing pressure anisotropy that alters the reconnection geometry by driving electron currents in extended layers. Hybrid simulations that incorporate electron pressure anisotropy have shown the formation of more elongated CSs near X-points when compared with those found when adopting a pressure-isotropic closure, although the reconnection rate in both cases was found to be similar (Le et al. Reference Le, Daughton, Karimabadi and Egedal2016). Because we are considering the case with zero guide field, we do not anticipate such effects to play a role; without a guide field, the electrons are not adiabatically trapped in the inner layer and so can stream freely across the sheet and phase mix with other electrons (Le et al. Reference Le, Egedal, Ohia, Daughton, Karimabadi and Lukin2013). They can, however, develop significant agyrotropy in their pressure tensor as they partially demagnetise in the electron diffusion layer and bounce in the field-reversal region (Hesse, Birn & Kuznetsova Reference Hesse, Birn and Kuznetsova2001). In collisionless reconnection, this agyrotropy contributes to the reconnecting electric field (Vasyliunas Reference Vasyliunas1975; Hesse et al. Reference Hesse, Neukirch, Schindler, Kuznetsova and Zenitani2011) and can be used as a proxy for identifying separatrices and X-points (Scudder & Daughton Reference Scudder and Daughton2008; Swisdak Reference Swisdak2016). Another way in which pressure anisotropy can develop during reconnection is due to Fermi acceleration within contracting islands, which increases the parallel pressure and, at sufficiently high plasma $\beta$, becomes regulated by the firehose instability (Schoeffler et al. Reference Schoeffler, Drake and Swisdak2011). However, all of these effects – particle trapping within the reconnecting layer and Fermi acceleration during island growth and merging – feature in CSs within which tearing has already onset and gone nonlinear. In the context of this paper, the more relevant missing physics is electron pressure anisotropy generated during the CS formation itself, in a manner analogous to that discussed in § 2.3 for the ions. In this case, the electron pressure anisotropy could contribute to destabilising the forming CS at high $\beta$, either by adding to the total pressure anisotropy that factors into the mirror instability threshold (Basu & Coppi Reference Basu and Coppi1982; Hellinger Reference Hellinger2007; Hellinger & Štverák Reference Hellinger and Štverák2018) or by triggering its own kinetic instabilities (e.g. electron whistler instability; Kennel & Petschek Reference Kennel and Petschek1966; Gary & Wang Reference Gary and Wang1996; Riquelme, Quataert & Verscharen Reference Riquelme, Quataert and Verscharen2016). All of these effects could be taken into consideration in future work by adopting a fully kinetic approach.

3.2 Initial and boundary conditions

For all of our simulations, we construct a doubly periodic, two-dimensional (2D) domain with initial size $L_x\times L_y$, in which we initialise a magnetic field having a double Harris-sheet profile with no guide field,

(3.10)\begin{equation} \boldsymbol{B}_{\rm r}(t=0) = B_{{\rm r}0} \left[ \tanh\left(\frac{x - x_{{\rm cs}, 1}}{a_0}\right) - \tanh\left(\frac{x - x_{{\rm cs}, 2}}{a_0}\right) - 1 \right] \hat{\boldsymbol{y}} . \end{equation}

The locations of the two CSs are set to be $x_{{\rm cs}, 1}\doteq L_x/4$ and $x_{{\rm cs}, 2}\doteq 3L_x/4$. As we vary the initial CS thickness $a_0$ in our parameter study, we keep $L_x/a_0=48$, a value large enough to minimise interactions between the two sheets.

Using this numerical setup eliminates some complications that can arise from adopting the more common (and, arguably, more appropriate) open boundary conditions. In addition to the challenge of runaway particles in the open boundary setup, it is difficult to choose an appropriate particle distribution function in the upstream. On the other hand, by disallowing reconnected magnetic flux from leaving the domain, periodic boundaries lead to unphysical behaviour at long times. Fortunately, with our focus being solely on the effect of pressure anisotropy and mirrors on the onset of tearing modes, and not on the long-time evolution of reconnection, this should not pose any problem. Having two CSs in each simulation has the added benefit of reducing statistical noise when computing mean quantities.

We normalise the magnetic field to the initial strength of the reconnecting field far away from the CS, $B_{{\rm r}0}$, and all velocities to the initial Alfvén velocity, $v_{{\rm A}0} \doteq B_{{\rm r}0}/ \sqrt {4{\rm \pi} m_i \langle n_{i0} \rangle }$, where the initial ion density is averaged over the entire simulation domain and set equal to unity. The simulation time is normalised to the initial ion gyrofrequency, $\varOmega _{i0} \doteq e B_{{\rm r}0}/ m_i c$, and all lengthscales are normalised to the initial ion skin depth, $d_{i0} \doteq v_{{\rm A}0}/\varOmega _{i0}$. For our fiducial case with initial sheet width $a_0 = 125$, we choose $L_x \times L_y = 6000 \times 1500$ with $N_x \times N_y = 2688 \times 672$ cells. The hyper-resistivity is set so that $L^{3}_y v_{{\rm A}0}/\eta _4 = 5.625 \times 10^{8}$. The CS length is chosen such that the smallest available parallel wavenumber $k_{y,{\rm min}} \doteq 2{\rm \pi} / L_y$ is just barely unstable to the tearing instability discussed in § 2.2, that is, $k_{y,{\rm min}}a_0 \approx 0.5$. Additional simulations using a wider sheet with $a_0 = 250$ and $L_x=12\,000$, for which $\varDelta '(k_{y,{\rm min}}a_0)<0$, are also presented for comparison. We vary $\tau _{\rm cs} \in [1.25, 2.5, 5, 10, 20] \times 10^{2}$, which implies characteristic inflow velocities $a_0/\tau _{\rm cs}$ that are sub-Alfvénic for all but the fastest compression ($M_{{\rm A}0}\lesssim 1$).

Initially, we draw $N = 10^{4} N_x N_y$ ion macro-particles from a stationary Maxwell distribution corresponding to an initial ion temperature $T_{i0}$ and distribute them spatially to satisfy pressure balance within the double-Harris sheet. The (massless, fluid) electrons are set to be in thermal equilibrium with the ions, $T_e/T_{i0} = 1$. We initialise the ions to be pressure-isotropic with $\beta _{\parallel,i0} = \beta _{\perp,i0} = \beta _{i0}=100$; the initial ion-Larmor radius in code units is then given by $\rho _{i0} \doteq \beta _{i0}^{1/2}=10$. For these physical and numerical parameters, the ion-Larmor scale and the thickness of the diffusive inner layer are resolved at all times. Note, however, that neither the production of ion pressure anisotropy and its regulation by the mirror instability nor the evolution of the hyper-resistive tearing modes triggered by them within the forming CS require the resolution of the ion skin depth.

3.3 Watershed segmentation: identifying magnetic islands and X-points

Several useful diagnostics for quantifying the onset of reconnection require the identification and subsequent tracking of the X-points that form along the neutral line of the CS. In 2D reconnection, these points may be identified by tracing isocontours of the magnetic flux function $\psi$ or, equivalently, of the plane-perpendicular component of the magnetic vector potential. In this case, X-point locations are defined as the saddle points of $\psi$, where the first-order derivatives vanish and the second-order derivatives in two orthogonal directions have opposing signs. The latter can be quantified by computing the Hessian matrix of $\psi (x,y)$, defined as $H_\psi (x,y)$, and looking for cells with negative determinants (e.g. Servidio et al. Reference Servidio, Matthaeus, Shay, Cassak and Dmitruk2009). This saddle-point identification method relies heavily on the ability to perform accurate finite differences at the grid level; this may be difficult in data from PIC simulations due to the grid-scale noise coming from the (second-order-accurate) deposition of the ion density and momentum on the grid using a finite number of simulation particles. Typically, this ‘PIC noise’ masks the precise locations of zero crossings in the first-order-derivatives test (e.g. Haggerty et al. Reference Haggerty, Parashar, Matthaeus, Shay, Yang, Wan, Wu and Servidio2017).

To circumvent this issue, we have tried a different method discussed in Zhdankin et al. (Reference Zhdankin, Uzdensky, Perez and Boldyrev2013) that takes three differently sized loops around each cell and considers the cell to contain an X-point if the values of $\psi (x,y)$ rise above twice and fall below twice of the tested cell's value on each loop. Unfortunately, we have also found this procedure to be sensitive to fluctuations at the grid scale from PIC noise, resulting in spuriously identified X-points.

We have instead developed a novel method to locate X-points based on the watershed segmentation algorithm (Beucher & Meyer Reference Beucher and Meyer2018). This algorithm has been studied and utilised extensively to separate out regions of a scalar field (e.g. Mangan & Whitaker Reference Mangan and Whitaker1999). It treats the value of the field at any given point as though it were the topographic height of a relief and segments the region by its flood basins around local minima (hence, the name of the algorithm).

For X-point identification, $\psi (x,y)$ acts as the scalar variable of interest. Depending on the directional change of the reconnecting field, it might be necessary to consider $-\psi (x,y)$ instead, such that locations of the O-points on the neutral line can be the basins’ local minima. Following application of the algorithm, each magnetic island behaves similarly to the rainfall basin that surrounds an O-point. Preprocessing by applying a low-pass filter on $\psi (x,y)$ can further eliminate any erroneously segmented regions on sub-$\rho _i$ scales. We found that a Gaussian filter with a radius of size $\rho _{i0}$ works best to minimise the misidentification of basins. Compared with the other two methods discussed previously, this method is much more robust to the grid-scale noise coming from the PIC deposition scheme.Footnote ³

Given the watershed algorithm properties, it is possible that some of the segmented regions do not contain any magnetic island. These regions are typically the result of spuriously identified valleys located away from the neutral line and can easily be discarded. For two neighbouring regions that contain magnetic islands, an X-point should exist on their boundary and correspond to the local minimum value along that boundary. In addition, locations of O-points can be easily identified as the local minimum value of the flux function within the regions containing a magnetic island.

In figure 1 we present a direct comparison between these three algorithms for detecting X-point locations. We select a small region centred about the neutral line from our fiducial simulation (see figure 2c) that clearly exhibits three X-points. Figure 1($a$) demonstrates that the simple saddle-point method tends to under-count X-points, identifying only one out of the three X-point locations. The reason for this is that the constraints on the first-order derivative must be simultaneously satisfied for gradients in both the $x$ and $y$ directions. In contrast, the loop comparison algorithm shown in figure 1($b$) tends to over-count, or even mistakenly identify, X-points. Cells satisfying the test conditions can also be disconnected, complicating the interpretation of identified X-point locations. The result from our proposed watershed segmentation algorithm is shown in figure 1($c$), and demonstrates a robust detection, identifying all X-point locations.

Figure 1. Comparison of X-point detection methods for a given flux function $\psi (x,y)$ taken from one of our simulations: ($a$) simple saddle point, ($b$) loop comparison and ($c$) watershed segmentation. Values of $\psi (x,y)$ are represented as the greyscale shading and the interpolated blue contours. Cell locations determined to contain X-points are shaded in red, whereas the purple shaded cells in panel ($a$) show locations of other zero crossings of the first-order derivative. Yellow shaded cells in panel ($c$) show the boundary cells obtained from watershed segmentation.

Figure 2. Time slices centred about the CS at $x=x_{{\rm cs},1}$ from our fiducial simulation, showing properties of the magnetic field at ($a$) $t = 50$, ($b$) $t = 120$, ($c$) $t = 160$ and ($d$) $t = 250$. All quantities have been transformed back into the stationary lab frame; note the geometric thinning and lengthening of the CS. The size of an undisturbed Harris sheet of width $a(t)=a_0/\varGamma (t)$ is marked for reference by the green lines. The quantity $\delta B_y \doteq B_y - \varGamma (t)\tanh [x/a(t)]$ is represented by the colour contours, with the overlaid black lines tracing levels of the flux function $\psi$. In the bottom two panels, locations of the inferred X-points, found using the method discussed in § 3.3, are marked by the purple crosses and surrounded by finer levels of the flux function.