4.1. Numerical set-up

The OpenMHD code is described by Zenitani & Miyoshi Reference Zenitani and Miyoshi2011 and has been used to study instabilities, shocks, turbulence and magnetic reconnection in low $\beta$ plasmas, where $\beta$ is the ratio of the thermal to magnetic pressure (Zenitani Reference Zenitani2015; Shimizu, Kondoh & Zenitani Reference Shimizu, Kondoh and Zenitani2017; Hosseinpour, Chen & Zenitani Reference Hosseinpour, Chen and Zenitani2018; Zenitani & Miyoshi Reference Zenitani and Miyoshi2020). OpenMHD is open-source and available for download at https://sci.nao.ac.jp/MEMBER/zenitani/openmhd-e.html. The resistive MHD equations used in this framework are made dimensionless with $\mu _{0}=1$ , enforcing the condition (Pontin & Galsgaard Reference Pontin and Galsgaard2007) that with dimensionless quantities of mass density and magnetic field strength $|\rho |=|\boldsymbol {B}|=1$ , respectively, an Alfvén wave will travel the full domain of length $L$ in an Alfvén time $\tau _{A} = L \sqrt {\mu _{0} \rho _{0}} / B_{0}$ as calculated using the dimensional quantities of $\mu _{0}$ , $\rho _{0}$ and $B_{0}$ . OpenMHD is written in the Cartesian coordinate system ( $x, y, z$ ), but to achieve cylindrical coordinates ( $r, \theta , z$ ) we claim the equivalence of the axisymmetric ( $r, z$ ) plane to ( $x, z$ ) at $y=0$ .

To the authors’ knowledge, this study represents the first implementation of physical BCs in using OpenMHD. Upon request, the authors will make the source code available to other researchers. We base the CPG model code on the example case $\mathsf {2D\_reconnection}$ available in the source code due to its inclusion of plasma resistivity and turbulence. The physical dimensions of the experimental device were duplicated in the MHD domain, as shown in figure 1 and described in § 2. The authors used the results of the global shutter imaging in figure 5, top left frame at 0.0 $\unicode{x03BC}$ s, to approximate the initial conditions of the plasma as detailed quantitatively in the system of equations (4.1) to (4.6). Since there are no diagnostics deployed on the CPG to quantify any plasma parameters, equations (4.1) to (4.6) are considered dimensionless. Successfully approximating the physical morphology of the plasma propagation with respect to the form of the initial conditions is considered a primary result of this study. Previous results using OpenMHD have been considered dimensionless with respect to the plasma parameters, with the evolution of the fields and density profiles being the significant takeaways therein (Zenitani & Miyoshi Reference Zenitani and Miyoshi2011; Zenitani Reference Zenitani2015; Shimizu et al. Reference Shimizu, Kondoh and Zenitani2017; Zenitani & Miyoshi Reference Zenitani and Miyoshi2020). The goal of this study is to do the same during an application of the OpenMHD framework to an operational experimental device.

In OpenMHD, modifications are made to the following source files. The lead file $\mathsf {main.f90}$ sets integer values for the size of the ( $x, y$ ) mesh (axisymmetric $r, z$ ), the time step and duration of the simulation and the output interval of time-stamped data sets. In this file we define the integer values of the CPG electrodes geometry (figure 1) and the initial plasma volume.

Table 1.

The velocity $v$ , magnetic field $B$ , density $\rho$ and pressure $P$ BCs and definitions implemented for physical electrodes and walls in OpenMHD. The horizontal $x$ -direction (axisymmetric $z$ ) uses index $i$ , and the vertical $y$ -direction (axisymmetric $r$ ) uses index $j$ . The iterative index $n$ is included to define the parametric geometry of the circular nose cone shape. Note that Cartesian coordinates ( $x, y$ ) in OpenMHD correspond to ( $r, z$ ) in the axisymmetric CPG geometry.

Table 1. Long description

The table lists boundary conditions (BCs) for velocity, magnetic field, density, and pressure in OpenMHD. It includes conditions for top, bottom, left, right, nose cone right, nose cone top, and electrode volumes. Each BC is defined using specific equations involving indices and variables. The conditions are set to ensure proper implementation of physical electrodes and walls in the simulation.

BCs for all electrodes, walls and initial plasma states are set within the file $\mathsf {bc.f90}$ . The boundaries for the physical electrodes are ghost cells that reflect the velocity $v$ (momentum) conditions normal to the surface. This type of free-slip BC where the momentum vectors normal to the surface interaction are reversed is typical of electrode boundaries (Subramaniam & Raja Reference Subramaniam and Raja2017; Datta et al. Reference Datta, Meier and Shumlak2024). The top, bottom and right sides of the domain are momentum outflow BCs. The top of the domain is a vacuum wall but we imposed the momentum outflow BC to eliminate reflective motions that shroud the primary plasma morphology. Table 1 lists the BCs imposed for the CPG electrode definitions shown in figure 7. We note that the imposed electrode BCs are made to the conserved variables and not to the interpolated spatial profile states (Zenitani & Miyoshi Reference Zenitani and Miyoshi2011). To exclude magnetic flux from the electrodes, the magnetic field is set to zero within the electrode boundaries, i.e. $B_{x}=B_{y}=B_{z}=0$ throughout the simulation within the electrode boundaries. In a practical sense, the exclusion of flux represents a physical conductor, and the flee-slip condition mimics a physical barrier.

Figure 7.

Electrode BCs of the CPG geometry shown in figure 1. Ghost cell BCs are implemented on top, bottom, right and nose cone edges, with constant values within. The wall boundaries each allow plasma to pass through. Plasma BCs, including initial velocity, are defined at the $t=0$ time step only.

The initial conditions (ICs) of the simulation are set within the file $\mathsf {model.f90}$ and define the initial density, pressure, velocity and magnetic field throughout the simulation volume. The ICs are chosen to empirically match the observed morphological plasma profile in the first ionisation frame (0.0 $\unicode{x03BC}$ s in figure 5) and to be mathematically consistent with typical ICs in MHD modelling of CPGs. Viewing the initial ionisation plasma sheath from global shutter imaging (figure 5, 0 $\unicode{x03BC}$ s, top left), the mass density profile is observed to have a relatively wide profile and is numerically represented by the $\psi (r, z)$ function in (4.1). This stand-alone $\psi (r, z)$ function has a constant value along $r$ at any discrete axial $z$ value. Along any axial cross-section, it forms a normal distribution centred about small $z$ values ( $Z_{s}$ variable below). Multiplying $\psi (r, z)$ by an $r^{-\gamma }$ term in the subsequent profiles then gives each IC an additional radial dependence on top of $\psi (r, z)$ to qualitatively match the spatial profile of the initial ionisation frame. This $\psi (r, z)$ profile is then used to define the initial spatial distribution of the other plasma profiles in (4.2) to (4.6). The initial spatial distribution of the $\psi (r, z) r^{-\gamma }$ profile can be seen in the first row ( $t=0$ frames in landscape format) of figure 8 described in § 4.2

(4.1) \begin{align} \psi (r, z) &= \frac {{\tanh}(\sigma (r + w_{s}) - Z_{s}) - {\tanh}(\sigma (r - w_{s}) - Z_{s})}{2 \, {\tanh}(\sigma Z_{s})},\\[-10pt]\nonumber \end{align}

(4.2) \begin{align} \rho (r, z) &= \rho _{b} + \psi (r, z) r^{-\gamma },\\[-10pt]\nonumber \end{align}

(4.3) \begin{align} P(r, z) &= 0.5 \beta \rho (r),\\[-10pt]\nonumber \end{align}

(4.4) \begin{align} v_{z}(r, z) &= v_{z,0} \psi (r, z),\\[-10pt]\nonumber \end{align}

(4.5) \begin{align} B_{\theta }(r, z) &= B_{0} \psi (r, z) r^{-1},\\[-10pt]\nonumber \end{align}

(4.6) \begin{align} \eta (r, z) &= \eta _{0} + \eta _{1} \psi (r, z) r^{-\gamma }. \end{align}

In these ICs, $\rho _{b} = 0.01$ is the uniform mass density of the background, $\sigma = 3$ regulates the flatness of the profile, $w_{s} = 0.5$ is the width of the sheath, $Z_{s} = 2.421$ is the axial position of the sheath and the exponent $\gamma = 1$ determines the radial profile distribution which is key to the emergence of blow-by (Cassibry et al. Reference Cassibry, Thio and Wu2006). The plasma beta parameter $\beta$ is also set within this file, which is the ratio of thermal to magnetic pressure. Using the numerical value $\beta =0.2$ is identical to other OpenMHD studies (Zenitani & Miyoshi Reference Zenitani and Miyoshi2011; Zenitani Reference Zenitani2015). The numerical result of this initial OpenMHD sheath in comparison with the observed sheath is shown in the top row (left and middle) of figures 8 and 9. These images were split across two figures for clarity.

Figure 8.

First half of CPG imaging and OpenMHD results for times $t^{\prime }=0$ to $t^{\prime }=1.8$ presented in landscape format. (Left column) Red–blue numerical scaling of the raw global shutter imaging shown in figure 5 (left). These images are a 3-D morphology projected onto a 2-D plane. (Middle column) Axisymmetric OpenMHD profiles of density and azimuthal magnetic field $B_{\theta }$ (colour bars, respectively), which are identical in topologies and which act as cross-sectional analogues of imaging. (Right column) OpenMHD pressure profiles at even time increments of $\Delta t^{\prime } = 0.6$ .

Figure 9.

Second half of CPG imaging and OpenMHD results for times $t^{\prime }=2.4$ to $t^{\prime }=4.2$ presented in landscape format. (Left column) Red–blue numerical scaling of the raw global shutter imaging shown in figure 5 (left). These images are a 3-D morphology projected onto a 2-D plane. (Middle column) Axisymmetric OpenMHD profiles of density and azimuthal magnetic field $B_{\theta }$ (colorbars, respectively), which are identical in topologies and which act as cross-sectional analogues of imaging. (Right column) OpenMHD pressure profiles at even time increments of $\Delta t^{\prime } = 0.6$ .

Figure 9. Long description

The image contains multiple panels showing the evolution of density, pressure, and magnetic field in a plasma experiment over time. Panel A: The first column shows global shutter frames at different times (70.0 µs, 87.5 µs, 105.0 µs, and 122.5 µs) with the radial distance (R) on the x-axis and an unspecified unit on the y-axis. Panel B: The second column displays density profiles at times 2.4, 3.0, 3.6, and 4.2 with the spatial coordinate (Z) on the x-axis and radial distance (R) on the y-axis. Panel C: The third column shows pressure profiles at the same times as the density profiles, with the same axes. Panel D: The fourth column presents magnetic field profiles at the same times, with the same axes. Each panel uses a color scale to represent the intensity of the respective quantities, with the color bar indicating the range of values.

The choice of $v_{z, 0}=5$ is made so that the longest extent of the plasma simulation best matches the greatest extent of the plasma imaging. While the simulations are axisymmetric in the ( $r, z$ ) plane, as shown in figures 1 and 7, the results described in § 4.2 are plotted by reflecting the ( $+r, z$ ) simulation time steps into the ( $-r, z$ ) plane to match the full CPG geometry captured in the imaging (figure 5). The computational grid is set within both $\mathsf {main.f90}$ and $\mathsf {model.f90}$ , with a set size of ( $r, z$ ) = (4.60, 22.86) corresponding to the cgs units of the CPG domain, and with a mesh size of (604, 3000). The time step is set to $\Delta t = 0.10$ and is advanced using a second-order Runge–Kutta method (Zenitani & Miyoshi Reference Zenitani and Miyoshi2011). OpenMHD uses a divergence cleaning approach using a using a generalised Lagrange multiplier method (Dedner et al. Reference Dedner, Kemm, Kröner, Munz, Schnitzer and Wesenberg2002) and tracks $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {B}$ errors, which are shown to be ${\sim}10^{-6}$ at each time step in our CPG simulations.

The result of imposing the initial density radial dependence $r^{-1}$ in (4.1) and (4.2) is the simulated blow-by in frame 2 (17.5 $\unicode{x03BC}$ s) of figure 5. An axial sheath velocity asymmetry emerges with a higher velocity at smaller radii, despite the density profile remaining high at these smaller radii. It was found that the numerical factor of five times the Alfvén speed, $v_{z,0} = 5.0$ , in (4.4) provides the best match to the imaged dynamics of the plasma propagating across the domain, and this factor will be further discussed in § 5. While the simulated blow-by does not match precisely with the axial extent seen in plasma imaging, the position of the peak density in such blow-by profiles will be shown to be an accurate estimation (see below).

The initial azimuthal magnetic field $B_{\theta }$ of the plasma sheath follows the typical coaxial profile $B \sim \mu _{0} I / 2 \pi r$ , and where $B_{0} = 1$ is chosen for (4.5) for this run. We note that the magnetic field profile does not include a circuit current $i(t)$ term embedded within OpenMHD that is useful for more robust MHD modelling of coaxial accelerators (Datta et al. Reference Datta, Meier and Shumlak2024). This may become a future OpenMHD upgrade, but it out of the scope of this study. Since the goal of this paper is to produce an initial, first approximation solution for an open source resistive MHD model of a CPG, the general solution of this magnetic field with $B_{0}$ of unity allows the plasma sheath to evolve to match the observed dynamics (see § 5). We note that increasing the $B_{0}$ order-of-magnitude value to mimic the ${\sim}$ kA-level currents destroyed the plasma sheath and introduced far more reflection off of the outer electrode than is seen to be physical.

In the file $\mathsf {flux\_resistive.f90}$ , we embedded a plasma resistivity profile $\eta _{r, z}$ , where a background value $\eta _{0} = 1/1000$ has a localised peaked value distribution $\eta _{1} = 1/60$ (Zenitani & Miyoshi Reference Zenitani and Miyoshi2011) that is coincident to the plasma density and identical to the other $r^{\gamma }$ profiles in (4.6). Taking the step of imposing a localised resistivity profile $\eta _{1}$ is important for considerations of turbulence, the Hall effect, high atomic number $Z$ effects, electron kinetics and the emergence of a possible magnetic reconnection dynamics (Zenitani & Miyoshi Reference Zenitani and Miyoshi2011) in future studies that are outside the scope of this paper (Chesny et al. Reference Chesny, Hatfield, Cassibry and Moffett2025). However, it is important that the plasma resistivity is editable in OpenMHD to demonstrate its robustness as a resistive MHD architecture applicable to CPGs. Finally, we note that most of the profiles using (4.1) are thus proportional to $r^{-\gamma }$ and that the effects of blow-by, or preventing blow-by, can be observed and investigated in future studies.

4.2. Numerical results

The resistive OpenMHD modelling results are shown in figures 8 and 9. The information is presented in landscape format to show the largest format possible. From bottom-to-top in each image (and left-to-right viewing in landscape format), respectively, the columns are global shutter imaging from figure 5, density $\rho$ , pressure $p$ , azimuthal magnetic field $B_{\theta }$ . Vertically in landscape format, the model timings follow equal spacings in increments of $t^{\prime }=0.6$ .

While an asymmetry in the first ionisation frame is apparent in the radial direction (top and bottom) of the actual CPG (figure 8, top left), the initial density, pressure and magnetic field profiles (figure 8, top row, middle and right), show the $t=0$ simulation plasma profiles matching the general density behaviour. Note that the resistivity profile matches the same $r^{\gamma }$ gradient (where $\gamma = -1$ ), while the full sheath velocity profile is a constant in $r$ that follows the axial $\psi (r, z)$ distribution function. At $t=0.6$ , the plasma has propagated roughly halfway down the electrodes (figure 8, second row, second column) and the blow-by effect has manifested with the plasma at smaller radii running past the plasma at the outer electrode. This was predicted to be the blow-by effect due to the imprecise mass loading in the radial direction, where Cassibry et al. (Reference Cassibry, Thio and Wu2006) demonstrated that a $\gamma = -2.5$ would prevent blow-by and result in a radially uniform sheath velocity. This improper mass loading is a direct result of our CPG gas puff being located at the centre of the inter-electrode gap. If the gas puff inlets were located closer to the inner electrode radius, we anticipate that this blow-by effect would be mitigated. The authors stress that it is this initial behaviour of plasma initiation and the manifestation of blow-by during the rundown phase that was the goal of the resistive MHD modelling using OpenMHD. Understanding the plasma profiles and especially the internal magnetic field are the primary goals accomplished here. It is noted that a leading pressure wave begins to emerge at $t=0.6$ that precedes the denser blowby. This may be a numerical manifestation of the shock heated deflagration front (Cheng Reference Cheng1970).

The runover and pinch phases are captured by the simulated time stamps from $t=1.2$ and $t=1.8$ with the internal plasma at small radii becoming the densest portions into $t=2.4$ . We note that these axisymmetric simulations are only a cross-section of the plasma behaviour, while the global shutter imaging captures the full 3-D morphology as a projection onto a 2-D plane. Due to the axisymmetric nature of the simulation, an observer looking along the radial direction would see a 3-D projection of the radial plane all along the $Z$ axis. The concavity of the simulation would form the leading edge of the imaged plasma bulk. The leading deflagration front begins to precede the dense plasma and is most evident in the heated pressure shock front shown in the third column. As this plasma bulk propagates axially at $t=3.0$ and later, the more diffuse plasma density connecting the sheath back to the electrodes begins to fade, as it does in the imaging at $t=87.5$ and $t=105.0$ $\unicode{x03BC}$ s. This is indicative of the plasmoid formation, separation and propagation as the magnetic field is carried along with the bulk plasma density (figure 5). Here, the bulk plasma density becomes more concave as it propagates, and it elongates and diffuses.

An important implication from this modelling is that the magnetic field component $B_{\theta }$ embedded in the initial plasma sheath (fourth column) is dragged along with the plasmoid and both density and magnetic field peak values occur at the same locations. More insight into this implication will be provided in § 5 when these results are compared with the classic snowplough model, but this shows that imaging the location of bright density enhancements implies the location of the strongest azimuthal magnetic field (Stepanov et al. Reference Stepanov, Shumlak, McLean, Nelson, Claveau, Forbes, Weber and Zhang2020). When future TMR experiments are conducted and the plasma must be forced to interact with an external magnetic field in the fan-spine topology, the locations of the plasma’s perturbing magnetic field can be known to a high degree of fidelity.

The background density is $\rho _{b} = 10^{-2}$ of the plasma slug IC, and the pressure profiles (third column) demonstrate the interaction of the sheath with a substantial background (figures 8 and 9, right column). The CPG chamber is at a background pressure of 0.5 Torr, considered a medium vacuum level, and collisions with the neutral background are significant. This effect likely contributes to the lack of propagation of the plasma bulk in imaging frames at $t=87.5$ to $t=122.0$ $\unicode{x03BC}$ s.

The OpenMHD modelling results can be extrapolated to elucidate unknown experimental plasma parameters. If we consider dimensional quantities of typical CPGs of the order of $B_{0} = 0.5$ T (Black, Mayo & Caress Reference Black, Mayo and Caress1997) and the average molecular weight of air $A = 28.97$ (the working gas), an Alfvén travel time of $v_{z, 0}=v_{A}=5$ at $\tau _{A} = 54$ $\unicode{x03BC}$ s across a distance of $L=10$ cm occurs when the electron number density is $n_{e} = 3 \times 10^{21}$ m $^{-3}$ . This travel time is nearly identical to the empirical propagation extent of the plasma at the global shutter time frame of $52.5$ $\unicode{x03BC}$ s in figures 5 and 8 (landscape, first column, fourth row). The $L=10$ propagation is demonstrated in the simulation run of the density propagation at $t=1.8$ in figure 8 (landscape, second column, fourth row). This dimensionless MHD time is computed to be $\tau _{A}=2$ with the $|\rho |=|\boldsymbol {B}|=1$ condition at $L=10$ and $v_{z, 0}=v_{A}=5$ , which is very close to the simulated time of $t=1.8$ .