2. Analytic derivation

We start with the equations for ideal, isentropic, cylindrically symmetric MHD flow. For the density $\rho$, radial velocity $v_r$, azimuthal velocity $v_\phi$, pressure $p$, azimuthal magnetic field $B_\phi$, axial magnetic field $B_z$, we write the continuity, radial momentum, angular momentum, specific entropy and Faraday's law as

(2.1)\begin{gather} \frac{\partial \rho}{\partial t}+\frac{1}{r} \frac{\partial}{\partial r}(r \rho v_r)=0, \end{gather}

(2.2)\begin{gather} \frac{\partial(\rho v_r)}{\partial t}+\frac{1}{r} \frac{\partial}{\partial r}\left(r \rho v_r^{2}\right)+\frac{\partial p}{\partial r}+\frac{\partial}{\partial r}\left(\frac{B^{2}_\phi+B^{2}_z}{8 {\rm \pi}}\right)+\frac{B^{2}_\phi}{4 {\rm \pi}r}-\frac{\rho v^2_\phi}{r}=0, \end{gather}

(2.3)\begin{gather} \frac{\partial(\rho v_\phi)}{\partial t}+\frac{\partial}{\partial r}\left(\rho v_r v_\phi\right)+ 2\frac{\rho v_r v_\phi}{r}=0, \end{gather}

(2.4)\begin{gather} \frac{\partial p}{\partial t}+\frac{\gamma p}{r} \frac{\partial}{\partial r}(rv_r)+v_r \frac{\partial p}{\partial r}=0, \end{gather}

(2.5)\begin{gather} \frac{\partial B_\phi}{\partial t}+\frac{\partial}{\partial r}(v_r B_\phi)=0, \end{gather}

(2.6)\begin{gather} \frac{\partial B_z}{\partial t}+\frac1r\frac{\partial}{\partial r}(r v_r B_z)=0. \end{gather}

In the case of cylindrical symmetry, the radial component of the magnetic field, $B_r$ is zero due to the divergence-free condition for the magnetic field. We will assume that the z-velocity component is zero because, due to Galilean invariance it does not participate in the dynamics, but is rather passively advected by the flow. Later in the paper, we will study the stability of this problem with respect to the azimuthal perturbations, in which case most components will be non-zero. However, while doing this, we will keep the origin at rest, which is implied by our cylindrical coordinates.

We use the equations above to describe a self-similar imploding cylindrical plasma, which stagnates in a shock front expanding at a constant velocity. We assume that at the shock the entropy conservation should be replaced with a shock jump condition, as described below. We define the self-similar variable as

(2.7)\begin{equation} \eta \equiv \frac{r}{t} \Rightarrow r=\eta t, \end{equation}

we also assume without loss of generality that the outward propagating shock velocity is unity, i.e. the shock front is always located at $\eta =1$. We seek the self-similar MHD solution in the form

(2.8)$$\begin{gather} v_r(r,t)=\eta U(\eta), \end{gather}$$

(2.9)$$\begin{gather}v_\phi(r,t)=\eta W(\eta), \end{gather}$$

(2.10)$$\begin{gather}\rho(r, t)=t^{2 \chi} N(\eta), \end{gather}$$

(2.11)$$\begin{gather}p(r,t)=t^{2 \chi} P(\eta), \end{gather}$$

(2.12)$$\begin{gather}B_\phi(r,t)=t^{\chi} \sqrt{4{\rm \pi}} H_\phi(\eta), \end{gather}$$

(2.13)$$\begin{gather}B_z(r,t)=t^{\chi} \sqrt{4{\rm \pi}} H_z(\eta). \end{gather}$$

Here $\chi >0$ is a dimensionless parameter characterizing the self-similar scaling of physical quantities. The initial conditions for our problem will be

(2.14)\begin{equation} \left. \begin{gathered} v_r(r, t=0) =v_{r0} ; \quad v_\phi(r, t=0)=v_{\phi 0} ; \quad \rho(r, t=0)=\rho_{0} r^{2 \chi},\\ B_{\phi,z}(r, t=0) =B_{\{\phi,z\} 0} r^{\chi} ; \quad p(r,t=0)=0, \end{gathered} \right\} \end{equation}

we also restrict ourselves to zero initial pressure, as we seek to generalize the Noh problem. With finite density this corresponds to zero temperature of inflowing gas. These initial values are related to the limit $\eta \to \infty$ as

(2.15)\begin{equation} \left. \begin{gathered} v_{r0} =\lim_{\eta\to\infty} U(\eta) \eta,\quad v_{\phi 0}=\lim_{\eta\to\infty} W(\eta) \eta,\\ \rho_0 =\lim_{\eta\to\infty} N(\eta) \eta^{{-}2\chi},\quad B_{\{\phi,z\} 0}=\sqrt{4{\rm \pi}} \lim_{\eta\to\infty} H_{\phi,z}(\eta) \eta^{-\chi}. \end{gathered} \right\} \end{equation}

To transform variables $(r,t) \to (\eta, t)$, we use

(2.16a,b)\begin{equation} \frac{\partial}{\partial r} \rightarrow \frac{1}{t} \frac{\partial}{\partial \eta}, \quad \frac{\partial}{\partial t} \rightarrow \frac{\partial}{\partial t}-\frac{\eta}{t} \frac{\partial}{\partial \eta}. \end{equation}

Then, in the MHD equations the time derivative drops out and we arrive at

(2.17)$$\begin{gather} (U-1) \frac{{\rm d} \ln N}{{\rm d} \ln \eta}+\frac{{\rm d} U}{{\rm d} \ln \eta}+2 U+2 \chi=0, \end{gather}$$

(2.18)$$\begin{gather}(U-1)\left[\frac{{\rm d} U}{{\rm d} \ln \eta}+U\right]-W^2+\frac{1}{\eta^{2} N}\left[\frac{{\rm d} P}{{\rm d} \ln \eta}+\frac{H_\phi}{\eta} \frac{{\rm d}(\eta H_\phi)}{{\rm d} \ln \eta}+\frac12\frac{{\rm d}H_z^2}{{\rm d} \ln \eta}\right]=0, \end{gather}$$

(2.19)$$\begin{gather}(U-1) \frac{{\rm d} \ln W}{{\rm d} \ln \eta}+2U-1=0, \end{gather}$$

(2.20)$$\begin{gather}(U-1) \frac{{\rm d} \ln P}{{\rm d} \ln \eta}+\gamma \frac{{\rm d} U}{{\rm d} \ln \eta}+2 \gamma U+2 \chi=0, \end{gather}$$

(2.21)$$\begin{gather}(U-1) \frac{{\rm d} \ln H_\phi}{{\rm d} \ln \eta}+\frac{{\rm d} U}{{\rm d} \ln \eta}+U+\chi=0, \end{gather}$$

(2.22)$$\begin{gather}(U-1) \frac{{\rm d} \ln H_z}{{\rm d} \ln \eta}+\frac{{\rm d} U}{{\rm d} \ln \eta}+2U+\chi=0. \end{gather}$$

Note the following symmetry of (2.17)–(2.22): they are invariant under arbitrary scaling transforms of $\eta$, $P$, $N$, $H_\phi$ and $H_z$, provided that $P/\eta ^2N$, $H_\phi ^2/\eta ^2N$, $H_z^2/\eta ^2N$ stay constant. Physically this means rescaling units for the pressure, density and magnetic field while keeping sonic and Alfvén speed the same. This symmetry can be utilized to simplify the system (2.17)–(2.22) further by introducing new variables

(2.23a,b)\begin{equation} S=\frac{\gamma P}{\eta^{2} N(1-U)} \quad {\rm and} \quad A_{\{\phi,z\}}=\frac{H^{2}_{\{\phi,z\}}}{\eta^{2} N(1-U)}, \end{equation}

as done in Liberman & Velikovich (Reference Liberman and Velikovich1986), Felber et al. (Reference Felber, Malley, Wessel, Matzen, Palmer, Spielman, Liberman and Velikovich1988) and Venkatesan & Choe (Reference Venkatesan and Choe1992). We define Mach numbers as

(2.24a,b)\begin{equation} M_{s}^{2}=\frac{(v_r-\eta)^{2}}{c_{s}^{2}}=\frac{1-U}{S}, \quad M_{A\{\phi,z\}}^{2}=\frac{(v_r-\eta)^{2}}{v_{A\{\phi,z\}}^{2}}=\frac{1-U}{A_{\{\phi,z\}}}, \end{equation}

where $c_{s}=\sqrt {\gamma p/ \rho }$ and $v_{A_{\{\phi,z\}}}=B_{\{\phi,z\}} / \sqrt {4 {\rm \pi}\rho }$ are the local values of the speed of sound and Alfvén speed, $\eta$ is the ‘phase’ velocity of the self-similar profile. These are the local Mach numbers of the radial flow velocity relative to the surfaces of constant $\eta$, i.e. $v_r-\eta =\eta (U(\eta )-1)$. At $\eta =1$, the Mach numbers (2.24a,b) characterize the flow velocity with respect to the shock front. It is also useful to introduce the fast magnetosonic Mach number as

(2.25)\begin{equation} M_F^{2}=\frac{1-U}{S+A_\phi+A_z}. \end{equation}

Introducing new variables, we are left with the following five equations:

(2.26)$$\begin{gather} (U+S+A_\phi+A_z-1) \frac{{\rm d} U}{{\rm d} \ln \eta}+U(U-1)-W^2 \nonumber\\ \hskip6pc +\,2 S\left(U+\frac{\chi}{\gamma}\right)+(\chi+1)A_\phi+(\chi+2U)A_z=0, \end{gather}$$

(2.27)$$\begin{gather}(U-1) \frac{{\rm d} \ln W}{{\rm d} \ln \eta}+2U-1=0, \end{gather}$$

(2.28)$$\begin{gather}(U-1) \frac{{\rm d} \ln S}{{\rm d} \ln \eta}+\gamma \frac{{\rm d} U}{{\rm d} \ln \eta}+2(\gamma U-1)=0, \end{gather}$$

(2.29)$$\begin{gather}(U-1) \frac{{\rm d} \ln A_\phi}{{\rm d} \ln \eta}+2 \frac{{\rm d} U}{{\rm d} \ln \eta}+2U-2=0, \end{gather}$$

(2.30)$$\begin{gather}(U-1) \frac{{\rm d} \ln A_z}{{\rm d} \ln \eta}+2 \frac{{\rm d} U}{{\rm d} \ln \eta}+4U-2=0. \end{gather}$$

By combining (2.27)–(2.30) and (2.17) it is possible to write several explicit integrals of our system, namely

(2.31)$$\begin{gather} S^{\chi+1}N^{1-\gamma}\left[\eta^2(1-U)\right]^{\gamma\chi+1}={\rm const.}, \end{gather}$$

(2.32)$$\begin{gather}A_\phi(1-U)^2\eta^2=v^2_{A\phi 0}, \end{gather}$$

(2.33)$$\begin{gather}A_z^{\chi+1}N^{{-}1}\left[\eta^2(1-U)\right]^{2\chi+1}=v^{2(\chi+1)}_{Az0}\rho_0^{{-}1}, \end{gather}$$

(2.34)$$\begin{gather}W^2 A_z^{{-}1}(1-U)^{{-}2}=(v_{\phi 0}/v_{Az0})^2, \end{gather}$$

where $v_{A\{\phi,z\}0}=v_{A\{\phi,z\}}(t=0)$. These correspond to the frozen-in condition for specific entropy, $\phi$ and $z$ components of the magnetic field and conservation of angular momentum. Note that the first integral is broken at the shock, being zero for $\eta >1$ and non-zero for $\eta <1$ because the specific entropy does not conserve in shocks. Its value in the post-shock fluid can be obtained from the shock jump condition that we describe in the next section. Numerically, it is often easier to solve the system (2.17) and (2.26)–(2.30) directly, without involving integrals.

We look for regular solutions which have $U\leq 0$, so $(U-1)$ is always non-zero. The prefix $U+S+A_\phi +A_z-1$ in the (2.26) can be represented as $(U-1)(1-1/M_F^2)$. The solution becomes singular at the fast magnetosonic point where the local velocity with respect to the shock is equal to the group velocity of the fast magnetosonic mode. At large $\eta$ the fast magnetosonic Mach number $M_F$ is larger than unity, while after crossing the fast shock it becomes smaller than unity. A regular solution is obtained when the prefix in the (2.26) changes its sign only at the shock. (The solution can, in principle, pass the fast magnetosonic point regularly, if we require that the last four terms of the (2.26) sum to zero at the fast magnetosonic point, but this will require one extra constraint imposed on the initial conditions. We will not seek these special solutions and they will be considered elsewhere.)

3. Evaluating self-similar non-rotating solutions

We solve the main system (2.26)–(2.30) in a natural fashion, starting with initial conditions (2.14) that impose boundary conditions for our dynamic variables at infinity (2.15). Numerically, we start with a very large value of $\eta$ and integrate the system (2.26)–(2.30) towards $\eta =1$ where we assume the shock is located. We do so without a lack of generality, because if the shock is located at a different $\eta$, i.e. the shock has a different velocity, we can rescale all quantities to renormalize back to the shock speed of unity. For example, we can keep $N$ constant and rescale initial $A_{\phi,z}$ by a shock speed squared. At large values of $\eta$ our numerical solution is consistent with the analytical asymptotic formulae presented in Appendix B.

At the shock we apply the standard MHD shock jump conditions (Landau & Lifshitz Reference Landau and Lifshitz1960; Liberman & Velikovich Reference Liberman and Velikovich1986),

(3.1)$$\begin{gather} U \to 1-\mu(1-U), \end{gather}$$

(3.2)$$\begin{gather}A_{\{\phi,z\}} \to A_{\{\phi,z\}}/\mu^2, \end{gather}$$

(3.3)$$\begin{gather}S \to (1-U)/M_S^2, \end{gather}$$

where the density compression parameter $\mu =\rho _1/\rho _2$ and $M_S$ are determined by the solution of a quadratic equation for $\mu$,

(3.4)$$\begin{gather} \mu=\frac{\gamma_2+1/M_{A}^2+\sqrt{(\gamma_2+1/M_{A}^2)^2+4(1-2\gamma_2)(2-\gamma_2)/M_{A}^2}}{4-2\gamma_2}, \end{gather}$$

(3.5)$$\begin{gather}\frac1{M_S^2}=\gamma\left(1-\mu-\frac{1/\mu^2-1}{2M_{A}^2}\right), \end{gather}$$

where $\gamma _2=1-1/\gamma$ and $M_A^{-2}=M_{A\phi }^{-2}+M_{Az}^{-2}$. Keep in mind that $M_S^2$ is different from $M_s^2$ in (2.24a,b) in that the former is the ratio of pre-shocked $1-U$ to post-shock $S$. After the shock, we continue to integrate the system (2.26)–(2.30) keeping an eye on the sign of $U$ and the sign of $1-M_F$. We want to obtain the regular solution that stagnates on axis, i.e. have $\eta U \to 0$ as $\eta \to 0$. In Appendix A we analytically derive boundary conditions for $U$ at the origin, $U_0$. For $\gamma >2$, $U_0=-\chi /\gamma$, for $\gamma <2$ and the non-rotating case, $U_0=-\chi /2$. For the rotating case, the expression is more complicated and provided in Appendix A.

We can always obtain the solution that stagnates at larger $\eta$ by choosing sufficiently small $-v_0$. At the same time, if $-v_0$ is too large, we have the solution that either crosses the fast magnetosonic point or has $U-U_0<0$ at $\eta =0$. Between these two extreme cases is our desired solution, which has $U=U_0$ at the origin.

Numerically, we can not integrate down to $\eta =0$, so we choose some small but finite $\eta _{min}$ to stop integration. We seek $v_0$ that minimizes $U-U_0$ at $\eta _{min}$ by the bisection method, assuming that crossing the fast magnetosonic point means that $-v_0$ is too large while stagnating before $\eta _{min}$ means that $-v_0$ is too small. Once we obtained the converged solution, at small values of $\eta$, our numerical solution is consistent with the analytical asymptotic formulae presented in Appendix A.

Thus, we obtained a self-similar solution that has two parameters, $\chi$ and $\gamma$, and two initial conditions for $A_{\{\phi,z\}}$ determined from (2.15). The solution for velocity and magnetic field expressed in physical units, have two more degrees of freedom – we can rescale the flow and the shock velocity to an arbitrary value, which will rescale pressure by the same factor squared and we can rescale density and the magnetic field while keeping $B^2/\rho$ constant. The physical solution of magnetized Noh, therefore, may be classified by six parameters: $\chi$, $\gamma$, shock velocity, initial density and the two components of the initial magnetic field. We can also obtain solutions with finite initial pressure, but will not analyse them here. In Appendix A we derive asymptotic behaviour of density, pressure, $B_\phi$ and $B_z$ on axis. For the most physically relevant case of $\gamma <2$, we find that density, pressure and $B_\phi$ go to zero, while $B_z$ goes to a constant. The azimuthal current density $j_\phi$ that creates $B_z$ is regular if $\gamma <3/2$ and $\chi >2/(3-2\gamma )$. Otherwise, $j_\phi$ is singular, and since $B_z$ is energetically a dominant term on axis, this will create sizable dissipation in non-ideal plasmas.

We found it is convenient to classify Mag Noh problems without rotation by four dimensionless numbers: the ratio of the initial axial to the initial azimuthal field $r_b$; initial Alfvén speed in units of the shock speed $M_a=v_A/V_s$ and $\chi$ and $\gamma$. The rotation adds another dimensionless parameter, which is a $v_{\phi 0}/v_{r0}$ ratio.

For a given set of parameters, we calculated the self-similar solution in a manner described above with a simple Euler ordinary differential equation solver, starting at $\eta _{max}=10^5$ and ending at $\eta _{min}=10^{-5}$. Typically, $2\times 10^5$ steps was enough to produce an accurate solution. We used the Athena code (Stone et al. Reference Stone, Tomida, White and Felker2020; Athena++ development team 2021) in cylindrical geometry to verify that our solutions are realized in numerical calculation. Athena is a finite volume Godunov code and we used its one-dimensional solver with first-order spatial reconstruction. (Higher-order reconstruction schemes often produce oscillations next to the shock. In this work we did not focus on the fast convergence of the numerical solution and, for the sake of simplicity, we used first order. Our publicly released code can be trivially modified to test higher-order convergence.) In Athena simulations we typically assume a shock speed of $10^7\ \textrm {cm}\ \textrm {s}^{-1}$, domain size of 8 cm with 4000 grid points and we evolve the solution to 30 ns, so that the shock is located at 0.3 cm. We use initial conditions (2.14), reflecting boundary at the origin and inflow boundary at the outer radius. The inflow boundary is not the correct boundary for our problem, but we assume that the domain size of 8 cm is large enough so that the perturbation set by the boundary does not propagate into the useful solution range 0–0.6 cm in 30 ns. Larger domain sizes may be needed for a large inflow velocity. Simulation parameters are tunable and we release our code to the public for an easy adoption by the community; see § 9. Figure 2 presents two example solutions, one of them having a regular and the other singular $B_z$ on axis. We also release this verification code to the public (see data release). We scanned our four-dimensional parameter space to determine the initial inflow speed normalized by the shock speed $v_0/V_s$ as a function of the four main parameters: $\chi, \gamma, M_a$ and $B_z/B_\phi$. The results are presented in figure 3. The strongest dependence is on $\gamma$, with $\gamma$ approaching unity the compressibility increases and the shock Mach number can be very high. The dependence on the ratio of the magnetic components is rather weak, similarly, the index of similarity only moderately affects the strength of the shock. The dependence on the Alfvén Mach number upstream is interesting in that it reveals a solution that has $v_0=0$. This is a rather special case of magnetic Noh because unlike classic Noh it does not have a flow hitting the wall, the initial condition is stationary.

Figure 2.

Two example solutions of magnetized Noh (red) verified with Athena++ in a one-dimensional cylindrical simulation (green circles). Both solutions initially have $B_z/B_\phi =1$. Top solution: $\chi =2.6, \gamma =1.1, v_A/V_s=0.84$ and $v_0/V_s=20$ have regular $j_\phi$ on axis; bottom solution: $\chi =1.5, \gamma =5/3, v_A/V_s=0.52$ and $v_0/V_s=6.1$ have singular $j_\phi$ on axis. Note that Athena++ is an MHD code that solves the initial value problem and is not in any way aware of either the solution self-similarity or the expected shock speed.

Figure 3.

Dependence of $-v_0/V_s$ vs various parameters of the Mag Noh problem, in all but the last plot $B_z/B_\phi =1$. Also note that in the upper left the $v_0$ axis is logarithmic.

4. Stationary initial conditions

Let us look at this interesting special case with zero initial velocity in more detail. This case deserves closer scrutiny because it is closer to practical applications such as Z-pinches, where the gas, initially at rest, is driven towards the axis and is shocked and heated. As far as code verifications are concerned, stationary initial conditions are easier to set up in some codes, e.g. the moving mesh codes.

In the case of fluid at rest initially, the flow develops due to the magnetic tension. One may expect that it would take some time to accelerate the flow before it develops a shock, this is not the case with our self-similar solution, however. In our solution the shock is present for all times $t>0$, this is perhaps because density is very small near the origin and it takes very small tension to create a shock. In order to obtain these solutions, we modified the solver described above. Now we do not change the initial velocity to find the stagnating solution, but rather we change $B_\phi$. For sufficiently small $B_\phi$, we expect the solution to stagnate before the origin, while at sufficiently large $B_\phi$ the solution either hits the origin with non-zero velocity or passes the fast magnetosonic point. Again, we seek $B_\phi$ that minimizes $|U-U_0|$ at $\eta _{min}$ by the bisection method. All solutions we found have a sizable $v_A/V_s$, between 2.0 and 4.8, which is not surprising since some minimum magnetic tension is needed to create a shock.

While studying the parameter space we found that solutions with different $\chi$ and $\gamma$ look very similar. Changing the $B_z/B_\phi$ ratio, however, makes quite a bit of difference. In figure 4 we show several solutions for $\chi =1.5$ and $\gamma =5/3$ and different field ratios. Generally, the lower $B_z$ results in a stronger shock with larger compression of the gas and both components of the field. The smaller the $B_z/B_\phi$ ratio is, the more singular $j_\phi$ becomes at the origin. In figure 4, bottom, we show the verification of the $B_z/B_\phi =0.3$ solution with Athena++. Note that close to the origin the gas pressure goes to zero and is replaced by the magnetic pressure of $B_z$. The field ratio $B_z/B_\phi$ is usually small in realistic Z-pinches and such solutions are potentially of practical interest since singular $j_\phi$ can provide extra heating near the low-density origin and create a hot spot. Such heating is, of course, not present in our solutions of ideal MHD.

Figure 4.

Solutions with zero initial velocity with $\chi =1.5$ and $\gamma =5/3$. Plots (a–d) show several solutions with different $B_z/B_\phi$ ratios. (e–i) Verification with Athena++ of the above $B_z/B_\phi =0.3$ solution with $v_A/V_s=3.45$, singular $j_\phi$ on axis, self-similar (red) overplotted on Athena++ (green circles). Note that the top plot is log-scale, while the bottom plot is linear scale.

5. Rotating solutions

We can obtain rotating solutions as well, in a manner similar to described above. For convenience, we introduce the rotation parameter, which is the ratio of initial rotation velocity to radial velocity. Rotating initial conditions are unusual in that we have to set rotating velocity to constant, even at the origin, which produces infinite vorticity. Our self-similar solver, however, does not have an issue with that, because initial conditions are set at infinity. It is, however, very interesting how numerical codes will handle such a situation. Figure 5 shows verification of our self-similar solution with the solution produced by Athena++ and the agreement is pretty good. One feature of the rotating solution is a singularity of rotation velocity at the origin, which is a consequence of the conservation of angular momentum. The other interesting feature is that adding rotation makes $B_z$ go to zero on axis. This is easy to understand as the radial electric field $E_r$, which is absent in non-rotating solutions is now present and is equal to the product of $B_z$ and $v_\phi$. Thus, $B_z v_\phi$ must go to zero at the origin and since $v_\phi$ diverges due to conservation of angular momentum, $B_z$ must go to zero. In Appendix A we derived asymptotics of our solution close to the origin and found that stagnating solutions, shocked or not shocked, can only be obtained if $v_{\phi 0} < \chi ^{1/2} v_{Az0}$.

Figure 5.

Rotating solution produced by self-similar solver (red) overplotted on the solution produced by Athena++ (green dots).

6. Empirical findings of stability

Developing an analytical theory of stability of our self-similar Mag Noh problem requires significant time and effort, so for the time being we employed Athena++ in cylindrical geometry to investigate stability in an empirical fashion. We modified boundary conditions in such a way as to perform a simulation with five different domains in azimuthal angle $\phi$, each domain covering $2{\rm \pi}$ in angle. In the first domain, we did not introduce perturbations, while four other domains were seeded with a sinusoidal perturbation along $\phi$, corresponding to mode $m=\{1,2,4,8\}$ in one series of simulations and $m=\{3,5,6,10\}$ in the other. In the first series, the resolution was 64 grid points in each domain, and in the second, 60 points. We chose the resolution to be a common multiple of all $m$. Since we did not know the $m$ eigenmodes dependence in the radial direction, we have chosen the perturbation to be confined within the shocked region of the solution, for this reason, we introduced the perturbation, which was $10^{-3}$ in amplitude and present only in $B_\phi$ and $B_r$ (the latter to preserve the divergence-free condition) starting with some initial time $t_i$ after which the shocked region has been already well established. Had we chose to perturb multiple variables, this would complicate the matter since we do not know the phase between different vectors in eigenmode space.

Our first numerical experiments have demonstrated some unstable solutions, many of them grew very rapidly at the maximum $m$ possible (i.e. grid scale in $\phi$), seeded by the truncation error. To counteract this, we used a module to project the solution to the $\exp (-\textrm {i} m \phi )$ at each time step. Projection was also needed due to cross-contamination of modes. Even if we set the initial perturbation in energy to $10^{-6}$, due to rapid growth of instability it may not stay small to the end of the simulation. In this case it will produce other modes due to nonlinearity, which will grow on their own. The projection made sure that we always study a pure $\exp (-\textrm {i} m \phi )$ mode. The loss of perturbation energy due to the projection was small if the amplitude was small. If not, a sizable fraction of energy can be transferred to other modes and removed by the projection, which had a stabilizing effect. We should assume, therefore, that the growth to large amplitudes is not accurate; see also the discussion below.

Figure 6 demonstrates the ratio of the perturbation energy to the total energy of the shocked region as it evolves in time starting with $t_i=4\times 10^{-9}\ \textrm {s}$ till the end of the simulation at $t_e=3\times 10^{-8}\ \textrm {s}$. We had no way to guarantee that perturbation will stay small in all points; however, our projection module will eliminate any harmonics growing out of nonlinearity. Due to this the growth of the fastest growing modes could have been somewhat suppressed. Despite the fact that their evolution was not accurate at later times, we can still conclude that these quickly growing modes are unstable based on their evolution at earlier times. In figure 6 we have chosen two solutions from our parameter space, one of which we evaluate to be stable for all modes up to $m=10$ (except for 7 and 9, which we did not test), the other we judged to be unstable for all $m$ up to 10. To evaluate stability, we compared maximum relative perturbation at earlier times, $0< t<5\times 10^{-9}\ \textrm {s}$ and later times $6\times 10^{-9}\ \textrm {s}< t<3\times 10^{-8}\ \textrm {s}$. Needless to say, such evaluation should be viewed as preliminary, until the actual eigenmodes and eigenfrequencies are found.

Figure 6.

Relative perturbation energies for different azimuthal mode numbers $m$ for one stable and one unstable solution of magnetized Noh.

We scanned our parameter space on the $\gamma,\chi$ plane, looking for stability. Figure 7 shows minimum unstable $m$, if it exists. Primarily, self-similar solutions with low $\gamma$, high $\chi$ and high inflow velocity tend to be unstable. We also checked a few solutions with zero inflow velocity and found all of them stable.

Figure 7.

Minimum unstable $m$ mode in the $\gamma,\chi$ place with $M_a=0.52$ and $B_z/B_\phi =1$.

We did not study stability analytically and did not find exact eigenmodes of azimuthal perturbation; however, the possibility that we are seeing numerical rather than physical instability is rather remote. First of all, we are using Athena in cylindrical coordinates, so that the azimuthal perturbation is exactly along the coordinate axis. Numerical instability to azimuthal perturbations would require Athena to generate rather large fluxes of conservative quantities in the $\phi$ direction. Given that perturbations in the $\phi$ direction are small and fluxes are generated based on the solution of the Riemann problem, this is unlikely. Secondly, conservative codes such as Athena are less likely to produce numerical instabilities. Thirdly, in many cases we see that the instability grows in a very well resolved $m=1$ mode, while numerical instabilities are usually generated on a grid scale (large $m$). Fourthly, there is a clear dependence of the instability on physical parameters, such as $\gamma$, as shown in figure 7. There is no reason to find such dependence in numerical codes, especially those based on the solution of the Riemann problem.

7. Two-dimensional Cartesian simulations

In all the simulations above we used Athena in cylindrical geometry with first-order spatial reconstruction. Sometimes the physics of the problem necessitates using Cartesian numerical simulations. For example, when an instability becomes nonlinear and transitions to turbulence, we prefer to study such turbulence on a uniform grid to avoid numerical effects associated with the coordinate system's singular point. Despite extra errors, it is interesting to see how numerical codes in the Cartesian case will deal with our self-similar solutions, including the unstable solutions. Cartesian coordinates will introduce discretization errors, C4-symmetric or periodic with ${\rm \pi} /2$ in azimuthal angle, i.e. harmonics of $m=4$, in addition, there may be truncation errors, not necessarily symmetric. Figure 8 shows $\rho$ and $B_\phi$ on a 2-D plain obtained in Athena++ simulation with different orders of spatial reconstruction. We initiated the problem with parameters $\chi =2.6, \gamma =1.1, v_A/V_s=0.84, v_0/V_s=20$, zero rotation. Earlier in this study we empirically found this solution unstable to all $m$ modes up to 10. We see from the figure that the instability seeded by the C4-symmetric discretization error is dominant with the first and second examples, but in the third, the instability was perhaps seeded by the truncation error as well. This is consistent with our cylindrical simulation where the discretization error did not produce azimuthal perturbations, nevertheless, the instability grew from the truncation error. In example four in figure 8 we used a more diffusive solver, the Harten–Lax-van Leer–Einfeldt (HLLE) solver, and this solver produced broadly sensible, but inaccurate solutions, especially for $B_\phi$. Finally, we found that the stabilizing code with viscosity and magnetic diffusivity and going to higher resolution provides accurate solutions. To test this, we performed a convergence study with diffusivities, scaling as grid scale, and observed convergence to a self-similar solution, albeit slow. On the bottom of figure 8 we show a similar simulation where we initially seeded $m=7$ perturbation in the inflow and observed that it is dominant almost everywhere, except near the axis where the remaining $m=4$ perturbation is still noticeable. Our interest in simulating accurate unstable solutions is primarily to determine if one can obtain a perturbed solution with the physically seeded perturbation. Collapsing Z-pinches are often unstable, however, one often encounters numerical studies where the instability is seeded by numerical error. In this case the post-instability evolution and transition to turbulence is code-dependent and, therefore, meaningless. In order to properly study disruption due to instability one has to seed perturbations explicitly and provide an accurate solution.

Figure 8.

Two-dimensional Cartesian simulations of the same problem as presented on top of figure 2. Plots (a–d) show density from four simulations using different numerical schemes: HLLD with first-, second- and third-order spatial reconstructions and HLLE with a second-order spatial reconstruction. Plots (e–h) show $B_\phi$ for the same schemes as (a–d). Plots (i–k) show higher resolution HLLE second-order simulations, stabilized with grid-scale viscosity and resistivity. In addition, an $m=7$ perturbation with a $10^{-2}$ amplitude was introduced in the $B_\phi$ and $B_r$ in the inflow.

9. Data release

The code needed to reproduce figures in this paper is available at https://github.com/beresnyak/magnoh.

This release includes a code to produce self-similar solutions, to estimate accuracy, to verify solutions with Athena++, also it includes a code for Athena++ to study stability to azimuthal perturbations.