2.1 Field structure

We use the 3-D magnetic field and electric potential solutions to the steady-state, kinematic, resistive MHD equations given by Wyper & Jain (Reference Wyper and Jain2010) and Wyper & Jain (Reference Wyper and Jain2011). In the torsional spine case, the twist in the magnetic field is localized to the fan plane and centred about the null point (figure 1b). All field structures are given in cylindrical coordinates $[r, \theta, z]$. The torsional spine magnetic field is of the form

(2.1)\begin{equation} \boldsymbol{B} = B_{0} [ r, j r^{\alpha} ( z r^{2} )^{\beta} e^{ -( 1 / l^{2} ) [ r^{2} + c^{2} ( z r^{2} )^{2} ] }, -2 z ], \end{equation}

with distances scaled to a characteristic length $L_{0}$ (Chesny & Orange Reference Chesny and Orange2020; Chesny, Orange & Dempsey Reference Chesny, Orange and Dempsey2021a; Chesny et al. Reference Chesny, Orange and Hatfield2021b). The free parameters, $j$, $l$ and $c$ govern the degree and spatial extent of the twist, and $\alpha$ and $\beta$ are positive integers. Throughout this paper, we set $j = 3$, $\alpha = 4$, $l = 1$ and $\beta = c = 0$, corresponding to a generic radial perturbation, and consistent with the TSR simulations of Hosseinpour et al. (Reference Hosseinpour, Mehdizade and Mohammadi2014) and Hosseinpour (Reference Hosseinpour2014b). Using these values, (2.1) reduces to

(2.2)\begin{equation} \boldsymbol{B} = B_{0} [ r, 3 r^{4} e^{{-}r^{2}}, -2z ]. \end{equation}

Figure 2(a,b) displays $XZ$ and $XY$ streamlines of (2.2), with the distance scaled to a characteristic length (see below). The location of greatest twist is at the radial distance $R = \sqrt {2}$ defined by setting ${\rm d}B_{\theta }/{\rm d}r = 0$ in (2.2). To approximate a controlled injection of plasma, we apply a localized resistivity profile $\eta$ depending only on radius $r$ as

(2.3)\begin{equation} \eta = \eta_{0} r^{\lambda} \mathrm{e}^{{-}r^{2}/l^{2}}, \end{equation}

where $\eta _{0} = m_{e} \nu _{e,i} / n_{e} e^{2}$ is the plasma resistivity calculated from the characteristic parameters of electron mass $m_{e}$, electron–ion collision frequency $\nu _{e,i}$, electron number density $n_{e}$ and the elementary charge $e$ (see § 2.3). We set $\lambda = 2$ to localize the resistivity profile to the radius $R = \sqrt {2}$ of maximum magnetic field twist predicted by (2.2). This choice, as will be shown in § 2.3, inherently confines and localizes the injected plasma sheath to the region of stressed, non-potential magnetic field. This choice also satisfies the constraints for realistic solutions when the sum of parameters $\alpha + \lambda \ge 4$ and is even (Wyper & Jain Reference Wyper and Jain2011). The electric potential in TSR with the above resistivity profile is of the form

(2.4)\begin{equation} \varPhi = \varPhi_{0} \left[ - \frac{4}{l^{2}} F(\alpha + \lambda - 1) + 2 (\alpha + 1) F(\alpha + \lambda - 3) \right] z r^{2}, \end{equation}

with the recurrence relation for $F(r, A)$ being

(2.5)\begin{equation} F(A + 2) = \frac{l^{2}}{2} \left( A F(A) - \frac{r^{A}}{B_{0}} e^{{-}r^{2}/l^{2}} \right), \end{equation}

and solutions for all values of $A$ found using the error function (erf) as

(2.6)\begin{equation} F(r, 1) = \frac{l \sqrt{\rm \pi}}{2 B_{0}} \mathrm{erf}\left( \frac{r}{l} \right), \end{equation}

and

(2.7)\begin{equation} F(r, 2) ={-}\frac{l^{2}}{2 B_{0}} e^{{-}r^{2}/l^{2}}. \end{equation}

Figure 2.

Vector field quantities of resistive TSR (Wyper & Jain Reference Wyper and Jain2010). The blue markers denote the radius of maximum magnetic perturbation curvature at $[r, \theta, \phi ] = [\sqrt {2}, 0, 0]$. Panels show the (a) $XZ$ magnetic field (Y=0), (b) $XY$ magnetic field (Z=0), (c) $XZ$ electric field (Y=0), (d) $XY$ electric field (Z=0). (e) $XZ$ drift velocity (Y=0) computed from $\boldsymbol {E} \times \boldsymbol {B}$, (f) $XY$ drift velocity (Z=0).

As in Wyper & Jain (Reference Wyper and Jain2011), the initial electric potential $\varPhi _{0}$ is in units of volts and depends on the current density $j_{0}$ and magnetic field as $\varPhi _{0} = j_{0} \eta _{0} B_{0}$ (Wyper & Jain Reference Wyper and Jain2010; Chesny et al. Reference Chesny, Orange and Hatfield2021b). The electric field can then be derived using $\boldsymbol {E} = - \boldsymbol {\nabla } \varPhi$. Using our numerical parameters, these equations reduce to

(2.8)\begin{equation} \varPhi = \varPhi_{0} \left[ \sqrt{\rm \pi} \, \mathrm{erf} \left( \frac{r}{l} \right) + 2 z (r^{5} - r^{3}) e^{{-}r^{2}/l^{2}} \right], \end{equation}

and in cylindrical units $[r, \theta, z]$

(2.9)\begin{equation} \boldsymbol{E} = E_{0} \left[ 4 z r^{6} - 14 z r^{4} + 6 z r^{2} - 2, 0, 2 ({-}r^{5} + r^{3}) \right] e^{- r^{2} / l^{2}}. \end{equation}

Figure 2 displays cross-sections of the vector fields of the torsional spine magnetic field ((a,b), (2.2)), electric field ((c,d), derived from (2.8)) and the drift velocity ((e,f), (2.10)) calculated from

(2.10)\begin{equation} v_{d} = \frac{\boldsymbol{E} \times \boldsymbol{B}}{B^{2}}. \end{equation}

It is important to understand that, in torsional reconnection modes, we do not consider the particle dynamics close to the null point where $|B|=0$ and $v_{d} \approx 0$ within which the gyroradius would be large. Vector maps of the drift velocity (figure 2e,f) demonstrate that $v_{d} \rightarrow 0$ close to the null, whereas significant magnitudes of $v_{d}$ occur away from the null, which is where we consider particle dynamics and acceleration. Blue contours denote the radius of maximum perturbing magnetic field twist applied to the fan plane at $r = \sqrt {2}$. Figure 3 shows the scalar fields of electric potential ((a,b), (2.8)), and plasma resistivity profile ((c,d), (2.3)) localized to the region of maximum magnetic twist. The localized resistivity profile demonstrates the spine-aligned, tubular current sheet of the TSR mode (Pontin & Galsgaard Reference Pontin and Galsgaard2007).

Figure 3.

Scalar field quantities of resistive TSR (Wyper & Jain Reference Wyper and Jain2010). The blue markers denote the radius of maximum magnetic perturbation curvature at $[r, \theta, \phi ] = [\sqrt {2}, 0, 0]$. Panels show (a) $XZ$ electric potential (Y=0), (b) $XY$ electric potential (Z=0), (c) $XZ$ resistivity (Y=0), (d) $XY$ resistivity (Z=0).

The time-dependent dynamics leading to the breakdown of ideal MHD conditions and the slipping of stressed magnetic field lines through the plasma to generate the electric potential, electric field and resistivity profiles of figures 2 and 3 are details that have been explored by the resistive MHD simulations of Pontin et al. (Reference Pontin, Al-Hachami and Galsgaard2011) and Pontin & Galsgaard (Reference Pontin and Galsgaard2007), while the equations of § 2.1 are solutions assuming that torsional reconnection has already occurred (Wyper & Jain Reference Wyper and Jain2010, Reference Wyper and Jain2011). Therefore, it is of interest to the experimental plasma physics community as to how to set up and initiate these conditions in a laboratory setting. Section 2.3 introduces one possible method by which the initial global magnetic field and a plasma with an embedded magnetic field can superimpose to form the general TSR magnetic field of figures 1(b) and 2(a,b), consistent with the laboratory TFR parameters of Chesny et al. (Reference Chesny, Orange and Hatfield2021b). While this manuscript is focused on addressing the efficacy of TSR as a particle accelerator under typical laboratory settings, the time-dependent plasma and magnetic field dynamics of scaled laboratory hardware component models are left for more advanced simulation methods, including particle-in-cell (PIC).

2.2 Particle trajectories

Particle trajectories are computed using the relativistic equations of motion for velocity and the Lorentz force in a collisional plasma, respectively,

(2.11)\begin{gather} \frac{{\rm d} \boldsymbol{r}}{{\rm d}t} = \frac{\boldsymbol{p}}{\gamma A m_{p}}, \end{gather}

(2.12)\begin{gather} \frac{{\rm d} \boldsymbol{p}}{{\rm d}t} = e \left( \boldsymbol{E} + \frac{\boldsymbol{p}}{\gamma A m_{p}} \times \boldsymbol{B} \right) - \nu_{in} \boldsymbol{p}. \end{gather}

In these equations, $\boldsymbol {p}$ is momentum, $\boldsymbol {r}$ is the 3-D Cartesian position, $\gamma = (1 - \boldsymbol {p} \boldsymbol {\cdot } \boldsymbol {p} )^{-1/2}$ is the relativistic Lorentz factor, $e$ is the charge of a singly ionized particle, $A$ is the atomic mass, $m_{p}$ is the proton rest mass in MeV c$^{-2}$ and $\nu _{in}$ is the ion–neutral collision frequency (see below). Equations (2.11) and (2.12) can be solved in dimensionless form by defining the following characteristic relations:

(2.13a)\begin{equation} \boldsymbol{r} = L_{0} \boldsymbol{x} ; \quad \boldsymbol{p} = p_{0} \boldsymbol{n}, ; \quad \boldsymbol{B} = B_{0} \boldsymbol{b}, ; \quad t = \tau_{0} t^{\prime}, \end{equation}

where $\boldsymbol {x}$ is a 3-D position unit vector, $\boldsymbol {n}$ is a unit vector in the direction of motion, $\boldsymbol {b}$ is a unit vector defining the local magnetic field direction and the time step $t^{\prime }$ is scaled to the characteristic non-relativistic gyroperiod $\tau _{0} = 2 {\rm \pi}m_{0} / e B_{0}$ (Dalla & Browning Reference Dalla and Browning2005, Reference Dalla and Browning2006, Reference Dalla and Browning2008). Using these relations, the dimensionless equation of motions for velocity and the Lorentz force can be written as respectively,

(2.14)\begin{gather} \frac{{\rm d} \boldsymbol{x}}{{\rm d}t^{\prime}} = \frac{c_{1}}{\gamma} \boldsymbol{n}, \end{gather}

(2.15)\begin{gather} \frac{{\rm d} \boldsymbol{n}}{{\rm d}t^{\prime}} = 2 {\rm \pi}\left( c_{2}\boldsymbol{E} + \frac{\boldsymbol{n} \times \boldsymbol{b}}{\gamma} \right) - c_{3} \boldsymbol{n}, \end{gather}

with the unitless velocity, electric field and collisional constants represented as, respectively,

(2.16a)\begin{equation} c_{1} = \frac{p_{0} \tau_{0} |c|}{A m_{p} L_{0}} ; \quad c_{2} = \frac{A m_{p}}{p_{0} B_{0} |c|} ;\quad c_{3} = \nu_{in} \tau_{0}, \end{equation}

and $|c| = 3\times 10^{8}$ is the magnitude of the speed of light in m s$^{-1}$. The initial momentum $p_{0}$ in MeV c$^{-1}$ is calculated from the initial particle kinetic energy $K_{0}$ in eV as $p_{0} = (K_{0}^{2} + 2 A m_{p} K_{0} )^{1/2}$. Our choices of $L_{0}$ and $B_{0}$ will be described in § 2.3.

Since full orbit particle tracing in plasmas is not limited by inelastic collisions (e.g. Marchand Reference Marchand2010; Homma et al. Reference Homma, Hoshino, Tokunaga, Yamoto, Hatayama, Asakura, Sakamoto and Tobita2018), the effects of collisions are accounted for in the final term on the right-hand side of (2.12). The addition of this collisional term has been shown to be effective in particle tracking studies of magnetic reconnection in high density laser plasmas up to ${\sim }10^{25}$ m$^{-3}$ (Zhong et al. Reference Zhong, Lin, Li, Wang, Li, Zhang, Yuan, Ping, Wei and Wang2016). The magnitude of the collision frequency is solved from the plasma parameters in § 2.3. Most laboratory plasmas are weakly ionized, where the degree of ionization is dependent on the ion ($n_{i}$) and neutral gas ($n_{g}$) number densities $\chi = n_{i} / ( n_{g} + n_{i} )$ and $\chi \ll 1$. Weakly ionized plasmas are further defined by the electron–ion collision frequency $\nu _{ei}$ being less than the electron–neutral collision frequency $\nu _{en}$ ($\nu _{ei} < \nu _{en}$; Park et al. Reference Park, Choe, Moon and Yoo2019). In the case of low-energy electrons (a few eV), the main collisional process in a weakly ionized plasma is dominated by polarization scattering of ions against ambient neutral particles $\nu _{in}$ (Lieberman & Lichtenberg Reference Lieberman and Lichtenberg2005). Electron and ion polarization scattering rate constants in units of cm$^{3}$ s$^{-1}$ are calculated from

(2.17)\begin{equation} K_{Le} = 3.85 \times 10^{{-}8} \alpha_{R}^{1/2}, \end{equation}

and

(2.18)\begin{equation} K_{Li} = 8.99 \times 10^{{-}10} \left( \frac{\alpha_{R}}{A} \right)^{1/2}, \end{equation}

respectively, where $\alpha _{R}$ is the species-dependent relative polarizability ($\alpha _{R, He} = 1.384$ for helium; Schwerdtfeger & Nagle Reference Schwerdtfeger and Nagle2019). The charged particle collision frequencies against neutrals are then calculated using

(2.19)\begin{equation} \nu_{e,i-n} = n_{g} K_{L-e,i}. \end{equation}

Parameters using typical laboratory values given in table 1 demonstrate numerically that the chosen helium plasma is weakly ionized ($\chi = 1$% and $\nu _{ei} < \nu _{en}$) and polarization scattering is the dominant collisional process. We note that (2.19) is independent of the background ion temperature of the plasma and that temperature-dependent collision frequencies are left for more advanced plasma simulation beyond particle tracking, such as PIC. In fact, considering the temperature-independent collision frequency as in this work may represent the worst-case scenario for particle acceleration, since particles with higher kinetic energies have a larger mean free path between collisions. In the TSR case here, and in the TFR case of Chesny et al. (Reference Chesny, Orange and Hatfield2021b), the final particle temperatures could actually become substantially higher by considering a temperature-dependent collision frequency.

Table 1.

Laboratory parameters used for the CPG helium plasma regime. The model assumes a capacitor bank discharge of 6 kV over 10 cm ($E=60$ kV m$^{-1}$) (Chesny et al. Reference Chesny, Orange and Hatfield2021b).

The dimensionless equations (2.14) and (2.15) are solved numerically using the scipy.integrate.odeint package which uses LSODA from the FORTRAN library to dynamically monitor data and reduce errors by automatically switching between stiff and non-stiff methods. Numerical validation of this routine first verified the relativistic ion gyro-radius $r_{g,i} = \gamma A m_{p} v_{\bot } / e B$ within uniform magnetic fields ($\boldsymbol {E} = 0$ V m$^{-1}$). Then, summing the kinetic and electric potential energies throughout test runs with collisionless single particle trajectories when $\boldsymbol {E} \neq 0$ V m$^{-1}$ demonstrated total energy conservation to five significant figures, of the order of $0.1$ eV, using (Dalla & Browning Reference Dalla and Browning2005, Reference Dalla and Browning2006, Reference Dalla and Browning2008; Gascoyne Reference Gascoyne2015)

(2.20)\begin{equation} W = K + e \varPhi. \end{equation}

Monitoring the total energy in the collisional plasma regime will elucidate the energy loss and energy loss rate due to collisions and scattering during TSR.

2.3 Laboratory plasma scaling parameters and particle injection

In our simulations, the laboratory plasma parameter space follows the typical profile of a coaxial plasma gun (CPG; Marshall Reference Marshall1960), which is a common plasma generation device across a range of fields from plasma jet magneto-inertial fusion (PJMIF, e.g. Yates et al. Reference Yates, Langendorf, Hsu, Dunn, Brockington, Case, Cruz, Witherspoon, Thio and Cassibry2020) to in-space electric propulsion. Previous numerical investigations of torsional reconnection (Chesny et al. Reference Chesny, Orange, Oluseyi and Valletta2017, Reference Chesny, Orange and Hatfield2021b) also used CPG plasma sheath parameters to investigate feasible experimental particle acceleration by exploiting their general azimuthal magnetic field profile as a source of magnetic field twist (2.2). The CPGs are cylindrically symmetric, pulse power devices consisting of a central anode rod surrounded by an annular outer electrode and separated by an insulator (e.g. Larson, Liebing & Dethlefsen Reference Larson, Liebing and Dethlefsen1966; Witherspoon et al. Reference Witherspoon, Case, Messer, Bomgardner, Phillips, Brockington and Elton2009; Subramaniam & Raja Reference Subramaniam and Raja2017). The CPG plasma profiles follow the snowplough model in which a capacitor bank discharge forms a thin plasma sheath along the insulator that then propagates down the electrodes due to the self $\boldsymbol {J} \times \boldsymbol {B}$ force (e.g. Hart Reference Hart1964; Lee & Serban Reference Lee and Serban1996; Gonzalez et al. Reference Gonzalez, Clausse, Bruzzone and Florido2004). This CPG is theorized to be located within a global fan-spine magnetic field, as enabled by conventional pulse power conducting coil technology (Chesny et al. Reference Chesny, Orange and Dempsey2021a, Reference Chesny, Moffett, Cole, Baptiste and Orange2022a), centred at the null point, and axially aligned with the spine axis.

In this study, individual particle injection positions vary radially due to the localized perturbation in the fan plane taken to be the embedded azimuthal magnetic field of a CPG-based plasma sheath (Moffett et al. Reference Moffett, Chesny, Cole, Hatfield and Rusovici2022). Due to the drift velocity calculations shown in figure 2(e), the only asymmetry in the axial vector components occur within the cylindrical $r = \pm L_{0}$ boundary and point in the $-Z$ direction only. Thus, we inject each particle with an identical initial momenta in the $-Z$ direction (kinetic energy $K_{0} = 10$ eV; e.g. Schaer Reference Schaer1994; Chesny et al. Reference Chesny, Orange and Hatfield2021b). This approximates to the bulk plasma velocity of a CPG-based sheath due to the $\boldsymbol {J} \times \boldsymbol {B}$ force in the snowplough mode (e.g. Kwek, Tou & Lee Reference Kwek, Tou and Lee1990). The authors note that other pitch angles, including along $+Z$, were attempted and which demonstrated that peak morphological behaviour and peak kinetic energy gain are invariant. Thus, the choice of displaying solutions to $-Z$ injection presents no loss of generality. Due to the cylindrical symmetry of (2.2) and (2.9), we solved for individual particle trajectories at five representative input radii along the $x$ axis $r_{x} = [ 0.50 \sqrt {2}, 0.75 \sqrt {2}, \sqrt {2}, 1.25 \sqrt {2}, 1.50 \sqrt {2} ]$. In § 3, we divide the results into three representative ‘regions:’ region 1 ($r < \sqrt {2} L_{0}$), region 2 ($r = \sqrt {2} L_{0}$) and region 3 ($r > \sqrt {2} L_{0}$). This choice of three representative regions showing five discrete trajectories, similar to the analyses of Hosseinpour (Reference Hosseinpour2014a,Reference Hosseinpourb, Reference Hosseinpour2015) and Hosseinpour et al. (Reference Hosseinpour, Mehdizade and Mohammadi2014) at solar scales rather than a bulk particle input (Dalla & Browning Reference Dalla and Browning2005, Reference Dalla and Browning2006, Reference Dalla and Browning2008; Chesny et al. Reference Chesny, Orange and Hatfield2021b), is necessary due to the dynamics of collisional energy losses. It was found that the particles injected into different regions required different total simulation times to properly diagnose the energy balance. Regardless, we will discuss in § 4 how bulk particle behaviour can be inferred from these representative particles.

Table 1 shows the experimentally derived numerical input plasma parameters for our relativistic particle tracking algorithm. Following the TFR analysis of Chesny et al. (Reference Chesny, Orange and Hatfield2021b), the parameters of a 6 kV capacitive discharge over 10 cm using helium fuel (Schaer Reference Schaer1994) and a magnetic field strength of 1 T follow achievable experimental values (Chesny et al. Reference Chesny, Moffett, Cole, Baptiste and Orange2022a,Reference Chesny, Moffett, Hatfield, Cole, Landers, Shokrollahi and Egeb). These are the parameters used in the scalar and vector field plots in figures 2 and 3. By inspecting the drift velocity vector plots of figure 2(e,f), particle trajectory behaviour can be anticipated when the injection positions are considered. For region 1 ($r < \sqrt {2}$), particles should be accelerated in the $-Z$ direction. The $XY$ plane shows that the particles may be directed toward and rotationally around the central spine axis, and the $XZ$ plane shows that the drift may rapidly end when the particles approach the spine axis before $Z \approx -1$. In region 2 ($r = \sqrt {2}$), the drift velocity pushes out along the $XY$ fan plane and down toward the fan in the $XZ$ plane. In region 3 ($r > \sqrt {2}$), the $XY$ drift velocity seems to push particles out along the fan while still spiralling around the central spine axis. The $XZ$ drift suggests an additional mild force along the $+Z$ direction. The highest acceleration due to the peak drift velocity profile should occur with particles that enter the vicinity $x, z \approx [\pm 0.50, -0.50]$.