Scaling Laws for Plasma Focus Machines from Numerical Experiments

Numerical experiments carried out systematically using the Lee Model code unveil insightful and practical wide-ranging scaling laws for plasma focus machines for nuclear fusion energy as well as other applications. An essential feature of the numerical experiments is the fitting of a measured current waveform to the computed waveform to calibrate the model for the particular machine, thus providing a reliable and rigorous determination of the all-important pinch current. The thermodynamics and radiation properties of the resulting plasma are then reliably determined. This paper provides an overview of the recently published scaling laws for neutron (Yn) and neon soft x-ray, SXR (Ysxr) yields: Yn = 3.2x10 Ipinch; Yn = 1.8x10 Ipeak; Ipeak (0.3 to 5.7), Ipinch (0.2 to 2.4) in MA. Yn~E0 at tens of kJ to Yn~E0 at MJ level (up to 25MJ) and Ysxr = 8.3x10 Ipinch; Ysxr = 6x10 Ipeak; Ipeak (0.1 to 2.4), Ipinch (0.07 to1.3) in MA. Ysxr~E0 (kJ range) to Ysxr~E0 (towards MJ).


Introduction
Plasma focus machines of various energies are increasingly being studied as sources of neutrons and soft x-rays.The most exciting prospect is for scaling the plasma focus up to regimes relevant for fusion energy studies.However, even a simple machine such as the UNU/ICTP PFF 3 kJ machine consistently produces 10 8 neutrons when operated in deuterium [1].Plasma focus machines operated in neon have also been studied as intense sources of soft x-rays with potential applications [2][3][4].Whilst many recent experiments have concentrated efforts on low energy devices [2][3][4] with a view of operating these as repetitive pulsed sources, other experiments have looked at x-ray pulses from larger plasma focus devices [5,6] extending to the MJ regime.Numerical experiments simulating x-ray pulses from plasma focus devices are also gaining more interest in the public domain.For example, the Institute of Plasma Focus Studies [7] conducted a recent international Internet Workshop on Plasma Focus Numerical Experiments [8], at which it was demonstrated that the Lee model code [9] not only computes realistic focus pinch parameters, but also absolute values of neutron yield Y n and soft x-ray yield Y sxr which are consistent with those measured experimentally.A com-parison was made for the case of the NX2 machine [4], showing good agreement between computed and measured Y sxr as a function of P 0 [8,10].This gives confidence that the Lee model code gives realistic results in the computation of Y n and Y sxr .
In this paper, we show the comprehensive range of numerical experiments conducted to derive scaling laws on neutron yield Y n [11,12] and neon Y sxr , in terms of storage energy E 0 , peak discharge current I peak and peak focus pinch current I pinch obtained from studies carried out over E 0 varying from 0.2 kJ to 25 MJ for optimised machine parameters and operating parameters.It is worth mentioning that the scaling laws in terms of I pinch and I peak have also been obtained for numerical experiments using the Lee model code fitted with the actual machine parameters and operating parameters and the difference from that obtained for the optimised conditions are within the order of 0.1 in the scaling laws power factor for neutrons and no change for neon SXR yield with I pinch .
We also wish to point out that the distinction of I pinch from I peak is of basic importance [13][14][15].The scaling with I pinch is the more fundamental and robust one; since obviously there are situations (no pinching or poor pinching however optimized) where I peak may be large but Y n is zero or small; whereas the scaling with I pinch is certainly more consistent with all situations.In these works the primary importance of I pinch for scaling plasma focus properties including neutron yield Y n , has been firmly established [11][12][13][14][15].

The Lee Model Code
The Lee model code couples the electrical circuit with plasma focus dynamics, thermodynamics and radiation, enabling realistic simulation of all gross focus properties.The basic model, described in 1984 [16] was successfully used to assist several projects [17][18][19].Radiationcoupled dynamics was included in the five-phase code leading to numerical experiments on radiation cooling [20].The vital role of a finite small disturbance speed discussed by Potter in a Z-pinch situation [21] was incorporated together with real gas thermodynamics and radiation-yield terms.Before this 'communication delay effect' was incorporated, the model consistently overestimated the radial speeds.This is serious from the point of view of neutron yields.A factor of two in shock speeds gives a factor of four in temperatures leading to a difference in fusion cross-sections of~1000 at the range of temperatures we are dealing with.This version of the code assisted other research projects [22][23][24][25][26][27] and was web-published in 2000 [28] and 2005 [29].Plasma self-absorption was included in 2007 [27] improving SXR yield simulation.The code has been used extensively in several machines including UNU/ICTP PFF [1,17,22,23,[25][26][27]30,31], NX2 [24,27,32], NX1 [3,32] and adapted for the Filippov-type plasma focus DENA [33].A recent development is the inclusion of the neutron yield Y n using a beam-target mechanism [11,12,14,15,34], incorporated in recent versions [9] of the code (versions later than RADPFV5.13),resulting in realistic Y n scaling with I pinch [11,12].The versatility and utility of the model are demonstrated in its clear distinction of I pinch from I peak [13] and the recent uncovering of a plasma focus pinch current limitation effect [14,15].The description, theory, code and a broad range of results of this 'Universal Plasma Focus Laboratory Facility' are available for download from [9].
A brief description of the code is given below.The five phases are summarised as follows: 1) Axial Phase: Described by a snowplow model with an equation of motion coupled to a circuit equation.The equation of motion incorporates the axial phase model parameters: mass and current factors f m and f c respectively.The mass swept-up factor f m accounts for not only the porosity of the current sheet but also for the inclination of the moving current sheet-shock front structure and all other unspecified effects which have effects equivalent to increasing or reducing the amount of mass in the moving structure during the axial phase.The current factor f c accounts for the fraction of current effectively flowing in the moving structure (due to all effects such as current shedding at or near the back-wall and current sheet inclination).This defines the fraction of current effectively driving the structure during the axial phase.
2) Radial Inward Shock Phase: Described by four coupled equations using an elongating slug model.The first equation computes the radial inward shock speed from the driving magnetic pressure.The second equation computes the axial elongation speed of the column.The third equation computes the speed of the current sheath, also called the magnetic piston, allowing the current sheath to separate from the shock front by applying an adiabatic approximation.The fourth is the circuit equation.Thermodynamic effects due to ionization and excitation are incorporated into these equations, these effects being important for gases other than hydrogen and deuterium.Temperature and number densities are computed during this phase.A communication delay between shock front and current sheath due to the finite small disturbance speed is crucially implemented in this phase.The model parameters, radial phase mass swept-up and current factors f mr and f cr respectively are incorporated in all three radial phases.The mass swept-up factor f mr accounts for all mechanisms which have effects equivalent to increasing or reducing the amount of mass in the moving slug during the radial phase.The current factor f cr accounts for the fraction of current effectively flowing in the moving piston forming the back of the slug (due to all effects).This defines the fraction of current effectively driving the radial slug.
3) Radial Reflected Shock (RS) Phase: When the shock front hits the axis, because the focus plasma is collisional, a reflected shock develops which moves radially outwards, whilst the radial current sheath piston continues to move inwards.Four coupled equations are also used to describe this phase, these being for the reflected shock moving radially outwards, the piston moving radially inwards, the elongation of the annular column and the circuit.The same model parameters f mr and f cr are used as in the previous radial phase.The plasma temperature behind the RS undergoes a jump by a factor of approximately two.
4) Slow Compression (Quiescent) or Pinch Phase: When the out-going reflected shock hits the in-coming piston the compression enters a radiative phase, in which for gases such as neon radiation emission may actually enhance the compression, where we have included energy loss/gain terms from Joule heating and radiation losses into the piston equation of motion.Three coupled equations describe this phase; these being the piston radial motion equation, the pinch column elongation equation and the circuit equation, incorporating the same model parameters as in the previous two phases.Thermodynamic effects are incorporated into this phase.The duration of this slow compression phase is set as the time of transit of small disturbances across the pinched plasma column.The computation of this phase is terminated at the end of this duration.
4) Expanded Column Phase: To simulate the current trace beyond this point, we allow the column to suddenly attain the radius of the anode, and use the expanded column inductance for further integration.In this final phase the snowplow model is used, and two coupled equations are used; similar to the axial phase above.This phase is not considered important as it occurs after the focus pinch.

Computation of Neutron Yield
The neutron yield is computed using a phenomenological beam-target neutron generating mechanism described recently by Gribkov et al [34] and adapted to yield the following equation.A beam of fast deuteron ions is produced by diode action in a thin layer close to the anode, with plasma disruptions generating the necessary high voltages.The beam interacts with the hot dense plasma of the focus pinch column to produce the fusion neutrons.The beam-target yield is derived [11,12,14,28] as: where n i is the ion density, b is the cathode radius, r p is the radius of the plasma pinch with length z p , σ the cross-section of the D-D fusion reaction, n-branch [35] and U, the beam energy.C n is treated as a calibration constant combining various constants in the derivation process.
The D-D cross-section is sensitive to the beam energy in the range 15-150kV; so it is necessary to use the appropriate range of beam energy to compute σ.The code computes induced voltages (due to current motion inductive effects) V max of the order of only 15-50 kV.However it is known, from experiments that the ion energy responsible for the beam-target neutrons is in the range 50-150 keV [34], and for smaller lowervoltage machines the relevant energy could be lower at 30-60 keV [31].Thus in line with experimental observations the D-D cross section σ is reasonably obtained by using U=3V max .This fit was tested by using U equal to various multiples of V max .A reasonably good fit of the computed neutron yields to the measured published neutron yields at energy levels from sub-kJ to near MJ was obtained when the multiple of 3 was used; with poor agreement for most of the data points when for example a multiple of 1 or 2 or 4 or 5 was used.The model uses a value of C n =2.7x10 7 obtained by calibrating the yield [9,13,14] at an experimental point of 0.5 MA.
The thermonuclear component is also computed in every case and it is found that this component is negligible when compared with the beam-target component.

Computation of Neon SXR Yield
We note that the transition from Phase 4 to Phase 5 is observed in laboratory measurements to occur in an extremely short time with plasma/current disruptions resulting in localized regions of high densities and temperatures.These localized regions are not modelled in the code, which consequently computes only an average uniform density, and an average uniform temperature which are considerably lower than measured peak density and temperature.However, because the 4 model parameters are obtained by fitting the computed total current waveform to the measured total current waveform, the model incorporates the energy and mass balances equivalent, at least in the gross sense, to all the processes which are not even specifically modelled.Hence the computed gross features such as speeds and trajectories and integrated soft x-ray yields have been extensively tested in numerical experiments for several machines and are found to be comparable with measured values.
In the code [9], neon line radiation Q L is calculated as follows: where for the temperatures of interest in our experiments we take the SXR yield Y sxr = Q L .Z n is the atomic number.
Hence the SXR energy generated within the plasma pinch depends on the properties: number density n i , effective charge number Z, pinch radius r p , pinch length z f and temperature T. It also depends on the pinch duration since in our code the Q L is obtained by integrating over the pinch duration.This generated energy is then reduced by the plasma self-absorption which depends primarily on density and temperature; the reduced quantity of energy is then emitted as the SXR yield.These effects are included in the modelling by computing volumetric plasma self-absorption factor A derived from the photonic excitation number M which is a function of Z n , n i , Z and T.However, in our range of operation, the numerical experiments show that the self absorption is not significant.It was first pointed out by Liu Mahe [23] that a temperature around 300 eV is optimum for SXR production.Shan Bing's subsequent work [24] and our experience through numerical experiments suggest that around 2x10 6 K (below 200 eV) or even a little lower could be better.Hence unlike the case of neutron scaling, for SXR scaling there is an optimum small range of temperatures (T windows) to operate.

Numerical Experiments
The Lee code is configured to work as any plasma focus by inputting the bank parameters, L 0 , C 0 and stray circuit resistance r 0 ; the tube parameters b, a and z 0 and operational parameters V 0 and P 0 and the fill gas.The standard practice is to fit the computed total current waveform to an experimentally measured total current waveform [11,[13][14][15]28,29] using the four model parameters representing the mass swept-up factor f m , the plasma current factor f c for the axial phase and factors f mr and f cr for the radial phases.
From experience it is known that the current trace of the focus is one of the best indicators of gross performance.The axial and radial phase dynamics and the crucial energy transfer into the focus pinch are among the important information that is quickly apparent from the current trace.
The exact time profile of the total current trace is governed by the bank parameters, by the focus tube geometry and the operational parameters.It also depends on the fraction of mass swept-up and the fraction of sheath current and the variation of these fractions through the axial and radial phases.These parameters determine the axial and radial dynamics, specifically the axial and radial speeds which in turn affect the profile and magnitudes of the discharge current.The detailed profile of the discharge current during the pinch phase also reflects the Joule heating and radiative yields.At the end of the pinch phase the total current profile also reflects the sudden transition of the current flow from a constricted pinch to a large column flow.Thus the discharge current powers all dynamic, electrodynamic, thermodynamic and radiation processes in the various phases of the plasma focus.Conversely all the dynamic, electrodynamic, thermodynamic and radiation processes in the various phases of the plasma focus affect the discharge current.It is then no exaggeration to say that the discharge current waveform contains information on all the dynamic, electrodynamic, thermodynamic and radiation processes that occur in the various phases of the plasma focus.This explains the importance attached to matching the computed current trace to the measured current trace in the procedure adopted by the Lee model code.

Scaling Laws for Neutrons from Numerical Experiments over a Range of Energies from 10kJ to 25 MJ
We apply the Lee model code to the MJ machine PF1000 over a range of C 0 to study the neutrons emitted by PF1000-like bank energies from 10kJ to 25 MJ.First, we fitted a measured current trace to obtain the model parameters.A measured current trace of the PF1000 with C 0 =1332 μF, operated at 27 kV, 3.5 torr deuterium, has been published [34], with cathode/anode radii b=16 cm, a=11.55 cm and anode length z 0 =60cm.In the numerical experiments we fitted external (or static) inductance L 0 = 33.5 nH and stray resistance r 0 =6.1 mΩ (damping factor RESF=r 0 /(L 0 /C 0 ) 0.5 =1.22).The fitted model parameters are: f m =0.13, f c =0.7, f mr =0.35 and f cr = 0.65.The computed current trace [11], [15] agrees very well with the measured trace through all the phases; axial and radial, right down to the bottom of the current dip indicating the end of the pinch phase as shown in Figure 1.
This agreement confirms the model parameters for the PF1000.Once the model parameters have been fitted to a machine for a given gas, these model parameters may be used with some degree of confidence when operating parameters such as the voltage are varied [9].With no measured current waveforms available for the higher megajoule numerical experiments, it is reasonable to keep the model parameters that we have got from the PF1000 fitting.
This series of numerical experiments is carried out at 35 kV, 10 torr deuterium, inductance L 0 = 33.5 nH, stray resistance r 0 =6.1 mΩ (damping factor RESF= r 0 / (L 0 /C 0 ) 0.5 =1.22).The ratio c=b/a is retained at 1.39.The numerical experiments were carried out for C 0 ranging from 14 µF to 39960 µF corresponding to energies from 8.5 kJ to 24 MJ [12].For each C 0 , anode length z 0 is varied to find the optimum.For each z 0 , anode radius a 0 is varied so that the end axial speed is 10 cm/µs.
For this series of experiments we find that the Y n scaling changes from Y n ~E0 2.0 at tens of kJ to Y n ~E0 0.84 at the highest energies (up to 25MJ) investigated in this series.This is shown in Figure 2. The scaling of Y n with I peak and I pinch over the whole range of energies investigated up to 25 MJ (Figure 3) is as follows: Y n = 3.2x10 11 I pinch 4.5 and where I peak ranges from 0.3 to 5.7 MA and I pinch ranges from 0.2 to 2.4 MA.This compares to an earlier study carried out on several machines with published current traces with Y n yield measurements, operating conditions and machine parameters including the PF400, UNU/ICTP PFF, the NX2 and Poseidon providing a slightly higher scaling laws: Y n ~Ipinch 4.7 and Y n ~Ipeak 3.9 . The slightly higher value of the scaling is because those machines fitted are of mixed 'c' mixed bank parameters, mixed model parameters and currents generally below 1MA and voltages generally below the 35 kV [11].

Scaling Laws for Neon SXR from Numerical
Experiments over a Range of Energies from 0.2 kJ to 1 MJ We next use the Lee model code to carry out a series of numerical experiments to obtain the soft x-ray yield in neon for bank energies from 0.2 kJ to 1 MJ [36].In this case we apply it to a proposed modern fast plasma focus machine with optimised values for c the ratio of the outer to inner electrode radius and L 0 obtained from our numerical experiments.
The following parameters are kept constant: 1) the ratio c=b/a (kept at 1.5, which is practically optimum according to our preliminary numerical trials; 2) the operating voltage V 0 (kept at 20 kV); 3) static inductance L 0 (kept at 30 nH, which is already low enough to reach the I pinch limitation regime [13,14] over most of the range of E 0 we are covering) and; 4) the ratio of stray resistance to surge impedance RESF (kept at 0.1, representing a higher performance modern capacitor bank).The model parameters [8][9][10][11][12][13][14] f m , f c , f mr , f cr are also kept at fixed values 0.06, 0.7, 0.16 and 0.7.We choose the model parameters as they represent the average values from the range of machines that we have studied.A typical current waveform is shown in Figure 4.The storage energy E 0 is varied by changing the capacitance C 0 .Parameters that are varied are operating pressure P 0 , anode length z 0 and anode radius 'a'.Parametric variation at each E 0 follows the order; P 0 , z 0 and a until all realistic combinations of P 0 , z 0 and a are investigated.At each E 0 , the optimum combination of P 0 , z 0 and a is found that produces the biggest Y sxr .In other words at each E 0 , a P 0 is fixed, a z 0 is chosen and a is varied until the largest Y sxr is found.Then keeping the same values of E 0 and P 0 , another z 0 is chosen and a is varied until the largest Y sxr is found.This procedure is repeated until for that E 0 and P 0 , the optimum combination of z 0 and a is found.Then keeping the same value of E 0 , another P 0 is selected.The procedure for parametric variation of z 0 and a as described above is then carried out for this E 0 and new P 0 until the optimum combination of z 0 and a is found.This procedure is repeated until for a fixed value of E 0 , the optimum combination of P 0 , z 0 and a is found.The procedure is then repeated with a new value of E 0 .In this manner after systematically carrying out some 2000 runs, the optimized runs for various energies are obtained.A plot Y sxr against E 0 is shown in Figure 5.
We then plot Y sxr against I peak and I pinch and obtain SXR yield scales as Y sxr ~Ipinch 3.6 and Y sxr ~Ipeak 3.2 .The I pinch scaling has less scatter than the I peak scaling.We next subject the scaling to further test when the fixed parameters RESF, c, L 0 and V 0 and model parameters f m , f c , f mr , f cr are varied.We add in the results of some numerical experiments using the parameters of several existing plasma focus devices including the UNU/ICTP PFF (RESF =0.We would like to highlight that the consistent behaviour of I pinch in maintaining the scaling of Y sxr ~ I pinch 3.6 with less scatter than the Y sxr ~Ipeak 3.2 scaling particularly when mixed-parameters cases are included, strongly support the conclusion that I pinch scaling is the more universal and robust one.Similarly conclusions on the importance of I pinch in plasma focus performance and scaling laws have been reported [11][12][13][14][15].It may also be worthy of note that our comprehensively surveyed numerical experiments for Mather configurations in the range of energies 0.2 kJ to 1 MJ produce an I pinch scaling rule for Y sxr not compatible with Gates' rule [37].However it is remarkable that our I pinch scaling index of 3.6, obtained through a set of comprehensive numerical experiments over a range of energies 0.2 kJ to 1 MJ, on Mather-type devices is within the range of 3.5to4 postulated on the basis of sparse experimental data, (basically just two machines one at 5 kJ and the other at 0.9 MJ), by Filippov [6], for Filippov configurations in the range of energies 5 kJ to 1 MJ.
It must be pointed out that the results represent scaling for comparison with baseline plasma focus devices that have been optimized in terms of electrode dimensions.It must also be emphasized that the scaling with I pinch works well even when there are some variations in the actual device from L 0 = 30 nH, V 0 = 20 kV and c = 1.5.However there may be many other parameters which can change and could lead to a further enhancement of x-ray yield.

Conclusions
Numerical experiments carried out using the universal plasma focus laboratory facility based on the Lee model code gives reliable scaling laws for neutrons production and neon SXR yields for plasma focus machines.The scaling laws obtained: ; I peak (0.1 to 2.4), I pinch (0.07 to1.3) in MA.Y sxr ~E0 1.6 (kJ range) to Y sxr ~E0 0.8 ( towards MJ).
These laws provide useful references and facilitate the understanding of present plasma focus machines.More importantly, these scaling laws are also useful for design considerations of new plasma focus machines particularly if they are intended to operate as optimized neutron or neon SXR sources.More recently, the scaling of Y n versus E 0 as shown above has been placed in the context of a global scaling law [38] with the inclusion of available experimental data.From that analysis, the cause of scaling deterioration for neutron yield versus energy as shown in Figure 2 (which has also been given the misnomer 'neutron saturation') has been uncovered as due to a current scaling deterioration caused by an almost constant axial phase 'dynamic resistance' interacting with a reducing bank impedance as energy storage is increased at essentially constant voltage.Solutions suggested include the use of ultra-high voltages and circuit enhancement techniques such as current-steps [39,40].It is suggested here that the deterioration of soft x-ray yield with storage energy as shown in Figure 5 could also be ascribed to the same axial phase 'dynamic resistance' effect as described in that reference [38].

Figure 1 .Figure 2 .
Figure 1.Current fitting of computed current to measured current traces to obtain fitted parameters f m = 0.13, f c = 0.7, f mr = 0.35 and f cr = 0.65

Figure 3 .Figure 4 .
Figure 3. Log(Y n ) scaling with Log(I peak ) and Log(I pinch ), for the range of energies investigated, up to 25 MJ

Figure 6 .
Figure 6.Y sxr is plotted as a function of I pinch and I peak .The parameters kept constant for the black data points are: RESF = 0.1, c = 1.5, L 0 = 30nH and V 0 = 20 kV and model parameters m , f c , f mr , f cr at 0.06, 0.7, 0.16 and 0.7 respectively.The unblackened data points are for specific machines which have different values for the parameters c, L 0 , V 0 and RESF