Hybrid MHD/PIC simulation of a deuterium gas puff z pinch

We present the hybrid MHD/PIC simulations of the time evolution of a deuterium gas puff z pinch. Recent experiments with 3-MA current pinches [6] made in a novel configuration have shown that the neutron yields can reach 3.6×1012. It was shown that the observed neutron spectra could be explained by a suprathermal distribution of deuterons with a power law fall off in the ion energy distribution function at large energy. In order to perform the numerical simulation of gas puff z pinch a new hybrid model was developed. The described hybrid model treats the electrons as a massless fluid and ions as macroparticles. The macroparticle dynamic is calculated with the use of PIC method. Ion-ion Coulomb collision is considered with the use of MC method. In the model simulation, in the configuration close to described in [6], it was obtained the neutron yields up to 1.2×1012. Most neutrons are not thermonuclear. This level of the neutron yield is reached only when a strongly nonuniform neck-like constriction of z-pinch plasma occurs. In this case, the obtained deuteron spectra (with energy up to several tens MeV) have suprathermal high energy tail. These simulations demonstrate the utility of the developed hybrid model for the z-pinch simulation.


Introduction
Z-pinches are known as efficient sources of x-ray radiation, fast electrons, ions and fast neutrons [1]. It is believed that as a neutron source the Z-pinch could be useful in radiation material science, radiobiology, nuclear medicine, recycling of nuclear waste and for a fusion-fission hybrid reactor. Significant efforts have been made to obtain neutron numbers sufficient for these applications.
Traditionally, the study of plasma liners is in line with the problems of inertial fusion [2][3][4]. Plasma liners along with the classic pinches and plasma foci can serve as sources of intense neutron radiation with a neutron yield in the range 10 9 -4×10 13 neutrons per pulse [5]. The highest neutron yield for pinches (~ 4×10 13 ) was received at the Z-machine (SNL, USA) in experiments on the compression of deuterium liners [5].
Recent experiments [6,7] made in a novel configuration have shown that the neutron yields can reach 3.6×10 12 at current about of 3 MA. Earlier, this neutron yield could be surpassed only on the 15-MA pinch. Energy neutron spectrum obtained in experiments [6,7] has a high energy non-Maxwellian tail. Possible reasons of this high energy tail origin are discussed, for example, in [8]. It was shown that the observed neutron spectra could be explained by a suprathermal distribution of deuterons with a power law fall off in the ion energy distribution function at large energy. These high energy deuterons with energy up to 40 MeV (highest achievement for Z-pinches experiments) were detected in [6,7]. The mechanism of generation of such deuterons is still not clear [8]. Computer simulation of Z-pinch implosion is a method to research the mechanism. MHD approach in this case is hardly applicable because the essential part of the deuterons has a suprathermal distribution. Various hybrid [9] and kinetic [10] models were developed to simulate different stages of Z-pinch dynamics. Here we describe a relatively new hybrid model, which has certain advantage in comparison with previously developed hybrid models, but is still simpler than the kinetic model [10].

Model description
A two-dimensional hybrid quasineutral model was developed to study the dynamic of the deuterium gas puff z pinch. According to the hybrid approach the electron subsystem is described as an inertialess fluid, but ions are described as macroparticles. The developed model is analogous to the model [11,12], but the described model is modified to account for the deuterium ionization and neutron production due to DD fusion. The governing equations in case of the 2D axial symmetrical geometry are the following: where V i -velocity of ion (of type i) macroparticles, r i -ion macroparticles radius vector, S-shape function of macroparticle (bilinear interpolation was used), H c -form factor of cell, N c -number of macroparticles in cell, n i -ion density, n e -electron density, u i -ion drift velocity, u e -the electron   conductivity, B-magnetic field (only  component exists in a given geometry), E-electric field, P e -electron pressure, -thermal conductivity, W r -energy loss due to bremsstrahlung radiation. The change of the ionic chaotic energy as a result of collisions with electrons is taken into account using the method described in [13]. The macroparticles undergo ion-ion Coulomb collisions. These collisions have been considered with the help of DSMC-like method developed in [14]. Processes of deuterium ionisation and DD fusion have been considered with the help of direct simulation Monte Carlo (DSMC) method. The geometric configuration of the computational domain has been chosen in order to match the geometry of the experiment [6,7]. Figure 1 shows sketch of the experimental setup and the computational domain. It is seen that the computational domain contains only the inner gas shell. The outer plasma shell is far away from the domain. The model does not consider the motion of the plasma shell. Our calculation starts at the instant when the plasma shell reaches the inner gas shell. The shell parameters at this point are the additional parameters of the problem. We varied velocity and linear mass distribution of the outer plasma shell in order to get the neutron yield as close to the results [6,7] as possible. Figure 2 shows an example of the initial density distribution. Deuterium density (n D ) in the inner shell is of about 10 18 cm -3 (linear mass 120 g/cm). Density of hydrogen plasma (n i ) in the outer shell is of about 10 19 cm -3 (linear mass 5 g/cm). Velocity of the incident plasma shell 610 7 cm/s.

Results and discussion
Possible mechanisms of neutron generation in DD fusion are the following [8]: thermonuclear mechanism (i.e. neutron production due to high ion temperature), ion acceleration due to the axial plasma outflow from a constriction, generation of fast ion in the low dense plasma surrounding the pinch column, target mechanism [15], coulomb explosion [16]. The last two mechanisms are not considered in our model because of assumption of plasma quasineutrality. However, we believe that the contribution of these mechanisms in the neutron generation is insignificant because in the case of their implementation the neutron pulse should be very short (less than 1 ns), while in the experiments it can be tens of nanoseconds. The highest neutron yield reached in our simulation at the mentioned liner masses of the shells and current 2.8 MA is 1.2410 12 neutrons that is in two time less than the yield obtained in the experiment. This result was obtained only at a certain density distribution of the incident shell. The shell density distribution has a pronounced maximum at the center, as shown in Figure 2. If the density distribution of incident shell has this form then the further compression of the inner shell remains stable almost up to the axis, as it shown in Figure 3. Two plasma streams appeared in the shell in this case, one directed to cathode and other directed to anode. Evidently, these streams stabilise the compression. However, the inner shell loses up to 25 percent in mass due to these streams. If the outer shell has for example a flat density distribution then instabilities quickly develop and the compression effectiveness decreases. Neutron yield in this case is reduced by an order.
Due to the particular form of the shell ( Figure 3) the compression at the axis starts approximately at the center (along Z) producing a neck-like constriction (Figure 4). Following the first constriction the next one or two constrictions appears during of about 20 ns ( Figure 5). Due to this no uniform constriction a plasma outflow from the neck appears ( Figure 6). As it seen from the figure, the ions escaping from the neck looks like beams propagating through relatively quiet plasma. The source of the most energetic deuterons in this case is the cave (hot spot) at Z of about 1 cm (compare Figures 4  and 6). This cave can be treated as the beginning of the neck disruption. Low density plasma in this cave can be accelerated up to several MeV with the help of Laurence force.
It is clear that this picture cannot be obtained with the help of MHD approach. Appearance of the high energy ion beams can in principle explain the suprathermal high energy tail of neutron energy distribution (Figure 7).
Energy distributions of deuterons come to the electrodes within the circle with radius 0.5 cm are shown in Figure 8. It is seen that deuterons can reach energy up to 100 MeV. Note that deuterons with energy up 50 MeV were detected in experiments [6,7]. Authors of paper [6] have shown that the   measured suprathermal neutron energy distribution can be explained if one assume that deuterons have an isotropic distribution function with power law fall E -1.8 . As it seen in Figure 8 the power low of high energy deuterons obtained in our calculation is rather E -1 . However, this contradiction with [6] can be easily explained by strongly anisotropic energy distributions of deuterons in plasma in our case ( Figure 6). Let us discuss briefly a possible mechanism of the production of suprathermal ions and neutrons. Average ion temperature in our simulation is of about 2 keV. Thus, the thermonuclear mechanism makes a certain contribution in total neutron yield but it is not a source of suprathermal neutrons. Plasma acceleration due to the axial plasma outflow from the constriction can accelerate ion up to 100 keV. As it shown in [8] the ion distribution function in this case should be proportional to E -2 . The main input in the production of suprathermal ion in our simulation is contributed by the ions accelerated in low dense plasma surrounding the pinch column. In our case the low dense plasma appears because the ion Larmor radius becomes more than the size of hot spot, and ions due to the Larmor rotation can escape from the hot spot area. In contrast to MHD models ( [8] in example) where the described mechanism is hardly possible, our hybrid model is able to consider this mechanism.