Full-scale modeling and experimental study of a gas neutron tube with a Penning ion source

Spectrometric pulsed neutron gamma-ray logging is commonly used for geophysical studies of oil wells with casing.^1,2 Neutron generators with gas neutron tubes (hereinafter referred to as NT)^1,3,4 are used to produce a neutron pulse of a specific intensity, duration, and frequency.

A neutron tube is a miniature linear accelerator consisting of an ion source, a system of accelerating electrodes, and a target node, assembled in a sealed frame. Traditionally, one of the main types of ion sources used to create NT is a Penning ion generator (hereinafter referred to as PIG source).^1,5,6 Despite a long history of research on Penning ion sources,^7–12 experiments are still being carried out to optimize geometric parameters, such as the radius/length of the anode, distance between electrodes, cathode diameters, extraction hole diameter, and anti-cathode shape. Research on optimizing PIG source operating modes (selection of anode voltage, gas pressure, magnetic field, and mass-charge composition) is also extremely relevant.^13–22 In Refs. 23–29, the authors studied the main NT parameters experimentally, namely, the neutron flux as a function of discharge current and accelerating voltage, as well as the amplitude-time characteristics of current pulses. Various designs and structures of PIG sources were discussed in those works, the experimental results however do not enable to extrapolate research on a specific structure to the entire class of similar devices. There are also a number of works devoted to numerical modeling of plasma behavior in separate parts of NT. PIC-MCC (particle-in-cell with Monte–Carlo collisions) methods are used to model gas-discharge processes in PIG.^30–40 These models typically study the dynamic behavior electrons and ions, as well as the dependence of discharge/extract currents on magnetic field or discharge voltage. Other models focus on the motion of deuterium ions in the NT accelerating system, using SIMION 8.0,⁴¹ CST Studio Suite,^42,43 Comsol Multiphysics,^44,45 and SUMA⁴⁶ software. The ion beam size on the target is most often studied in relation to the electrode geometry and the accelerating voltage.^41–44 In Refs. 45 and 46 the influence of gas medium on the transport of ions from the PIG source to target is taken into account, and the neutron flux is estimated in Refs. 47–50. Besides, authors used analytical formulas to study the effect of the atomic-molecular and isotopic composition of ions on the neutron yield from NT. Rather than using values based on calculated or experimental data, they set the atomic-molecular and isotopic composition as a variable parameter.

Conflicting experimental data and incomplete numerical modeling results in those works mentioned above encourage us to perform a full cycle of experimental and research studies of NT, supported by computer modeling. Therefore, the purpose of this work is to provide not only a full numerical modeling of all NT physical processes, but also to study experimentally and determine the performance characteristics of NT with the PIG source. This will enable a comparison of the numerical modeling with the experimental results.

Figure 1 shows a simplified diagram of the gas NT with the PIG source. The PIG source consists of an anode and two cathodes placed in a longitudinal magnetic field. The shape of the anode is a cylinder with a potential typically ranging from $1$ to 10 kV. Next, there are the accelerating electrode and the neutron-producing target under a high negative potential (about 100 kV), which are galvanically isolated from the ion source using a high-voltage insulator.

The main advantage of NT with the PIG source is a simple control and power supply system, as well as stable operation at a low gas pressure. To operate the tube, current is sent through a filament of the gas storage made from hydride-forming material, which heats up and causes neutral gases such as D $2$⁠, T $2$⁠, and DT to desorb. The tube pressure typically ranges from $1$ to 10 mTorr. At the initial moment, a voltage pulse $U a$ of several kV is applied to the anode, causing free electrons to fall into a quadrupole electric and axial magnetic field. In crossed $E×B$ fields, these free electrons gain enough energy to ionize a neutral gas, resulting in an avalanche-like ignition of a gas discharge. The positive ions produced by impact ionization are then accelerated toward the cathode and anti-cathode. Due to the short residence time of ions within the discharge region, compared to electrons, a negative volumetric charge is formed. This reduces the ion emission rate, compensating for the potential drop. At this point, the discharge reaches a quasi-stationary equilibrium and the discharge current stabilizes. These physical processes establish an independent gas discharge in the volume.

The accelerating electrode potential ions from the gas discharge plasma to be then accelerated and focused on the neutron-producing target. A nuclear reaction occurs, when hydrogen (H) isotopes collide with D/T atoms in the target, releasing neutrons:¹ D + T $→ 4$He [3.5 MeV] + $n$ [14.1 MeV]. Interaction of ions with target atoms, due to ion-electron emission, generates a counter beam of electrons directed toward the PIG source. To suppress this electron current, a transverse magnetic field is created in the target node.

A Penning gas discharge ignites at a relatively low pressure (⁠ $1 − 10$ mTorr). Electrons are delayed in crossed $E×B$ fields and produce significant ionization of the neutral gas so that the plasma is self-sustaining. To model this physical system under consideration, the particle-in-cell method with Monte–Carlo collisions (PIC-MCC) is a good approach at such pressures.^51–54 

The PIC-MCC enables direct calculation of gas discharge combustion in the PIG source, as well as ion acceleration and focusing onto different NT elements. To create a full-scale numerical and theoretical method, it is needed to model the processes occurring in the target itself, such as target filling with H isotopes, D/T diffusion and thermal desorption during NT operation, NT elements sputtering, and nuclear reaction modeling. The NT full-scale modeling was split into the following stages due to different physical processes:

• Gas discharge combustion in the PIG source.

• Ions drift in the ion-optical system (IOS).

• H isotope diffusion and thermal desorption.

• Target sputtering.

• Nuclear reactions with the release of neutrons.

We used a 3D electrostatic PIC method based on structured rectangular grids and implemented in the TechX VSim software package⁵⁵ for numerical modeling within this study. Application of three-dimensional geometry is due to the occurrence of instability of plasma discharge combustion under certain conditions (when changes the mixture pressure, anode voltage, or magnetic field), as shown in Refs. 34 and 56. The Monte–Carlo collision method was applied to simulate kinetic processes in a gas-discharge plasma.

This study considers the characteristic dimensions of the Penning sources on the order of several cm, while the gas discharge volume is $1−10 cm 3$⁠, the magnetic field ranges from 100 to 1000 G, and the anode voltage is between $1 and5 kV$⁠. H isotopes were used as the neutral gas at pressures from 1 to 10 mTorr. The experimental electron temperature in miniature Penning gas discharges depends on discharge characteristics and gas type. In Refs. 57–59, the measured electron temperature for hydrogen/deuterium plasma was in the range from 2 to 12 eV. The maximum electron density in these experiments reached values up to $3.7× 10 15 m − 3$⁠.

The following plasma characteristics can be assessed:⁶⁰ Debye radius (⁠ $λ D$⁠), plasma electron frequency (⁠ $ω p e$⁠), cyclotron electron frequency (⁠ $ω c e$⁠). These values limit the grid size and time step for PIC-MCC methods^61–64 (for estimate were used $T e=1 eV$⁠, $n e=1× 10 16 m − 3$⁠, and $B=1000 G$⁠),

• Maximum cell size: $Δ x m a x<3.4× λ D=3.4×7433× T e [ eV ] n e [ m − 3 ]$ $∼0.25 mm$⁠.

• Maximum time step conditioned by the plasma frequency: $Δ t m a x< 0.2 ω p e= 0.2 5.64 × 10 4 n e [ cm − 3 ]∼35 ps$⁠.

• Maximum time step conditioned by the cyclotron frequency: $Δ t m a x< 0.2 ω c e= 0.2 1.76 × 10 7 × B [ G ]∼11 ps$⁠.

There are generally accepted criteria in literature sources⁶⁴ for the correct operation of Monte–Carlo methods, which guide the selection of the real particle number in a single macroparticle (⁠ $N P I M$⁠). Typically, $N P I M$ is selected so that the number of macroparticles in the Debye cube^64,65 falls within the range from 10 to 100. Then, we get the following constraints on $N P I M$⁠:

• at $n e= 10 16 m − 3$⁠: $N P I M∼1× 10 3−1× 10 4$⁠.

The grid size over the entire computational domain was uniform (without thickening) and made up $0.2 mm$⁠. The time step was set to 5 ps. The calculations accounted for three types of macroparticles: electrons (⁠ $e −$⁠), molecular ions (⁠ $D 2 +$⁠), and atomic ions (⁠ $D +$⁠). At the beginning of the numerical modeling, electrons and ions $D 2 +$ with a density of $10 14 m − 3$ were placed in the cylinder region of anode. The initial velocities of macroparticles were determined using the Maxwell distribution at a temperature of 300 K.

In Refs. 66–73, a kinetic model of plasma-chemical reactions in molecular hydrogen occurring in a Penning discharge plasma is presented, as well as the results of numerical modeling for the component gas-discharge plasma composition in a zero-dimensional approximation. The model includes over 100 different plasma chemical reactions and defines limiting reactions that reproduce the chemical kinetics of the Penning gas discharge with a good accuracy of about 5%. Based on the studies mentioned above, as well as on the fact that the kinetics of interaction depends more on the electronic structure than on the mass of particles (modeling of a gas discharge plasma in H isotope atmosphere), were identified by the following plasma-chemical reactions used in the modeling:

• ionization: $e −+ H 2→ H 2 ++2 e −$⁠.

• ionization: $e −+ H 2→ H ++ H+2 e −$⁠.

• scattering: $e −+ H 2→ e −+ H 2$⁠.

• ionization: $H 2 ++ H 2→ H 2 ++ H 2 ++ e −$⁠.

• scattering: $H 2 ++ H 2→ H 2 ++ H 2$⁠.

• charge exchange: $H 2 ++ H 2→ H 2+ H 2 +$⁠.

• ionization: $H ++ H 2→ H ++ H 2 ++ e −$⁠.

• scattering: $H ++ H 2→ H ++ H 2$⁠.

The primary electron source in a Penning gas discharge is electron impact ionization of the neutral gas. In order to accurately model the electron temperature, it is essential to calculate the scattering of electrons on the neutral gas with high precision. As demonstrated in Ref. 83, the primary contribution to the scattering of electrons on the neutral gas up to energies of $20−30 eV$ is elastic scattering. Inelastic scattering of electrons becomes a significant factor at energies above 100 eV, where the electron concentration is very low.

The Tech-X VSim software package has an ion-electron emission model for stainless steel, verified by a variety of experiments.⁸⁴ This particular model was used when modeling the PIG source for electrodes made of stainless steel and kovar. However, some cathode parts are made from molybdenum to prevent melting caused by the electron beam from the target. The ion-electron emission coefficient for Mo parts was⁸⁵  $γ=0.4$ for $D +$ and $γ=0.5$ for $D 2 +$⁠.

Once the processes in the PIG source enter the quasi-stationary mode, the NT geometry, including the PIG source, is modeled in a full scale. To do this, all charged particle and field distributions in the Penning discharge are stored at some point in time and used as a start point for a new calculation of the entire NT geometry.

In full NT modeling, the energy of charged particles can reach up to 100 keV (the voltage between the anode and the target), while the maximum energy of charged particles in the PIG source is limited by the voltage between the anode and the cathode (which is on the order of 5 keV). To preserve the characteristic paths of charged particles and the correct operation of the Monte Carlo method, the time step was reduced to 2 ps (the characteristic time for the electron to transit with the maximum distance energy corresponding to the cell size).

The modeling also considered secondary ion-electron emission from the target surface (i.e., the electron ejection process from structural elements under ion hit). An ejected electron moves toward a “neutral” ion source (the cathode/anti-cathode potential is 0 kV, the anode potential is $2−5 kV$⁠). As a result, an emission current of electrons is created. To suppress the electron current in the target node, a uniform magnetic field with a strength of $B=1000 G$ perpendicular to the NT axis was created when modeling. The residual gas effect on the ion beam travel in IOS was also taken into account.⁴⁴ The main objective for this stage of modeling is to obtain the ion current distribution for all NT structural elements, but first of all, to determine the ion current density on the target.

The knowledge of the ion current density distribution on the target allows us to proceed with the consideration of processes in the target itself. Figure 3 shows the NT target operation diagram. It is generally a thin film of a hydride-forming metal (such as Ti, Sc, Er, and others), several $mkm$ thick, sputtered onto a substrate that does not dissolve hydrogen.^1,86 Therefore, the substrate within this study is a cylinder with a diameter of several centimeters, whose single end is sputtered with a target material. In our case, Ti sputtered onto a Mo substrate was used as the target.

When H isotopes are bombarding the target node, some ions interact with the target D/T, leading to a nuclear reaction with the release of neutrons. However, most of those ions accumulate in Ti over time, leading to the formation of Ti hydride. H isotopes can diffuse and desorb from the target surface over time, when exposed to ion irradiation. There is no diffusion of H into the Mo substrate, since H is poorly soluble in Mo. These processes directly affect the neutron yield because the neutron flux depends on the amount of D and/or T retained in the target.

The diffusion and thermal desorption of H isotopes from $TiD x T y$ depend on the crystal structure of the hydride, the mechanisms of H diffusion, and the boundary conditions on the free surface. The highest stoichiometric H concentration is observed in $δ -TiH 2$ with a face-centered cubic (fcc) lattice, 2 H atoms per one Ti atom.^87,88

In Refs. 89–91, a model was developed that accurately calculates the diffusion and thermal desorption of H in Ti hydrides under the effect of ion beam. The main difference from classical models, representing the ordinary diffusion equation and the thermal desorption condition on a free surface,^90,91 is the consideration of the H sublattice. Hydrogen can take two positions in $δ$-Ti: tetrahedral (⁠ $t$⁠) and octahedral (⁠ $o$⁠) (for more details see Refs. 90 and 91). The main positions are tetrahedral. Hydrogen solution energy in $t$ positions is negative, and positive in $o$ positions (the hydrogen atom energy in the H $2$ molecule is assumed to be zero). In other words, it is energetically advantageous for H to be absorbed in $t$ positions and desorbed from $o$ positions.

A constant influx of $D +/ T +$ ions into the target results in the occupation of all $t$-positions in $δ$-Ti. As the number of these positions is limited, H isotopes start occupying energetically unfavorable $o$-positions. Desorption from the $o$- and $t$-positions occurs at different rates. Therefore, it is important to consider the H sublattice of the target material. In our calculations, we used the diffusion and thermal desorption model described in Refs. 89–91. The modeling was performed using methods of continuum mechanics in the OpenFOAM hydrodynamic software package.⁹² 

During the NT operation, the target is bombarded with $D +$ and $T +$ ions with an impact energy of 100 keV. The thickness of the target hydride layer is on the order of several $mkm$⁠. Long-term bombardment of the NT target with $D +/ T +$ ions causes surface sputtering, which erodes the target hydride layer and directly affects the neutron yield.

While bombarding the target, $D +/ T +$ ions form a collision cascade. Such cascades can return and reach the surface, transferring energy to near-surface atoms greater than the binding energy to the target surface. In this case, the near-surface atom is ejected from the target. When the target is thin (the average path length of the incident ions is comparable to the target thickness), then the collision cascades can reach the rear surface, and sputter from the rear surface of the target. The sputtering coefficient is the average number of atoms ejected from the target per incident ion and depends on the angle of incidence, the energy and mass of ions, as well as on the target atoms and their binding energy to the surface. For crystalline targets, the orientation of the crystallographic planes relative to the target surface is also relevant.

When irradiating a multi-component target with ions, the sputtering of one of the components can proceed at a higher rate. This is due to more efficient energy transfer to that component and/or lower binding energy to the surface. If the A component is sputtered more efficiently in the AB alloy, then during a long-term bombardment, the solid surface is enriched with the B component, thereby increasing the probability of B sputtering.

The sputtering coefficients of Ti hydride (TiDT) by deuterium and tritium ions were determined using the SRIM/TRIM software package.⁹³ The SRIM software package implements the Thompson analytical model for cascade sputtering of amorphous flat films.⁹⁴ A Ti target with a concentration of deuterium and tritium equal to TiDT with a 2  $mkm$ thickness was set in SRIM to determine the sputtering coefficients. This target underwent bombardment with $D +/ T +$ ions of different energies at an angle close to the normal (the deviation from the normal was $1 °− 5 °$⁠). The ion energy ranged from 5 to 100 keV. The modeling studied the dependence of Ti atom sputtering coefficients on the energy of incident $D +/ T +$ ions. The sputtering target rate was defined by the sputtering coefficient for Ti atoms, since it forms the hydride and determines the target thickness.

The deuterium–tritium mixture is desorbed when a getter with an equal D/T content is heated. Since molecules are mainly desorbed (atomic desorption is not energetically advantageous), the resulting composition of the NT neutral gas is as follows: DT $∼$ 50%, $D 2$ $∼$ 25%, $T 2$ $∼$ 25%. DT gas is predominant. Molecular ions make up the main part of the current supplied to the NT target (⁠ $∼$90% of the total current coming to the target).¹⁴ So, when two $DT +$ molecular ions with a 100 keV energy collide, they produce a different number of neutrons compared to the collision of $D 2 +$ and $T 2 +$ molecular ions with the same energy. Since the DT ion binding energy is on the order of several eV, and the kinetic energy is $∼100 keV$⁠, when such an ion hits the target, it immediately decays. Further nuclear reactions involve D and T ions separately from each other. Energy between D and T is redistributed by their masses.

The experimental measurements can be divided into three blocks: experiments related to the PIG source, ion-optical system study, and neutron flux measurement. The main objective of the PIG source experimental studies was to determine the optimal gas discharge combustion modes to achieve the required values of the discharge/extraction currents and maximum ion source energy efficiency. The ion-optical system study aimed to achieve the most uniform ion current distribution on the target.

This work presents experiments on NT parts that were already studied on various experimental installations.^{14,16,19,56,97,98} Below, we briefly summarize the results and links to individual studies on the PIG source and IOS parameters. These studies describe the experimental installations. This section focuses on describing the experimental installation used to measure the main NT parameters, including amplitude-time (ATC) and current-voltage characteristics (CVC) of discharge and extraction currents, parameters for the neutron yield and operating time.

Measurements were taken of the CV characteristics, energy distribution, and mass-charge ion composition emitted from the discharge under various PIG source combustion modes using the energy-mass analysis equipment.¹⁴ The results show a correlation between surges of the discharge current and the increase in potential sag (up to 50% of the anode voltage). The concentration of atomic hydrogen ions was determined for various discharge parameters.

The experimental installation¹⁶ was used to study operating modes of the PIG source when implementing both DC and AC power supply modes. Effects of the voltage amplitude on the PIG anode, pulse repetition frequency (PRF), and pulse duration in relation to PIG ATC and CVC were also studied. The study further identified the ion source operating modes and the gas pressure ranges, at which various combustion modes for the PIG discharge are implemented.

The experimental installation and technique for measuring direct current distribution over the target surface under a negative potential are described in Refs. 56 and 98. Also, Refs. 56 and 98 show the visualization of both the ion beam and the ion flux trace on the target surface, which was carried out by the luminous beam method. Examples of detected current distributions on the target surface are provided, depending on the geometric and physical NT IOS parameters.

The experimental installation used in this study is presented in Figs. 4 and 5. It includes a power supply and measurement system and a bath with transformer oil [Fig. 5(1)], where the NT under study is installed [Fig. 5(13)]. The power supply and measurement system are similar to those outlined in Refs. 16 and 19 and are shown in Fig. 4. The pulse PSU for the discharge generation in this experiment is assembled using the HVS-6-10 solid-state switches. This unit is based on a push–pull configuration for generating voltage pulses on across the load with partial discharge of $C 1$ capacitive energy storage device (K75-14: 0.1 mkF, 4 kV), whereby the duration of voltage pulses can be adjusted within wide limits. Time switching modes for the switches, including pulse repetition frequency, duration of control pulses and pauses between them, are set by an external two-channel Tektronix AFG 3022 pulse generator [Fig. 5(3)]. The capacitive energy storage device is charged to a maximum voltage of 4 kV using an external high-voltage DC power supply Spellman SL 10P300 [Fig. 5(4)], controlled remotely via an analog signal generated by an ICP CON ET7026 digital-to-analog converter (DAC) [Fig. 5(10)]. The storage device is charged through a series inductive-resistive circuit (⁠ $L 1=0.7H$⁠, $R 1=820Ω$⁠). The ICP CON ET7026 is powered using a laboratory DC power supply with an output voltage of $+12 V$ [Fig. 4(9)].

The energy storage device is connected to the PSU output circuits, when the VT1 switch triggers, and the front and the peak of the output voltage pulse are formed. The VT2 switch is triggered after turning off the VT1 voltage switch, providing conditions for the formation of a short tail of the voltage pulse. The switches operate sequentially with a short pause in their switching (⁠ $∼1$ mks) to prevent through currents from flowing across them. The wiring diagram in Fig. 4 uses current-limiting resistors, specifically $R 5=0.47 kΩ$ and $R 6=3.3 kΩ$⁠, to reduce peak current amplitudes when switching switches. The low-voltage energy storage device output is connected to the ground bus through a resistive shunt, including a $0.6 kΩ$ resistor (⁠ $R 3$⁠), a varistor (⁠ $R 4$⁠), and a two-anode protective diode (VD1). The shunt signal is used to measure pulse parameters of the NT ion source anode current. To measure the voltage pulse at the anode, a Tektronix P6015A high-voltage divider was used [Fig. 4(11)]. Pulse signals were recorded with a Tektronix DPO 3014 oscilloscope [Fig. 4(7)].

The integrated control circuits of the HVS-6-10 switches are powered by a laboratory DC power supply with an output voltage of +12 V [Fig. 4(8)]. The NT getter is powered by PSU2. Resistance of $R 2=90Ω$ is necessary to limit the getter filament current. The accelerating voltage is supplied using a high-voltage power supply source Spellman SL 1200 (PSU3) [Fig. 4(7)]. The extraction current is measured using a capacitance–resistance divider $C 2=660 pF$⁠, $R 7=0.6 kΩ$⁠, and $R 8=1 MΩ$⁠. The signal from the $R 7$ shunt (installed in the grounding bus of the energy storage device), comes to the Tektronix DPO 3014 oscilloscope input. The measurements are synchronized according to the control signal of the Tektronix AFG 3022 generator, which forms the leading voltage pulse edge.

The INPA (automated neutron flux density meter) detection unit was used to measure neutrons (the unit is shown in Figs. 4 and 6). The measured neutron flux density was then converted into a neutron flux, considering the 150 mm distance between the target and detector. INPA measurements were calibrated following the activation technique with the application of calibrated copper plates. The INPA error in measuring the neutron flux was 12%.

The 8-channel D-link switch [Fig. 5(5)] was used for the communication of peripheral devices with a PC, except for the INPA unit (⁠ $RS−232$⁠) and AKIP 1129 PSU (⁠ $USB$/virtual $RS$-232). Automation and installation control were developed in the Labview (National Instruments) environment. The control program implements a virtual proportional integral-differentiating regulator to create a control circuit with feedback on the discharge current, the extraction current and the neutron flux.

Experimental data on discharge combustion modes in the PIG source, as well as CVC dependences on pressure, magnetic field, etc., are presented in Refs. 14 and 16. The experiments were carried out with a pulsed supply of anode voltage at a frequency of 10 kHz, a duration of 30 mks, and an amplitude of 2 kV. The obtained results^14,16 were compared with those previously obtained in Refs. 16 and 19, where the source inside pressure was controlled on a completely identical deuterium NT model.

The Penning source^14,16 was used in the simulation. The anode voltage was 2 kV. Neutral gas pressure was constant. Deuterium was used as the neutral gas. Simulations were performed for pressures ranging from 1 to 8 mTorr. The calculations were performed for an isotropic magnetic field (⁠ $B=700 G$⁠) and for an anisotropic field^14,16 (permanent magnets were used). The weight of the macroparticles was 2500 and 5000 to compare the results.

To reduce the estimated discharge ignition time during modeling, an electrically neutral plasma was initially placed in the central region. The discharge stabilizes and reaches a stationary mode over time. The stationary mode is achieved when the number of macroparticles in the calculation does not change over 1–2 mks. Figure 7 shows the dependence of the full particles number in the calculation on time for $B=700 G$ and $P=6 mTorr$⁠. The time required to reach a stationary state ranges from 2 to 10 mks, with variations contingent on the pressure and magnetic field. The full number of macroparticles in the different calculations were about $1−2× 10 6$⁠.

The energy particle distribution in a Penning gas discharge does not adhere to the Maxwellian type. Figure 8 illustrates the distributions of electrons (⁠ $e −$⁠) and deuterium ions (⁠ $D 2 +$⁠) in a simulation for a magnetic field $B=700$ G and a pressure 6 mTorr at time 10 mks. Most of molecular ions have energies in the range from 0 to 0.4 eV. As a result, it is challenging to ascertain the temperature of the ions. This behavior is typical for all calculations. The electron distribution function was employed, with energies ranging from 0 to 20 eV, and approximated with a Maxwell function to estimate the electron temperature in the simulations.

Figure 9 shows the dependence of electron temperature on time for $B=700$ G and pressure 6 mTorr. All calculations are characterized by a slightly higher electron temperature (⁠ $T e=10−14 eV$⁠) than observed in the experiment^58,59 (⁠ $T e=5−12$ eV). This discrepancy can be attributed to the fact that the modeling did not include inelastic scattering and excitation reactions for electrons and neutral particles. Accordingly, the electron temperature estimated in the calculations may exceed that observed in the experiment.

Figure 10 shows the discharge density and potential distributions in plasma at different time points when the discharge combustion reaches the stationary mode. Simulation data are presented for $B=700$ G and pressure 6 mTorr. In this mode, the main part of electrons forms a stable cloud rotating around the PIG source axis with a characteristic frequency of about $∼$3.5 MHz. However, the exact value for the rotation period may vary depending on the geometry and discharge parameters. The ions are distributed fairly evenly over a cylinder with a radius of $∼1−2 mm$⁠. The concentration of charged particles fluctuates slightly around the point of equilibrium and does not change over time (Fig. 10). Fluctuations in the discharge current are explained by the use of Monte–Carlo methods to calculate kinetic processes as well as a finite-difference scheme to solve the Poisson equation. The average electron concentration in the stationary mode is in the range from $6× 10 15$ to $8× 10 15 m − 3$ (maximum is $1.2× 10 16 m − 3$⁠) and it is slightly higher than the average molecular ion concentration $4−6× 10 15 m − 3$ (maximum is $8× 10 15 m − 3$⁠). This behavior is consistently observed across all calculations.

Figure 11 shows the discharge and extraction currents calculated for different macroparticle weights (⁠ $N P I M$ = 2500 and 5000). The discharge characteristics obtained in the modeling with a different number of macroparticles are in good agreement with the experimental data.

The primary objective of the PIC-MCC model was to achieve a stationary solution. Therefore, the experimental results in Fig. 11 were shifted along the time scale. The discharge ignition time in the experiment depends on factors such as pressure, voltage, and magnetic field, and is located within the time range of 2–20 mks from the start of the anode voltage supply. In Refs. 14 and 16, a more detailed analysis of the experimentally measured time for the flash delay (⁠ $t d$⁠) (which is defined as the time between voltage supply to the anode and the start of current registration), and pulse rise time (⁠ $t f r$⁠) is presented.

A series of Penning source calculations were carried out over 8 mks from a stationary state to study convergence over the grid size and the macroparticle number. The pressure was 4 mTorr, the macroparticle weight (⁠ $N P I M$⁠) varied from 2500 to 5000, and the grid size ranged from 0.16 to 0.25 mm. The results in Fig. 12 clearly show the independence of the number of macroparticles. The grid size does not affect the extraction current (⁠ $I e x$⁠), but it does have a slight effect on the discharge current (⁠ $I d$⁠). This is due to a better approximation of electric field gradients on a smaller grid during the plasma discharge.

Figure 13 shows the dependence of discharge and extract currents on pressure with both calculated and experimental data for various magnetic fields. The study was conducted in an isotropic magnetic field with a strength of 700 G (the PIG source was inside a solenoid), and in an anisotropic magnetic field (permanent magnets were used). The configuration of the anisotropic magnetic field was determined experimentally and is detailed in Ref. 19. Note that in an experiment with a constant magnetic field at a pressure below $1−2 mTorr$⁠, and with a pulsed power supply mode, the discharge does not ignite (Fig. 13). At pressures between 2 and 6 mTorr, stable discharge combustion is observed. At pressures above 8 mTorr, a sharp increase in the discharge current occurs regardless of the magnetic field configuration. A similar pattern is observed in modeling, and the obtained values for discharge and extraction current coincide in calculations with different weights of macroparticles.

Once quasi-stationary solutions are obtained for modeling the ion source, the next step is to proceed with numerical calculations for the particle motion in the NT IOS and their interaction with the target. As already mentioned, this approach can significantly reduce the modeling time. The complete neutron tube system was observed for a $1−2 mks$ until the ion current on the target reached a stable state. The total ion current on the target exhibited a range of 60%–90% of the extracted ion current from Penning, contingent upon the pressure and accelerating voltage. Figure 14 shows the numerical results of ions and electrons motion in the PIG source and IOS. Calculations at 60 and 120 kV of accelerating voltage are given as examples. To suppress the electron current due to the secondary ion-electron emission, a magnetic field is created in the target node. A uniform magnetic field with a strength of $B=1000$ G perpendicular to the NT axis was created in the target area (highlighted in yellow in Fig. 14) to account for this effect. This figure also shows the positions of $D 2 +$ in red and $e −$ in blue at some point. As can be seen, some of the ions fall on the accelerating electrode, which is a parasitic effect, since these ions do not contribute to the generation of neutrons. Figure 14 clearly shows that not all emission electrons are suppressed in the target node. The region where the magnetic field turns to zero, the concentration of electrons decreases as a result of their drift to the ion source. Some of the emitted electrons, as well as $e −$⁠, formed as a result of the residual gas ionization by deuterium ions, move toward the PIG source and focus into a narrow beam.

The dependences of the ion current distribution density on the target radius, as well as the total ion current, were obtained as a result of this modeling. Figure 15 illustrates the current distribution of $D 2 +/ D +$ on the target surface in both 1D and 2D cases. The $D 2 +/ D +$ distribution is axisymmetric, with a maximum concentration near the target center. At this point, the ion current distribution density can be divided into several areas: the central part with a radius of about 0.1 mm for D $2 +$ and 1 mm for D $+$ (here the current density practically does not depend on the radius, probably due to rounding errors) with the highest current density; next, the middle part is a circle with a radius of about 2 mm, where the ion current density decreases almost linearly from the radius; and finally the peripheral part is a circle with a radius of 4 mm, where there is a lower rate of drop in current density from the radius. Almost the entire ion current is within a circle with a radius of $5−6 mm$⁠. The ion distribution densities on the target obtained in the calculation correlate well with the previously measured ion current profile.^98,99

The target heating and the thermal desorption of H isotopes under the influence of an ion beam were calculated, considering the dependence of the ion current distribution density on the target radius, as well as the total ion current. In addition, the dependences of the neutron yield on time were obtained. The modeling considered the target material sputtering and changing the characteristics of the target node exposed to irradiation. The SRIM software was used to obtain the degradation coefficients. Target sputtering coefficients were determined based on the incident particle energy and their isotopic composition. Calculations were performed for hydrogen and its isotopes (D and T), as the sputtering coefficient is proportional to the incident ion mass. The OpenFOAM software was used to calculate the diffusion and thermal desorption of H isotopes over time. The numerical solution was obtained by solving systems of differential equations describing this process. The model of diffusion and thermal desorption of H isotopes was used, considering the hydrogen sublattice.^89–91  Figure 16 displays the target degradation/sputtering over time, as well as the distribution of H isotope concentrations (D and T) across the target as a result of diffusion/thermal desorption of H isotopes. Figure 16 shows the substrate in orange, and the titanium hydride layer in turquoise-purple. The concentration of H isotopes ranges from 0 to 2.

Figure 16 clearly shows that the distribution of H isotopes throughout the thickness is almost uniform. The central target part with a radius of about 2 mm becomes saturated with H isotopes after 10 h of NT operation. After 150 h of operation, the radius of the saturated area increases up to 3 mm. It is worth noting the very high rate of target degradation due to ion irradiation. After 100 h of NT operation, the target central part is sputtered onto the substrate.

Figure 16 also depicts the distribution of the neutron generation probability over thickness and radius, as obtained in the Geant4 software with the CHIPS+TPT library. Those already mentioned distributions of H isotopes by thickness and radius, considering the radiation sputtering, were used to calculate the neutron yield using Monte–Carlo methods. It should be noted that the thickness of the region, where nuclear reactions occur with the release of neutrons, is about $300−400 nm$⁠. The central part stops making a significant contribution to the neutron flux at a given focusing of the ion beam on the target in $100−150$ h of operation.

The purpose of these experiments was to verify the neutron flux calculation method by studying the neutron flux dependence at different accelerating voltages (70, 80, and 90 kV). To ensure accurate results, we used NT with a short operating time (⁠ $30−60$ h), since during this time the target degradation does not yet affect the neutron flux. When modeling, a “perfect” target was used (excluding sputtering). As shown previously, the actual target is close to the “perfect” one when the NT operation time is less than 100 h. Figure 17 shows the modeling results and experimental dependences of the neutron flux on the tube current at various accelerating voltages and they are in good agreement.

Figure 18 compares the experimental and calculated neutron yield with time using a serial sample. It is clearly seen that the calculated neutron flux almost does not change with time when the operating time is less than 100 h, and this is due to the fact that the main contribution to neutrons is made by the central part of the target node, which has not yet sputtered to the substrate. Additionally, the neutron flux reaches its maximum in 50 h. This is due to the fact that in shorter times, the radial diffusion of H isotopes occurs, i.e., the D/T emission from the central part. And finally in typical $300−350$ h, the neutron flux drops down to $10 8n/s$⁠. Experimental data show a similar neutron yield behavior.

This paper describes a full-scale numerical NT model with the PIG source, which is tested in experimental measurements. The amplitude-time parameters of current waveforms are compared. The charge density behavior over time in the PIG source in a quasi-stationary mode is shown. Based on previous works the following results were calculated, experimentally measured and validated: ion current densities on the target, energy and mass-charge ions spectra extracted from ions of the PIG source. The neutron yield dependences on the discharge parameters, extraction current and accelerating voltage, as well as the neutron flux dependence on the operating time were experimentally studied.

The present study shows that the developed full-scale NT numerical model describes the observed experimental data. This allows the calculation and theoretical results to be used for further design and development of various types of NT with the PIG source. Also, existing tubes can be optimized to improve technical characteristics.

The authors express their gratitude to S. E. Kuratov, M. G. Lobok, D. A. Storozhev, and A. S. Dikalyuk for their valuable consultations and useful recommendations at all work stages.

The authors have no conflicts to disclose.

A. Rokhmanenkov: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). N. Mamedov: Methodology (equal). I. Kanshin: Methodology (equal). S. Maslennikov: Methodology (equal). A. Solodovnikov: Methodology (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.