Particle-in-cell Simulations of Relativistic Magnetic Reconnection with Advanced Maxwell Solver Algorithms

High-energy radiation is observed from various astrophysical systems, such as pulsar wind nebulae, and from jets in active galactic nuclei (Giannios 2010), X-ray binaries (Tetarenko et al. 2017), and gamma-ray bursts (Piran 2004; Kumar & Zhang 2015). In particular, pulsars produce high-energy gamma-ray flares that evolve too rapidly to be explained by conventional particle acceleration theory (Abdo et al. 2011; Tavani et al. 2011). Magnetic reconnection is often invoked to explain the rapid nonthermal particle acceleration and emission in these systems (Lyutikov & Uzdensky 2003; Cerutti et al. 2012; McKinney & Uzdensky 2012; Nalewajko et al. 2015; Petropoulou et al. 2019; Philippov et al. 2019).

During magnetic reconnection, magnetic field energy is converted to particle kinetic energy in the form of both bulk motion and plasma heating. For the highly magnetized astrophysical systems, particle acceleration is caused by relativistic magnetic reconnection, where the unreconnected upstream plasma has a magnetic field energy density many times its enthalpy density. The plasma gains relativistic bulk and thermal velocities, and the particle energy distributions develop nonthermal high-energy power-law tails. This population of energized particles are thought to be a source of high-energy emission.

Particle-in-cell (PIC) is a well-established method for studying nonthermal plasma acceleration from first principles (Birdsall & Langdon 1991). Each species in the plasma is modeled with computational particles, which generate currents as they move in the domain. The current is deposited on a spatial grid and is a source term in the Maxwell equations. An electromagnetic field-solve step (also referred to as a Maxwell solve) calculates the electric and magnetic fields. The electromagnetic forces are interpolated to the particle positions, and their positions and velocities are advanced accordingly in time. The PIC method therefore fully captures kinetic particle acceleration, as well as the feedback of the accelerated plasma onto the fields.

A number of PIC studies have investigated relativistic reconnection for collisionless electron–positron (Zenitani & Hoshino 2001, 2007; Cerutti et al. 2012; Sironi & Spitkovsky 2014; Nalewajko et al. 2015; Werner et al. 2016) and electron-ion (Melzani et al. 2014) plasma in two-dimensions. Long thin current sheets become unstable to the tearing-mode instability, resulting in the formation and mergers of chains of trapped plasma, called plasmoids. Simulations have also indicated that reconnection progresses at a rate of approximately 0.1 in such systems, normalized to the reconnecting magnetic field and the Alfvén velocity (Guo et al. 2015; Cassak et al. 2017; Werner et al. 2018). Cerutti et al. (2014) investigated the dispersion relations of the tearing mode (in 2D) and the drift-kink mode that develops in 3D. The fastest-growing tearing mode in simulations has been shown to agree well with analytical expectations (Zenitani & Hoshino 2007). Recent work with PIC simulations has focused on the mechanisms underpinning the onset of reconnection and phases of particle energization (Guo et al. 2019; Hakobyan et al. 2021; Sironi 2022). The particle energy spectra that result from reconnection show hard power laws that extend to high energies (Guo et al. 2015; Werner et al. 2016, 2018; Petropoulou et al. 2019; Hakobyan et al. 2021). These spectra can then be combined with radiation models to predict observational signatures of reconnection in astrophysical systems (Cerutti et al. 2013; Nalewajko et al. 2018).

While substantial progress has been made in understanding particle acceleration in 2D, and more recently in 3D (e.g., Guo et al. 2015; Zhang et al. 2021; Schoeffler et al. 2023), algorithmic and computational innovation in PIC simulations of such systems has been limited. Virtually all studies employ a finite-difference time-domain Maxwell solver with a staggered Yee grid (Yee 1966; sometimes called FDTD), and only a few have explored the advantages of GPU acceleration for astrophysical PIC simulations (Bussmann et al. 2013; Chien et al. 2020; Xiong et al. 2023). The Yee approach is second-order in both space and time. For certain plasma systems, the numerical dispersion inherent to the method can lead to significant errors. With the Yee solver using a time step at the Courant limit, the numerical dispersion error is maximal along the axes and zero along the principal diagonals of the cells. To obtain solutions with less dispersion, we need alternate solvers, eventually based on higher-order methods. Cole and Kärkkäinnen proposed a nonstandard finite-difference approach to mitigate the effects of numerical dispersion along the principal axes when using the time step at the Courant limit (Cole 1997, 2002; Kärkkäinen et al. 2006), which Cowan et al. (2013) extended to noncubic cells. This combination is known as the Cole–Kärkkäinnen–Cowan (CKC) scheme. While numerical dispersion can be suppressed with CKC along the main axes, it remains at other angles.

Higher-order methods, including Fourier-based spectral methods, can be used to reduce dispersion even further. Pseudo-Spectral Analytical Time Domain (PSATD; Haber et al. 1973; Vay et al. 2013) is one such method, which enables arbitrary-order accuracy that can be set at runtime. Because they are finite-difference schemes, Yee and CKC are only numerically stable when the time step is below a value set by the Courant limit (see Section 2.3); on the other hand, PSATD is based on analytical integration in Fourier space and has no such constraint. In this paper, we compare the performance and accuracy of the two nonstandard approaches, CKC and PSTAD, for relativistic reconnection with the widely adopted Yee scheme. While PSATD does not have a Courant stability limit on the time step with regard to the Maxwell solve, an overly large time step may still reduce the accuracy of the simulation (as particles traveling close to the speed of light may travel over a cell size in a single time step). We therefore explore the performance and accuracy of PSATD with time steps above the light travel time across a cell.

The PIC algorithm captures reconnection physics accurately from first principles, but it can be computationally expensive, especially when performing high-resolution simulations with higher-order particle shape factors. Graphics processing units (GPUs) can offer remarkable acceleration over conventional CPU architectures for a number of scientific applications, including PIC (Bussmann et al. 2013; Germaschewski et al. 2016; Chien et al. 2020; Vay et al. 2020; Myers et al. 2021). We use the GPU-accelerated electromagnetic PIC code, WarpX (Vay et al. 2020; Myers et al. 2021). It has excellent full-machine scaling at leadership-class computing facilities, including Summit and Perlmutter (NVIDIA GPUs) and the world's first reported exascale machine, Frontier (AMD GPUs; Fedeli et al. 2022). The code is built on the AMReX (Zhang et al. 2019) framework, which supports MPI+X parallelism, where MPI enables inter-rank communications, and X corresponds to an interface such as OpenMP, CUDA, HIP, or SYCL for parallel programming on multicore CPUs or GPUs. PIConGPU (Huebl 2019), VPIC 2.0 (Bird et al. 2022), and the Plasma Simulation code (PSC; Germaschewski et al. 2016) also employ similar strategies that enable performance portability and allow scaling to multiple GPU nodes. Nonrelativistic magnetic reconnection has been used as a comparison case to validate multiple GPU-accelerated PIC codes, including PSC and sputniPIC (Chien et al. 2020), which can make use of a single node with multiple GPUs, and a CUDA Fortran single-GPU PIC code (Xiong et al. 2023).

In this paper, we perform first-of-their-kind 2D GPU simulations of relativistic reconnection with the advanced Maxwell solvers CKC and PSATD. Because these have never before been used for relativistic reconnection, we validate our results by comparing against simulations that use the conventional Yee solver, which has been well studied for 2D systems. In particular, we focus on the evolution of the current sheet structures, the particle-field energy balance, particle energy spectrum, and reconnection rate. We investigate the accuracy-based constraints on PSATD time steps by parameterizing the time step relative to the standard Courant limit of the finite-difference simulations. As with the simulations with different solvers, we compare our results with different time steps to the baseline Yee simulations. For the same cell size, the Courant limit for CKC admits a longer time step than Yee. PSATD is unconditionally stable with limitations on accuracy that may be imposed by other time-integration algorithms in the PIC simulation. Both CKC and PSATD may allow for faster simulations, so we compare the time to solution for these advanced solvers and assess the performance gains of a large time step with PSATD while holding the cell size constant.

The rest of the paper is organized as follows: in Section 2, we describe our simulations, including details about the initial configuration of the current sheets, the perturbation to trigger reconnection, and the algorithmic and numerical parameters. In Section 3, we discuss the accuracy and performance results from using different Maxwell solvers. In Section 4, we present the results from increasing the time step past the Courant limit with PSATD. We summarize and discuss future directions for our work in Section 5.

2.1. Harris-like Sheets

The simulations are initialized with two pair-plasma Harris-like current sheets (Harris 1962) in equilibrium. This section describes the configuration, which is summarized in the diagram in Figure 1. Full details of the derivation are in the Appendix, and summaries of important parameters can be found in Table 1 (scaled units) and Table 4 (SI units). Input files to replicate the simulations and analysis are available online. ³

Zoom In Zoom Out Reset image size

Table 1. Physical Parameters and Symbols Common to All of Our Simulations

Parameter	Symbol	Value
Background (cold) magnetization*	σ	30
Background temperature*	θb	0.15
Current sheet half-width*	δ	12.15 ρc
Current sheet skin depth	λe	2.45 ρc
Current sheet overdensity factor*	nd /nb	5
Current sheet velocity	β0	0.22 c
Current sheet temperature	θd	1.57
Domain half-width (x)*	Lx	2195 ρc
Domain half-width (z)*	Lz	1058 ρc

Note. Quantities marked with * are freely chosen; others are derived.

Download table as:  ASCIITypeset image

Far from the current sheets, in the upstream, the magnetic field is ${\boldsymbol{B}}=\pm {B}_{0}\hat{{\boldsymbol{z}}}$. Its sign changes at the current sheets. Our principal unit of time will be the inverse upstream electron gyrofrequency, ${\omega }_{{\rm{c}}}^{-1}\equiv {m}_{{\rm{e}}}/({{eB}}_{0})$, where c is the speed of light, e is the elementary charge, and m_e is the electron mass. Our base unit of length will be ${\rho }_{{\rm{c}}}=c{\omega }_{{\rm{c}}}^{-1}$, which is the nominal relativistic Larmor radius. The current sheets are chosen to have half-width δ = 12.15 ρ_c, and they extend in the y–z plane. They are centered at x = ±x_c ≡ ±L_x /2, where L_x is the half-width of the domain in the x-direction, which spans the interval [−L_x , L_x ].

We establish spatial distributions of number density n(x) and bulk velocity β (x):

$$\begin{eqnarray}&&n(x)={n}_{{\rm{b}}}+({n}_{d}-{n}_{{\rm{b}}})\left({\rm{sech}} \displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }+{\rm{sech}} \displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }\right)\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\boldsymbol{\beta }}(x)=\beta (x)\hat{{\boldsymbol{y}}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&=\ {\beta }_{0}\left({\rm{sech}} \displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }-{\rm{sech}} \displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }\right)\hat{{\boldsymbol{y}}}.\end{eqnarray} \tag{ 3 }$$

Each species has number density n(x). Positrons have velocity β (x), and electrons have velocity − β (x). The number density in the upstream region is n_b , and that in the current sheet is n_d . The bulk velocity at the center of the current sheet is β₀, and $\hat{{\boldsymbol{y}}}$ is the unit vector parallel to the y-axis. From Ampere's Law, ∇ × B = μ₀ J :

$$\begin{eqnarray}&&-\,\displaystyle \frac{\partial }{\partial x}{B}_{z}=2{\mu }_{0}{en}(x)\beta (x)c,\end{eqnarray} \tag{ 4 }$$

where B is the magnetic field, J is the current density, and μ₀ is the vacuum permeability. The factor of 2 in front of J_y is due to there being two species, positrons and electrons, that contribute to the total current density. Except where otherwise indicated, all quantities are given in the observer frame. We solve Equation (4) for B_z (x) with the boundary conditions

$$\begin{eqnarray}&&{B}_{z}(\infty )=-{B}_{z}(0)={B}_{0}.\end{eqnarray} \tag{ 5 }$$

This gives

$$\begin{eqnarray}\begin{array}{rcl}{B}_{z}(x) & = & -\displaystyle \frac{2{B}_{0}}{\left(\tfrac{\pi }{2}+\tfrac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)}\\ & & \times \,\left(\arctan \,\tanh \displaystyle \frac{x+{x}_{{\rm{c}}}}{2\delta }-\arctan \,\tanh \displaystyle \frac{x-{x}_{{\rm{c}}}}{2\delta }-\displaystyle \frac{\pi }{4}\right.\\ & & +\,\left.\displaystyle \frac{1}{2}\left[\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right]\left[\tanh \displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }-\tanh \displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }-1\right]\right),\end{array}\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&{B}_{0}=2{\mu }_{0}{{en}}_{{b}}{\beta }_{0}c\delta \left(\displaystyle \frac{\pi }{2}+\displaystyle \frac{{n}_{d}}{{n}_{{b}}}-1\right).\end{eqnarray} \tag{ 7 }$$

We choose the upstream "cold" magnetization $\sigma \,\equiv {B}_{0}^{2}/({\mu }_{0}{n}_{{b}}{m}_{{\rm{e}}}{c}^{2})=30$. This is somewhat different from the typical relativistic "hot" magnetization ${\sigma }_{{\rm{h}}}={B}_{0}^{2}/({\mu }_{0}h)$, with h the relativistic enthalpy density. The temperature in the upstream will be very mildly relativistic, with h ≲ 1.5n_b m_e c², so the hot magnetization is around 20, and σ ∼ σ_h. Based on either definition of magnetization, reconnection will proceed in the highly relativistic regime.

We choose the current sheet overdensity to be a factor of five, such that n_d = 5n_b . The skin depth in the current sheet is, by definition,

$$\begin{eqnarray}&&{\lambda }_{{\rm{e}}}\equiv \displaystyle \frac{c}{{\omega }_{{\rm{p}}}}=c\sqrt{\displaystyle \frac{{m}_{{\rm{e}}}{\epsilon }_{0}}{{n}_{d}{e}^{2}}},\end{eqnarray} \tag{ 8 }$$

where ω_p is the plasma frequency in the current sheet. From the quantities fixed thus far, λ_e = 2.45ρ_c. Combining Equation (7) with expressions for σ and λ_e gives an expression for the velocity at the center of the current sheet:

$$\begin{eqnarray}&&{\beta }_{0}=\displaystyle \frac{1}{\left(\tfrac{\pi }{2}+\tfrac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)}\displaystyle \frac{{\lambda }_{{\rm{e}}}\sqrt{\tfrac{{n}_{d}}{{n}_{{\rm{b}}}}\sigma }}{2\delta }=0.22\,c.\end{eqnarray} \tag{ 9 }$$

We calculate the temperature profile from pressure balance in the x-direction in the observer's (unprimed) frame:

$$\begin{eqnarray}&&{P}_{\mathrm{gas}}(x)+{P}_{\mathrm{mag}}(x)=C,\end{eqnarray} \tag{ 10 }$$

where C is a constant. The gas pressure P_gas = T^xx where T^{μ ν} is the stress-energy tensor. In the fluid (primed) frame, $T{{\prime} }^{{xx}}=2n^{\prime} \theta {{mc}}^{2}$, where θ = k_B T/(m_e c²) is the dimensionless temperature, and $n^{\prime}$ is the number density in the fluid frame. The fluid bulk velocity is in the y-direction, so $T{{\prime} }^{{xx}}={T}^{{xx}}$. The observer-frame number density is $n=\gamma n^{\prime}$, where $\gamma =1/\sqrt{1-{\beta }^{2}}$. Substituting,

$$\begin{eqnarray}&&{P}_{\mathrm{gas}}={T}^{{xx}}=2{{nmc}}^{2}\theta {\gamma }^{-1}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&=\,2{{nmc}}^{2}\theta \sqrt{1-{\beta }^{2}}.\end{eqnarray} \tag{ 12 }$$

The magnetic pressure,

$$\begin{eqnarray}&&{P}_{\mathrm{mag}}=\displaystyle \frac{{B}^{2}}{2{\mu }_{0}},\end{eqnarray} \tag{ 13 }$$

does not need to be transformed, given that B is already in the observer's frame.

We now have

$$\begin{eqnarray}&&2{{nmc}}^{2}\theta \sqrt{1-{\beta }^{2}}+\displaystyle \frac{{B}^{2}}{2{\mu }_{0}}=C.\end{eqnarray} \tag{ 14 }$$

Far from the current sheets,

$$\begin{eqnarray}&&n={n}_{{b}}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\theta ={\theta }_{{\rm{b}}}=\sigma /\eta \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\beta =0\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&B={B}_{0},\end{eqnarray} \tag{ 18 }$$

where we have chosen η ≡ 200, so θ_b = 0.15. Solving for θ(x):

$$\begin{eqnarray}&&\theta (x)=\displaystyle \frac{\sigma }{4}\displaystyle \frac{(4+\eta )/\eta -{[{B}_{z}(x)/{B}_{0}]}^{2}}{[n(x)/{n}_{{\rm{b}}}]\sqrt{1-\beta {\left(x\right)}^{2}}}.\end{eqnarray} \tag{ 19 }$$

The temperature at the center of the current sheet is calculated by evaluating Equation (19) at x = ±x_c, giving θ_d = 1.57.

Electron–positron pairs are initialized at the start of the simulation; they are arranged such that they are uniformly spaced, and their momenta are initialized by sampling from a Maxwell–Jüttner distribution at the local temperature θ (Zenitani 2015).

We note that the setup described here is slightly nonstandard. We have used spatially varying distributions of number density (Equation (1)), bulk velocity (Equation (3)), and temperature (Equation (19)) to represent both the upstream plasma and current sheets. More commonly, the current sheet plasma (with fixed hot temperature θ_d and drift velocity β₀) is overlaid on a domain filled with the upstream plasma (temperature θ_b, number density n_b ). Consequently, at initialization, the plasma at the center of the current sheets will have a two-temperature distribution. The number density of the drifting plasma is the only quantity that varies spatially. Ultimately, the resulting current distributions are very similar, as are the induced magnetic fields, and both configurations are in equilibrium. The primary difference is that the values of θ_d and β₀ used here differ from those of Werner et al. (2018), whose Harris sheets may otherwise appear identical to ours.

2.2. Perturbation

We add a perturbation to the equilibrium configuration in order to control the location and number of x-points. We model our perturbation after that of Werner et al. (2018), who add a one percent sinusoidal perturbation to the vector potential, which decreases the magnetic pressure at the z-axis just outside of each of the current sheets:

$$\begin{eqnarray}\begin{array}{rcl}{A}_{y} & = & \displaystyle \int {B}_{z}(x)\,{dx}\\ & & \times \,\left[1-0.01{\cos }^{51}\left(\displaystyle \frac{\pi }{{L}_{z}}z\right){\cos }^{2}\left(\displaystyle \frac{\pi (x-{x}_{c})}{{L}_{x}}\right)\right].\end{array}\end{eqnarray} \tag{ 20 }$$

We set the constant of integration equal to zero. In most studies that use this perturbation, the cosine in z is not raised to any power. We find that instead raising it to the 51st power narrows the region being perturbed and thus reduces the number of initial X-points in our configuration. This helps to ensure a more consistent comparison between the simulations. However, after the initial phases of reconnection, the form of the perturbation does not affect the results. The exact choice of power is arbitrary, though it must be odd in order to maintain the sign of the cosine term.

We initialize the fields in the simulation with ${\boldsymbol{B}}={\rm{\nabla }}\times ({A}_{y}\hat{{\boldsymbol{y}}})$. The simulations will therefore all start with ∇ · B = 0, which will be conserved by all numerical methods used here.

2.3. Particle-in-cell Simulations

3.1. Energy Conversion and Particle Acceleration

The PIC simulations of reconnecting Harris-like sheets (described in Section 2) were conducted until reconnection has completed and magnetic energy is no longer being converted to particle kinetic energy, at around t = 8000–$9000\,{\omega }_{{\rm{c}}}^{-1}\approx 4\times (2{L}_{z}/c)$. The qualitative current evolution of the Yee simulation is captured in Figure 2. The other simulations evolve similarly. Shortly after the start of the simulation, both the top and bottom current sheets (left and right columns in Figure 2) collapse due to the perturbation (Equation (20)) that lowers the magnetic pressure just above and below the current sheet (second row in Figure 2). Quasi-circular structures of trapped plasma and current form, called plasmoids. They inherit their average current from the current sheet where they form. Several magnetic X-points form between the plasmoids in the current sheets. The plasmoids move outward along the current sheet, and occasionally they merge, as highlighted in the blue box in the fourth row of Figure 2. This merger creates a new current sheet antiparallel to the bulk current in the plasmoids, and that extends perpendicularly from the original. This is the site of "secondary reconnection." At $t\approx 8000\,{\omega }_{{\rm{c}}}^{-1}$, both primary and secondary reconnection have ended, leaving a single plasmoid at z ≈ ± L_z (bottom row in Figure 2).

Zoom In Zoom Out Reset image size

In all three of the Maxwell solvers studied, magnetic reconnection proceeds at approximately the same rate and with the same structures. Figure 3 compares the plasma structures around the upper current sheet at $t=1470\,{\omega }_{{\rm{c}}}^{-1}$ for simulations that use the Yee, CKC, and PSATD solvers. Several small plasmoids have formed in each current sheet, with a single larger one forming at the edge of the domain. The exact positions and sizes of the smaller plasmoids differ between the solvers, but the current sheet fragments in a similar way in all cases.

Zoom In Zoom Out Reset image size

Energy conversion and particle acceleration also proceed similarly with all three solvers. When reconnection saturates at around $t=6000\,{\omega }_{{\rm{c}}}^{-1}$, about 40% of the energy in electromagnetic fields has been converted to particle kinetic energy (Figure 4). This includes energy associated with both thermal and bulk motion. Energy conversion proceeds nearly identically for the first $1800\,{\omega }_{{\rm{c}}}^{-1}$, after which there are slight differences between the numerical methods. This initial interval of identical evolution appears to be one of linear growth of the tearing-mode instability. The amplitudes of the fastest-growing spatial Fourier modes are approximately exponential in this interval. At around $t=1800\,{\omega }_{{\rm{c}}}^{-1}$, the exponential growth stops, suggesting the beginning of a nonlinear phase of evolution. This nonlinearity coincides with the appearance of small but noticeable differences between the energy conversions when using the different solvers, suggesting that this divergence is a consequence of nonlinear evolution. PSATD+Esirkepov shows a similar result.

Zoom In Zoom Out Reset image size

All simulations conserve energy to within one part in 2000, which is well within an acceptable level of energy nonconservation. The CKC simulation loses 4 × 10⁻⁴ of the initial total energy, slightly more than the 3.5 × 10⁻⁴ and 3 × 10⁻⁴ lost by the Yee and PSATD simulations, respectively.

The three methods and PSATD+Esirkepov also produce quantitatively similar particle acceleration. The highest particle γ at the start of the simulation is around 30, while by the end of reconnection at $t=8500\,{\omega }_{{\rm{c}}}^{-1}$, the fastest particles have γ = 500, over ten times higher. The majority of the energy in the domain is in particles with γ ≲ σ. Prior work has found that ${dN}/d\gamma \propto {\gamma }^{-\alpha }$ with α ≈ 1–3 (Guo et al. 2015; Werner et al. 2016, 2018). This corresponds to power laws with indices between 0 and −2 for the Lorentz factor distribution. Our energy distributions roughly follow this slope (Figure 5).

Zoom In Zoom Out Reset image size

3.2. Reconnection Rate

The dimensionless reconnection rate

$$\begin{eqnarray}&&\beta =-\displaystyle \frac{1}{{v}_{{\rm{A}}}{B}_{0}{L}_{x}}\displaystyle \frac{d{\rm{\Phi }}}{{dt}}\end{eqnarray} \tag{ 21 }$$

parameterizes much of the linear theory of reconnection (e.g., Werner & Uzdensky 2021). The Alfvén velocity is v_A, and Φ is the unreconnected flux. Directly measuring the amount of unreconnected flux is difficult, so we instead use the approximation β ≈ v_in/v_out, where v_in is the inflow velocity into the reconnection layer, and v_out is the terminal exhaust velocity downstream (Cassak et al. 2017; Liu et al. 2022). We calculate v_in by averaging ∣v_x ∣ within a region of size x_R × z_R = 245 ρ_c × 980 ρ_c = 0.11 L_x × 0.92 L_z centered on (x, z) = ( − x_c + x_R , 0). This region is marked by the green rectangle in Figure 2. The measured inflow velocity is relatively insensitive to the choice of x_R and z_R . The measurement is also symmetric across the current sheet, i.e., moving the box to x = − x_c − x_R does not change the measurement. The measurement is also similar on the other current sheet. For the sake of simplicity, we therefore only show reconnection rates on the +x side of the lower current sheet. The v_in average only includes cells in the upstream, i.e., where more than 70% of the plasma originated on the same side of the current sheet. Because reconnection mixes plasma across the current sheet, this excludes plasmoids and other reconnection exhaust. The value of this threshold does not strongly affect our results.

We measure the outflow velocity, v_out, by taking the median of the 10 highest cell-averaged z-velocities within δ = 12 ρ_c of the center of the current sheet. By $t=1000\,{\omega }_{{\rm{c}}}^{-1}$, this approaches the expected value of ${v}_{{\rm{A}}}=c\sigma /\sqrt{{\sigma }^{2}+1}\approx c$.

This estimates of the reconnection rate for Yee, PSATD, and CKC simulations are shown in Figure 6. The ratio v_in/v_out grows nearly identically in all simulations from 0 to 0.2 at a time of $1200\,{\omega }_{{\rm{c}}}^{-1}$. For the following $1200\,{\omega }_{{\rm{c}}}^{-1}$, the estimate of the rate holds relatively constant between 0.15 and 0.2 before dropping at around $t=2500\,{\omega }_{{\rm{c}}}^{-1}$. This coincides with 2L_z /c, the light-crossing time across the current sheet. After this point, the steady-state assumption under which β ≈ v_in/v_out breaks down as the reconnection fronts interfere with one another, due to the periodic boundary condition.

Zoom In Zoom Out Reset image size

The measurement of a rate between 0.15 and 0.2 matches the expectation from prior work. Typically, β for nonrelativistic reconnection is observed to be around 0.1, with higher rates for relativistic reconnection (Blackman & Field 1994; Guo et al. 2015; Comisso & Bhattacharjee 2016; Liu et al. 2017).

3.3. Performance

Current sheet structures (Figure 3), energy conversion (Figure 4), particle acceleration (Figure 5), and reconnection rate (Figure 6) are all similar between the Yee, PSATD, and CKC Maxwell solvers. This suggests that reconnection proceeds similarly in all of our simulations, demonstrating that the results produced by the PSATD and CKC solvers for the reconnection problem are comparable to those from the more conventional Yee solver. The PSATD+Esirkepov simulations also produce comparable results, despite the violation of Gauss's law. Adding this simulation to our performance comparison allows us to decouple the performance effects of the Maxwell solvers and current deposition schemes.

In these reconnection simulations, the majority of the runtime (∼60%) is spent in the current deposition routine, the performance of which is unaffected by the Maxwell solver. Consequently, the walltime per time step is approximately the same between the Yee, CKC, and PSATD+Esirkepov runs. CKC and PSATD permit a time step that is 40% longer than in Yee, reducing the time to solution for CKC and PSATD+Esirkepov by 40% and 30%, respectively, over the baseline Yee (see Table 2). PSATD+Esirkepov's steps take 10% longer than Yee(+Esirkepov) or CKC(+Esirkepov), indicating that the PSATD field-solve itself has only a small impact on the computational performance. However, PSATD(+Vay) has a 50% longer time per step, due to differences in the current deposition kernels. Thus, despite the longer time step permitted by the PSATD field solver, there is only a slight net increase in the time to solution. The implementation of Esirkepov deposition in WarpX is highly optimized, so it is possible that similar optimization in Vay deposition could make PSATD+Vay a more advantageous combination. This is an area for future work.

Table 2. Performance Comparison between Solvers

Algorithmic	Time Step	Walltime	Walltime to
Options	$({\omega }_{{\rm{c}}}^{-1})$	per Step (s)	$1470\,{\omega }_{{\rm{c}}}^{-1}$ (s)
Yee	0.411	0.077	274.6
CKC	0.581	0.077	193.5
PSATD (+ Vay)	0.581	0.115	290.0
PSATD+Esirkepov	0.581	0.083	209.9

Note. The walltime per time step is similar for Yee, CKC, and PSATD+Esirkepov, but 50% more expensive for PSATD(+Vay). Because PSATD and CKC allow longer time steps, CKC reaches a solution 40% faster, and PSATD(+Vay) is slightly slower than Yee. PSATD+Esirkepov performs comparably to CKC, reflecting that walltime per step is largely governed by the deposition scheme. These numbers are from simulations run to a final time of $t=1470\,{\omega }_{{\rm{c}}}^{-1}$ with dynamic load balancing and no I/O on 21 OLCF Summit nodes (126 GPUs).

Download table as:  ASCIITypeset image

A particular advantage of the PSATD method over either CKC or Yee is that it is not subject to a Courant stability criterion. Consequently, we are able to further increase the time step in the PSATD simulations by raising the CFL factor above 1. If other algorithmic choices are kept the same, then the time per time step is unlikely to change, reducing the time to solution. In this section, we study a sequence of simulations that are identical except for their CFL factors and therefore time steps. The CFL factors studied range from a baseline of 0.95–2.2. We refer to simulations in this sequence as CX.X where "X.X" is the CFL factor. C0.95 is the simulation labeled as "PSATD" in Section 3. As discussed in Section 2, we exclusively use the Vay current deposition scheme for the simulations in this sequence.

The statistical properties of the reconnecting plasma tend to remain the same early in the simulations and for lower values of the CFL factor. For values of the CFL factor ≲1.65, the energy conversion from magnetic fields to particles proceeds almost identically (Figure 7). We see the same behavior as with the baseline runs (Figure 4), where during the initial $\sim 1500\,{\omega }_{{\rm{c}}}^{-1}$, about 20% of the field energy is transferred to the particles. Following that phase, the energy conversion proceeds more slowly, saturating at a final distribution where about 40% of the energy is in particles and 60% remains in the fields.

Zoom In Zoom Out Reset image size

C1.7 matches the simulations with a lower CFL factor through the main period of interest, until reconnection saturates at around $7000\,{\omega }_{{\rm{c}}}^{-1}$. In the last $1000\,{\omega }_{{\rm{c}}}^{-1}$, it shows a slight increase in both electromagnetic field and particle energy, indicating that energy is not conserved. This is reflected in the relative energy conservation (see Figure 8). Nonconservation starts off small but increases sharply at around $t=7000\,{\omega }_{{\rm{c}}}^{-1}$, reaching 1% by the end of the simulation at $t=9000\,{\omega }_{{\rm{c}}}^{-1}$.

Zoom In Zoom Out Reset image size

The runs that all appear the same in the energy balance plot (CFL factor ≤1.65, Figure 7) show greater degrees of energy conservation. For runs with a CFL factor ≲1.6, nonconservation is comparable to the baseline, C0.95. C1.65 also shows low levels of nonconservation until the very end of the simulation, when it starts to increase. Left to run further, we expect that nonconservation would continue to increase like in the simulations with higher CFL factors. For the duration of the simulation, C1.65 is within a reasonable tolerance, only violating energy conservation by less than 1 part in 500.

For simulations with CFL factors ≥1.8, significant energy nonconservation develops during the last part of reconnection, much earlier than in C1.7. Errors in the energy have built up significantly by $4000\,{\omega }_{{\rm{c}}}^{-1}$ for C1.8, and by $2000\,{\omega }_{{\rm{c}}}^{-1}$ for both C2.0 and C2.2 (Figure 8). These are first apparent in the energy conservation plot, but continue to grow until they show in the plot of energy conversion (orange line in Figure 7). By $t=6000\,{\omega }_{{\rm{c}}}^{-1}$ the field and particle energies are visibly different from the baseline. C2.0 and C2.2 also show earlier runaway growth in field and particle energy, but they are omitted from Figure 7 for clarity.

The particle energy spectra further demonstrate agreement between the simulations with a CFL factor ≲1.7 throughout the majority of the simulation. Figure 9 compares the initial particle energy distributions (gray lines) with those at $t=2000\,{\omega }_{{\rm{c}}}^{-1}$ (mid-reconnection) and $t=8500\,{\omega }_{{\rm{c}}}^{-1}$ (after reconnection, almost at the end of the simulation). At the earlier time, all of the runs we study show close agreement with the baseline C0.95 out to γ ∼ 100. Those runs with a CFL factor less than 2 also agree closely with one another at the highest energies; there are few, if any, particles with Lorentz factors in excess of 300. C2.0 shows a slight overabundance of particles with γ ∼ 100, which grows as the simulation progresses. C2.2 shows an even larger overabundance of high-energy particles, reaching a maximum Lorentz factor of 400, 30% higher than the maximum achieved by C1.8 at this time. This also coincides with the beginning of the energy nonconservation seen in Figure 8. The overabundance of high-energy particles may directly cause the initial nonconservation, but shortly thereafter we also see artificial heating in the cold upstream. In C0.95 and the runs with lower (≤1.7) CFL factors, the cold upstream plasma appears as a peak at γ ∼ 1, and it remains largely unchanged even at the end of the simulation. However, in C1.8 and above, this peak broadens and moves out to γ ∼ 4. This is apparent near the end of C1.8 (bottom panel, Figure 9). This also occurs in C2.0 and C2.2, beginning shortly after the time of the top panel; we omit their distributions in the later panel for the sake of clarity.

Zoom In Zoom Out Reset image size

At the end of the simulations, the runs with CFL factors ≲1.7 have particle energy distributions that agree almost exactly for γ < 100. For the highest-energy particles, there are minor differences, comparable to the spread seen in the different electromagnetic solvers (see Figure 5). These highest-energy particles are also the rarest—there are ∼5000 times fewer particles at γ = 2 as there are at γ = 100, so the higher energies are more subject to statistical noise.

The nonconservation seen for CFL factors ≥1.65 in Figure 8 appears to be at least partially driven by numerical heating in the upstream region. The exact cause is uncertain, though we expect it is a direct consequence of the large CFL factor, rather than the under-resolution of a physical timescale. Coincidentally, C1.65 has a time step just over $1\,{\omega }_{{\rm{c}}}^{-1}$, meaning that the cyclotron motion of upstream particles cannot be captured. However, this heating does not occur if we obtain that same time step by, for example, doubling the cell size in each direction and reducing the CFL factor below 1. The effect is driven by the upstream plasma. Simulations of a domain with no current sheet, filled only with upstream (θ = 0.15, σ = 30) plasma, show the same numerical heating.

The dimensionless reconnection rate is remarkably similar for all runs except for C2.2 (Figure 10). As discussed in Section 3.2, we are primarily interested in the value of the estimated reconnection rate for $t\lt 2000\,{\omega }_{{\rm{c}}}^{-1}$. Within that time interval, the simulations all show a nearly identical rise to 0.2, followed by a slight decline to 0.15 over the following $1000\,{\omega }_{{\rm{c}}}^{-1}$. C2.2 largely follows this pattern, but declines much more quickly than the other runs at $t=2000\,{\omega }_{{\rm{c}}}^{-1}$. This highlights the robustness of the reconnection rate and indicates that it likely should not be used to diagnose whether the simulations are producing correct results. At $t=2000\,{\omega }_{{\rm{c}}}^{-1}$, C2.0 and C2.2 both show qualitative disagreement with the baseline models in the high-energy particle spectra (Figure 9). However, they agree with baseline models in the reconnection rate measurement.

Zoom In Zoom Out Reset image size

Of the diagnostics discussed here, energy conservation is the most sensitive to numerical problems that arise due to a high CFL factor. Because our simulations are closed systems, the total energy should remain the same throughout the evolution, providing a straightforward ground truth comparison. When increasing the CFL factor, we find that the particle acceleration and reconnection rate obtained from those simulations do not deviate from the baseline until after energy nonconservation begins to rise rapidly. Even C2.2 agrees with C0.95 for the first $1500\,{\omega }_{{\rm{c}}}^{-1}$ of the simulation, before its nonconservation starts to increase. While particle energy spectra may match the baseline even after energy nonconservation begins to run away (e.g., C1.7 at $t=8500\,{\omega }_{{\rm{c}}}^{-1}$, bottom panel of Figure 9), we would not be able to verify that particle acceleration was still correct in the absence of a known baseline. When increasing the CFL factor, we suggest using a runaway in energy nonconservation as a heuristic for closed systems to determine when a simulation result is unreliable.

In Table 3, we summarize the computational performance of simulations with CFL factors greater than 1. Again, we compare against the baseline C0.95 with the PSATD Maxwell solver. We find that, between C0.95 and C2.2, the walltime to simulate a fixed physical time interval decreases by half. Over this range of CFL factors, the physical time step increases by a factor of 2.3, and the walltime per time step increases by about 15%. Most of this walltime difference is because current deposition is slower (on a per step basis) at higher CFL factors. When particles move farther per step, they are more likely to move between cells, and this likely reduces the cache performance of the deposition routines.

Table 3. Performance of Simulations with CFL Factors between 0.95 and 2.2

CFL	Time Step	Walltime	Walltime to
$({\omega }_{{\rm{c}}}^{-1})$	Per Step (s)	$1470\,{\omega }_{{\rm{c}}}^{-1}$ (s)
0.95	0.582	0.115	290.0
1.5	0.919	0.126	200.9
1.6	0.978	0.127	190.9
1.65	1.010	0.128	186.4
1.7	1.041	0.128	180.9
1.8	1.102	0.129	172.5
2.0	1.225	0.131	157.2
2.2	1.347	0.131	143.0

Note. These numbers are from otherwise-identical simulations run to a final time of $t=1470\,{\omega }_{{\rm{c}}}^{-1}$ with dynamic load balancing and no I/O on 21 OLCF Summit nodes (126 GPUs).

Download table as:  ASCIITypeset image

The speedup obtained for a particular application is limited by how high one can increase the CFL factor while maintaining a reliable result throughout the time interval of interest. We expect that the interval over which energy is conserved for a given CFL factor will vary based on the problem setup. For the configuration described here, a CFL factor of 2.2 may be adequate if we were only interested in the onset of reconnection. In that case, the time to solution will be 50% of what it is in C0.95. In this study, we were interested in following reconnection from its onset until it saturates. For that, we needed a CFL factor of at most 1.7, and would have terminated the simulation at around $6500\,{\omega }_{{\rm{c}}}^{-1}$. At that CFL factor, we reach the solution 1.6 times faster than in the baseline C0.95.

We have performed first-of-their-kind particle-in-cell simulations of relativistic reconnection with a 16th-order pseudo-spectral Maxwell solver (PSATD) and a time step that exceeds the conventional CFL limit. We find that PSATD and the nonstandard finite-difference scheme CKC qualitatively and quantitatively produce the same results as the standard second-order finite-difference scheme Yee. We have verified that all three schemes produce the same qualitative plasmoid evolution, particle-field energy balance, particle acceleration, and reconnection rate (Figures 3–6). The particle energy distribution has a power-law tail with a spectral slope α of ${dN}(\gamma )/d\gamma \propto {\gamma }^{-\alpha }$. We measure α between 1 and 3, as expected (Guo et al. 2015; Werner et al. 2016, 2018). The reconnection rate is between 0.15 and 0.2 for most of the linear phase of the simulation, consistent with expectations for relativistic reconnection.

We also compare the performance of the solvers and measure efficiency in terms of time to solution. For the same CFL factor, CKC and PSATD allow for a time step that is longer than that of Yee by a factor of $\sqrt{2}$. The walltime taken per step, though, depends more on the current deposition scheme than on the Maxwell solver. The walltime per step is the same in CKC as in Yee (which both use Esirkepov deposition), giving a ∼40% speedup in time to solution. High-order PSATD is only charge-conserving when used with Vay deposition, which in its current WarpX implementation is not as optimized as Esirkepov and is therefore more computationally expensive. Thus, we decouple the comparison of the solvers and current deposition methods by performing an additional PSATD+Esirkepov simulation. In doing so, we verify that PSATD itself is not the primary cause of the more expensive time step. Despite not being charge-conserving, PSATD+Esirkepov gives the correct answer for this problem, and has a walltime per step only 10% higher than CKC or Yee. In conjunction with the $\sqrt{2}$ larger time step allowed by PSATD, this produces a net 30% reduction in time to solution. A time step in the charge-conserving PSATD+Vay takes 50% longer than in Yee, so it has a slightly longer time to solution than the most commonly used charge-conserving Yee algorithm.

Unlike either of the finite-difference schemes, PSATD is numerically stable at any time step, even one greater than the light travel time across a cell. We explore the accuracy and performance of CFL factors >1 in the relativistic reconnection problem. We find that the timescale of interest sets the maximum allowable time step parameterized by the CFL factor. Factors ≤1.65 conserve energy comparably well to the baseline PSATD case until $t=8000\,{\omega }_{{\rm{c}}}^{-1}$, well past the end of reconnection. In that interval, they show good agreement in both particle acceleration (Figure 5) and reconnection rate (Figure 10). A slightly higher CFL factor of 1.7 still shows good agreement in energy distribution (Figure 7) and particle acceleration, but it begins to show signs of runaway errors in energy at $t\approx 6500\,{\omega }_{{\rm{c}}}^{-1}$, right after reconnection saturates. This suggests that, for this particular problem and computational setup, the CFL factor could reach 1.7 without compromising the accuracy of the simulation results during reconnection. As we progressively increase the time step (up to a CFL factor of 2.2), runaway energy nonconservation begins earlier and earlier. In the interval where energy is conserved, the other diagnostics such as particle energization and reconnection rate agree with the baseline case, suggesting that energy nonconservation is a good diagnostic of when the simulation results are reliable for closed systems with periodic boundary conditions.

The walltime per time step only increases modestly with CFL factor, about 15% from the baseline C0.95 to C2.2. The increases in the time step outweigh the increases in walltime per step, reducing the walltime per physical time by a factor of about two between C0.95 and C2.2, if all other algorithmic choices are unchanged. Comparing against CKC+Esirkepov, the most efficient of the CFL = 0.95 simulations, we obtain a 25% reduction in time to solution by adopting a CFL factor of 2.2, and a <10% reduction by using the CFL factor of 1.7 that we have determined to be suitable for our scientific purposes. If a CFL factor below 1.7 were necessary, CKC+Esirkepov would be the most efficient option, due to the faster current deposition. Future work will include further optimization of the Vay deposition routine, which would shift the trade-offs to favor the PSATD+Vay combination in more situations. While PSATD+Esirkepov was accurate in these simulations, we do not recommend its use without an additional current correction or divergence-cleaning operation to guarantee charge conservation.

One of PSATD's strengths is that it reduces numerical dispersion that may appear when using Yee or CKC solvers. While an exploration into the mitigating effect of the PSATD solver on numerical dispersion in relativistic plasmas is outside the scope of this paper, it is indeed a known effect that can contaminate simulations of astrophysical jets, shocks, and magnetic reconnection (Godfrey 1974; Melzani et al. 2013; Godfrey & Vay 2014; Ikeya & Matsumoto 2015; Li et al. 2017; Nishikawa et al. 2021; Tomita et al. 2022). PSATD can therefore provide new opportunities to study highly relativistic plasmas, a numerically challenging regime that characterizes a number of astrophysical systems.

The numerical and algorithmic innovations leveraged in this study can be used to enable larger and more efficient simulations of astrophysical and laboratory plasmas (Ji et al. 2022). Our simulations are also some of the first GPU-accelerated astrophysical PIC simulations. As a first step, we have only focused on two-dimensional systems without additional kinetic and radiative physics. A third spatial dimension is dynamically important in reconnection because it makes the current sheet susceptible to an additional instability, called "drift-kink," which can suppress particle acceleration (e.g., Cerutti et al. 2014; Sironi & Spitkovsky 2014; Guo et al. 2015; Werner & Uzdensky 2017). CKC and PSATD may be especially efficient in 3D because their time steps are larger than Yee's by a factor of $\sqrt{3}$, rather than $\sqrt{2}$ in 2D. In many astrophysical reconnection environments, synchrotron emission and pair production play an important role. With WarpX, we will be able to take advantage of GPU-accelerated exascale computing resources to perform 3D simulations that include this radiative physics.

We thank Luca Comisso for the insightful discussions on magnetic reconnection and reconnection rate. This work was partially supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Exascale Computing Project (17-SC-20-SC) under contract DE-AC02-05CH11231, and Simulation and Analysis of Reacting Flow, FWP #FP00011940, funded by the Applied Mathematics Program in ASCR. This research used the open-source particle-in-cell code WarpX https://github.com/ECP-WarpX/WarpX, primarily funded by the US DOE Exascale Computing Project. The primary WarpX contributors are affiliated with LBNL, LLNL, CEA-LIDYL, SLAC, DESY, CERN, and Modern Electron. We acknowledge all WarpX contributors. This work used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. This research also used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231 using NERSC awards DDR-ERCAP m3881 for 2022 and ASCR-ERCAP mp111 for 2022 and 2023.

Software: WarpX (Fedeli et al. 2022; Vay et al. 2023), AMReX (Zhang et al. 2019; the AMReX Development Team et al. 2023), matplotlib (Hunter 2007), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), yt (Turk et al. 2011).

A.1. Equilibrium

In all of our simulations, we fix ρ_c = 4.1 × 10⁻³ m, which gives B = 0.42 T and fixes the other dimensional values in Table 4. These are used in the WarpX simulations, which require SI units. This choice of length scale is largely arbitrary. In this nonradiative regime, the magnetization determines most of the physics; the rest of the results will scale accordingly.

Table 4. Physical Parameters and Symbols Common to All of Our Simulations

Parameter	Symbol	Value
Background Larmor radius	ρc	4.1 × 10−3 m
Background Larmor frequency	ωc	7.3 × 10−10 s−1
Skin depth	λe	0.01 m
Current sheet half-width	δ	0.05 m
Background magnetization	σ	30
Background magnetic field	B0	0.42 T
Current sheet number density	nd	2.8 × 1017 m−3
Background number density	nb	5.6 × 1016 m−3
Current sheet velocity	β0	0.22 c
Background temperature	θb	0.15
Current sheet temperature	θd	1.57
Domain half-width (x)	Lx	8.96 m
Domain half-width (z)	Lz	4.32 m

Download table as:  ASCIITypeset image

From the relationship between ρ_c and λ_e and the definition of skin depth, we can calculate number density in the current sheet n_d :

$$\begin{eqnarray}&&{n}_{d}=\displaystyle \frac{m{\epsilon }_{0}{c}^{2}}{{e}^{2}{\lambda }_{{\rm{e}}}^{2}}=2.8\times {10}^{17}\,{{\rm{m}}}^{-3},\end{eqnarray} \tag{ A1 }$$

and in turn, the background number density n_b = 5.6 × 10¹⁶ m⁻³. At this density, a magnetization of 30 requires the background magnetic field B₀ = 0.42 T.

After integrating Ampere's Law (Equation (4)),

$$\begin{eqnarray}\begin{array}{rcl}{B}_{z}(x) & = & -4{\mu }_{0}{{en}}_{{\rm{b}}}c{\beta }_{0}\delta \left[\arctan \left(\tanh \left(\displaystyle \frac{x+{x}_{{\rm{c}}}}{2\delta }\right)\right)\right.\\ & & -\,\left.\arctan \left(\tanh \left(\displaystyle \frac{x-{x}_{{\rm{c}}}}{2\delta }\right)\right)\right]\\ & & +\,\displaystyle \frac{1}{2}\left(\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)\left(\tanh \left(\displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }\right)-\tanh \left(\displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }\right)\right)+C.\end{array}\end{eqnarray} \tag{ A2 }$$

To solve for the constant of integration C, we apply the condition from Equation (5):

$$\begin{eqnarray}&&{B}_{z}(\infty )=-{B}_{z}(0)={B}_{0}.\end{eqnarray} \tag{ A3 }$$

First, taking the limit as x → ∞ ,

$$\begin{eqnarray}\begin{array}{rcl}\mathop{\mathrm{lim}}\limits_{x\to \infty }{B}_{z}(x) & = & -4{\mu }_{0}{{en}}_{{\rm{b}}}{\beta }_{0}c\delta \left(\arctan \left(1\right)-\arctan \left(1\right)\right.\\ & & +\left.\,\displaystyle \frac{1}{2}\left(\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)\left(1-1\right)\right)+C\\ & = & C,\end{array}\end{eqnarray} \tag{ A4 }$$

because ${\mathrm{lim}}_{x\to \infty }\tanh (x)=1$. To double precision, $\tanh (x)=1$ if x ≥ 20, so the following is true if x_c/δ ≥ 40, which is the case in our simulations:

$$\begin{eqnarray}\begin{array}{rcl}{B}_{z}(0) & = & -4{\mu }_{0}{{en}}_{{\rm{b}}}{\beta }_{0}c\delta \left(\arctan \left(1\right)-\arctan \left(-1\right)\right.\\ & & +\left.\,\displaystyle \frac{1}{2}\left(\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)\left(1-(-1)\right)\right)+C\end{array}\end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray}&&=\ -4{\mu }_{0}{{en}}_{{\rm{b}}}{\beta }_{0}c\delta \left(\displaystyle \frac{\pi }{2}+\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)+C.\end{eqnarray} \tag{ A6 }$$

Combining the conditions on B_z as x → ∞ and at x = 0 (Equation (5)):

$$\begin{eqnarray}&&{B}_{0}=\mathop{\mathrm{lim}}\limits_{x\to \infty }{B}_{z}(x)=C=-{B}_{z}(0)\end{eqnarray} \tag{ A7 }$$

$$\begin{eqnarray}&&C=4{\mu }_{0}{{en}}_{{\rm{b}}}{\beta }_{0}c\delta \left(\displaystyle \frac{\pi }{2}+\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)-C\end{eqnarray} \tag{ A8 }$$

$$\begin{eqnarray}&&C=2{\mu }_{0}{{en}}_{{\rm{b}}}{\beta }_{0}c\delta \left(\displaystyle \frac{\pi }{2}+\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right).\end{eqnarray} \tag{ A9 }$$

This yields the expression for magnetic field in Equations (6) and (7):

$$\begin{eqnarray}\begin{array}{rcl}{B}_{z}(x) & = & -\displaystyle \frac{2{B}_{0}}{\left(\tfrac{\pi }{2}+\tfrac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)}\\ & & \times \,\left(\arctan \left(\tanh \left(\displaystyle \frac{x+{x}_{{\rm{c}}}}{2\delta }\right)\right)-\arctan \left(\tanh \left(\displaystyle \frac{x-{x}_{{\rm{c}}}}{2\delta }\right)\right)-\displaystyle \frac{\pi }{4}\right.\\ & & +\,\left.\displaystyle \frac{1}{2}\left(\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)\left(\tanh \left(\displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }\right)-\tanh \left(\displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }\right)-1\right)\right),\end{array}\end{eqnarray} \tag{ A10 }$$

with

$$\begin{eqnarray}&&{B}_{0}=2{\mu }_{0}{{en}}_{{\rm{b}}}{\beta }_{0}c\delta \left(\displaystyle \frac{\pi }{2}+\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right).\end{eqnarray} \tag{ A11 }$$

The full expressions for the equilibrium fields and plasma properties are

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{B(x)}{{B}_{0}} & = & \displaystyle \frac{2}{\tfrac{\pi }{2}+\tfrac{{n}_{d}}{{n}_{{\rm{b}}}}-1}\\ & & \times \,\left[\arctan \left(\tanh \left(\displaystyle \frac{x-{x}_{{\rm{c}}}}{2\delta }\right)\right)-\arctan \left(\tanh \left(\displaystyle \frac{x+{x}_{{\rm{c}}}}{2\delta }\right)\right)+\displaystyle \frac{\pi }{4}\right.\\ & & +\,\left.\displaystyle \frac{1}{2}\left(\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)\left(\tanh \left(\displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }\right)-\tanh \left(\displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }\right)+1\right)\right]\end{array}\end{eqnarray} \tag{ A12 }$$

$$\begin{eqnarray}&&\displaystyle \frac{n(x)}{{n}_{{\rm{b}}}}=1+\left(\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)\left({\rm{sech}} \left(\displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }\right)+{\rm{sech}} \left(\displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }\right)\right)\end{eqnarray} \tag{ A13 }$$

$$\begin{eqnarray}&&\beta (x)=\displaystyle \frac{{\lambda }_{{\rm{e}}}\sqrt{\sigma }}{2\delta \left(\tfrac{\pi }{2}+\tfrac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)}\sqrt{\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}}\left({\rm{sech}} \left(\displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }\right)-{\rm{sech}} \left(\displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }\right)\right)\end{eqnarray} \tag{ A14 }$$

$$\begin{eqnarray}&&\theta (x)=\displaystyle \frac{\sigma }{4}\displaystyle \frac{(4+\eta )/\eta -{[{B}_{z}(x)/{B}_{0}]}^{2}}{[n(x)/{n}_{{\rm{b}}}]\sqrt{1-\beta {\left(x\right)}^{2}}}.\end{eqnarray} \tag{ A15 }$$

A.2. Magnetic Field Perturbation

The magnetic field perturbation is based on the the vector potential (Equation (20)), which is an integral of B_z (x) (Equation (6)). We split B(x) into three types of terms, based on functional form: constants, $\arctan (\tanh )$ terms, and tanh terms. Constants are straightforward to integrate. The integrals of the other two types of terms are

$$\begin{eqnarray}&&\int \tanh \left(\displaystyle \frac{x\pm {x}_{{\rm{c}}}}{\delta }\right){dx}\equiv {\zeta }_{\pm }(x)=\delta \mathrm{log}\,\cosh \left(\displaystyle \frac{x\pm {x}_{{\rm{c}}}}{\delta }\right)\end{eqnarray} \tag{ A16 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \int \mathrm{arctan\; tanh}\left(\displaystyle \frac{x\pm {x}_{{\rm{c}}}}{2\delta }\right)\,{dx}\equiv {\xi }_{\pm }(x)\\ =\,2\delta \left(\displaystyle \frac{x\pm {x}_{{\rm{c}}}}{2\delta }\right)\left[\arctan \left(\exp \left(-\displaystyle \frac{x\pm {x}_{{\rm{c}}}}{\delta }\right)\right)+\arctan \left(\tanh \left(\displaystyle \frac{x\pm {x}_{{\rm{c}}}}{2\delta }\right)\right)\right]\\ \ +\,\displaystyle \frac{\delta i}{2}\left[{\mathrm{Li}}_{2}\left(-i\exp \left(-\displaystyle \frac{x\pm {x}_{{\rm{c}}}}{\delta }\right)\right)-{\mathrm{Li}}_{2}\left(i\exp \left(-\displaystyle \frac{x\pm {x}_{{\rm{c}}}}{\delta }\right)\right)\right].\end{array}\end{eqnarray} \tag{ A17 }$$

The value of ∫B_z (x) dx is

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \int {B}_{z}(x)\,{dx}=\displaystyle \frac{2{B}_{0}}{\tfrac{\pi }{2}+\tfrac{{n}_{d}}{{n}_{{\rm{b}}}}-1}\\ \times \,\left[{\xi }_{-}(x)-{\xi }_{+}(x)+\displaystyle \frac{\pi }{4}x+\displaystyle \frac{1}{2}\left(\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)\left({\zeta }_{-}(x)-{\zeta }_{+}(x)+x\right)\right]\end{array}\end{eqnarray} \tag{ A18 }$$

$$\begin{eqnarray}&&\begin{array}{l}=\,\displaystyle \frac{2{B}_{0}}{\tfrac{\pi }{2}+\tfrac{{n}_{d}}{{n}_{{\rm{b}}}}-1}\left[\displaystyle \frac{\pi }{4}x+(x-{x}_{{\rm{c}}})\left(\mathrm{arctan\; exp}\left(-\displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }\right)+\mathrm{arctan\; tanh}\left(\displaystyle \frac{x-{x}_{{\rm{c}}}}{2\delta }\right)\right)\right.\\ -\,(x+{x}_{{\rm{c}}})\left(\mathrm{arctan\; exp}\left(-\displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }\right)+\mathrm{arctan\; tanh}\left(\displaystyle \frac{x+{x}_{{\rm{c}}}}{2\delta }\right)\right)\\ +\,\displaystyle \frac{\delta i}{2}\left({\mathrm{Li}}_{2}\left(-i\exp \left(-\displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }\right)\right)-{\mathrm{Li}}_{2}\left(i\exp \left(-\displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }\right)\right)\right)\\ -\,\displaystyle \frac{\delta i}{2}\left({\mathrm{Li}}_{2}\left(-i\exp \left(-\displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }\right)\right)-{\mathrm{Li}}_{2}\left(i\exp \left(-\displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }\right)\right)\right)\\ +\,\left.\displaystyle \frac{1}{2}\left(\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)\left(\delta \mathrm{log\; cosh}\left(\displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }\right)-\delta \mathrm{log\; cosh}\left(\displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }\right)+x\right)\right]\end{array}\end{eqnarray} \tag{ A19 }$$

$$\begin{eqnarray}&&\begin{array}{l}=\,\displaystyle \frac{2{B}_{0}}{\tfrac{\pi }{2}+\tfrac{{n}_{d}}{{n}_{{\rm{b}}}}-1}\left[\displaystyle \frac{\pi }{4}(x-2{x}_{{\rm{c}}})+\displaystyle \frac{1}{2}\left(\displaystyle \frac{{n}_{d}}{{n}_{{\rm{b}}}}-1\right)\left(\delta \mathrm{log\; cosh}\left(\displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }\right)-\delta \mathrm{log\; cosh}\left(\displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }\right)+x\right)\right.\\ +\,\displaystyle \frac{\delta i}{2}\left({\mathrm{Li}}_{2}\left(-i\exp \left(-\displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }\right)\right)-{\mathrm{Li}}_{2}\left(i\exp \left(-\displaystyle \frac{x-{x}_{{\rm{c}}}}{\delta }\right)\right)\right)\\ -\,\left.\displaystyle \frac{\delta i}{2}\left({\mathrm{Li}}_{2}\left(-i\exp \left(-\displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }\right)\right)-{\mathrm{Li}}_{2}\left(i\exp \left(-\displaystyle \frac{x+{x}_{{\rm{c}}}}{\delta }\right)\right)\right)\right].\end{array}\end{eqnarray} \tag{ A20 }$$

Equation (A20) results from applying the simplification

$$\begin{eqnarray}&&\arctan \left({e}^{-2y}\right)+\arctan \left(\tanh (y)\right)=\displaystyle \frac{\pi }{4}.\end{eqnarray} \tag{ A21 }$$

For the purposes of symmetry and preserving the periodic boundary conditions, we have set the constant of integration equal to zero.

The perturbation to the magnetic field is

$$\begin{eqnarray}&&{\boldsymbol{B}}^{\prime} ={\boldsymbol{\nabla }}\times {A}_{y}-{\boldsymbol{B}}=-\displaystyle \frac{\partial {A}_{y}}{\partial z}\hat{{\boldsymbol{i}}}+\displaystyle \frac{\partial {A}_{y}}{\partial x}\hat{{\boldsymbol{k}}}-{\boldsymbol{B}}\end{eqnarray} \tag{ A22 }$$

$$\begin{eqnarray}\begin{array}{rcl} & = & -0.01\times \hat{{\boldsymbol{i}}}\left[\displaystyle \frac{51\pi }{{L}_{z}}\sin \left(\displaystyle \frac{\pi }{{L}_{z}}z\right)\right.\\ & & \left.\times \,{\cos }^{50}\left(\displaystyle \frac{\pi }{{L}_{z}}z\right){\cos }^{2}\left(\displaystyle \frac{\pi (x-{x}_{{\rm{c}}})}{{L}_{x}}\right)\displaystyle \int {B}_{z}(x)\,{dx}\right]\\ & & +\,0.01\times \hat{{\boldsymbol{k}}}\left[\displaystyle \frac{\pi }{{L}_{x}}{\cos }^{51}\left(\displaystyle \frac{\pi }{{L}_{z}}z\right)\sin \left(\displaystyle \frac{2\pi (x-{x}_{{\rm{c}}})}{{L}_{x}}\right)\right.\\ & & \times \,\displaystyle \int {B}_{z}(x){dx}-{\cos }^{51}\left(\displaystyle \frac{\pi }{{L}_{z}}z\right){\cos }^{2}\left(\displaystyle \frac{\pi (x-{x}_{{\rm{c}}})}{{L}_{x}}\right)\\ & & \left.\times {B}_{z}(x)\right].\end{array}\end{eqnarray} \tag{ A23 }$$