Numerical Experiments Providing New Insights into Plasma Focus Fusion Devices

Recent extensive and systematic numerical experiments have uncovered new insights into plasma focus fusion devices including the following: (1) a plasma current limitation effect, as device static inductance is reduced towards very small values; (2) scaling laws of neutron yield and soft x-ray yield as functions of storage energies and currents; (3) a global scaling law for neutron yield as a function of storage energy combining experimental and numerical data showing that scaling deterioration has probably been interpreted as neutron 'saturation'; and (4) a fundamental cause of neutron 'saturation'. The groundbreaking insights thus gained may completely change the directions of plasma focus fusion research.


Introduction
The Lee model code couples the electrical circuit with plasma focus dynamics, thermodynamics, and radiation, enabling a realistic simulation of all gross focus properties.The basic model, described in 1984 [1], was successfully used to assist several projects [2][3][4].Radiation-coupled dynamics was included in the five-phase code, leading to numerical experiments on radiation cooling [5].The vital role of a finite small disturbance speed discussed by Potter in a Z-pinch situation [6] was incorporated together with real gas thermodynamics and radiation-yield terms.This version of the code assisted other research projects [7][8][9][10][11][12] and was web published in 2000 [13] and 2005 [14].Plasma OPEN ACCESS self-absorption was included in 2007 [13], improving the SXR yield simulation.The code has been used extensively in several machines including UNU/ICTP PFF [2,7,8,10,11,[15][16][17], NX2 [9,12,18], and NX1 [18,19] and has been adapted for the Filippov-type plasma focus DENA [20].A recent development is the inclusion of the neutron yield Y n using a beam-target mechanism [21][22][23][24][25], incorporated in recent versions [26,27] of the code (versions later than RADPFV5.13),resulting in realistic Y n scaling with I pinch [21,22].The versatility and utility of the model are demonstrated in its clear distinction of I pinch from I peak [28] and the recent uncovering of a plasma focus pinch current limitation effect [23,24,29] as static inductance is reduced towards zero.Extensive numerical experiments had been carried out systematically resulting in the uncovering of neutron [21][22][30][31][32][33] and SXR [30][31][32][33][34][35][36][37] scaling laws over a wider range of energies and currents than attempted before.The numerical experiments also gave insight into the nature and cause of 'neutron saturation [31,33,38].The description, theory, code, and a broad range of results of this "Universal Plasma Focus Laboratory Facility" are available for download from [26].
A brief description of the 5-phase model is given in the following.

The 5-Phase Lee Model Code
The five phases (a-e) are summarized [26,27,31,36,39] as follows: a. Axial Phase (see Figure 1 left part): Described by a snowplow model with an equation of motion which is coupled to a circuit equation.The equation of motion incorporates the axial phase model parameters: mass and current factors f m and f c .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, boundary layer effects, 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.
b. Radial Inward Shock Phase (see Figure 1 right part, also Figure 2): 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, (magnetic piston), allowing the current sheath to separate from the shock front by applying an adiabatic approximation [6].The fourth is the circuit equation.Thermodynamic effects due to ionization and excitation are incorporated into these equations, these effects being particularly important for gases other than hydrogen and deuterium.Temperature and number densities are computed during this phase using shock-jump equations.A communication delay between shock front and current sheath due to the finite small disturbance speed [6,26,39] is crucially implemented in this phase.The model parameters, radial phase mass swept-up and current factors f mr and f cr 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.
Figure 1.Schematic of the axial and radial phases.The left section depicts the axial phase, the right section the radial phase.In the left section, z is the effective position of the current sheath-shock front structure.In the right section r s is the position of the inward moving shock front driven by the piston at position r p .Between r s and r p is the radially imploding slug, elongating with a length z f .The capacitor, static inductance and switch powering the plasma focus is shown for the axial phase schematic only.c.Radial Reflected Shock (RS) Phase (See Figure 2): 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 reflected shock undergoes a jump by a factor close to 2. Number densities are also computed using the reflected shock jump equations.2): When the out-going reflected shock hits the inward moving 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.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.
e. 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.
We note [36] that in radial phases b, c and d, axial acceleration and ejection of mass caused by necking curvatures of the pinching current sheath result in time-dependent strongly center-peaked density distributions.Moreover the transition from phase d to phase e 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 center-peaking density effects and localized regions are not modeled 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 four 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 modeled.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.

From Measured Current Waveform to Modeling for Diagnostics
The Lee model code [26,27] is configured [21−24,26−39] to work as any plasma focus by inputting:  Bank parameters, L 0 , C 0 and stray circuit resistance r 0 ;  Tube parameters b, a and z 0 and  Operational parameters V 0 and P 0 and the fill gas.
The computed total current waveform is fitted to the measured waveform by varying model parameters f m , f c , f mr and f cr one by one, until the computed waveform agrees with the measured waveform.
First, the axial model factors f m , f c are adjusted (fitted) until the features (1) computed rising slope of the total current trace and (2) the rounding off of the peak current as well as (3) the peak current itself are in reasonable (typically very good) fit with the measured total current trace (see Figure 3, measured trace fitted with computed trace).
Then we proceed to adjust (fit) the radial phase model factors f mr and f cr until features (4) the computed slope and (5) the depth of the dip agree with the measured.Note that the fitting of the computed trace with the measured current trace is done up to the end of the radial phase which is typically at the bottom of the current dip.Fitting of the computed and measured current traces beyond this point is not done.If there is significant divergence of the computed with the measured trace beyond the end of the radial phase, this divergence is not considered important.
In this case, after fitting the 5 features (1) to (5) above, the following fitted model parameters are obtained: f m = 0.1, f c = 0.7, f m r = 0.12, f cr = 0.68.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 [26,27,31].
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.There are many underlying mechanisms in the axial phase such as shock front and current sheet structure, porosity and inclination, boundary layer effects and current shunting and fragmenting which are not simply modeled; likewise in the radial phase mechanisms such as current sheet curvatures and necking leading to axial acceleration and ejection of mass, and plasma/current disruptions.These effects may give rise to localized regions of high density and temperatures.The detailed profile of the discharge current is influenced by these effects and 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 total current trace to the measured total current trace in the procedure adopted by the Lee model code.Once matched, the fitted model parameters assure that the computation proceeds with all physical mechanisms accounted for, at least in the gross energy and mass balance sense.

Diagnostics-Time Histories of Dynamics, Energies and Plasma Properties Computed from the Measured Total Current Waveform by the Code
During every adjustment of each of the model parameters the code goes through the whole cycle of computation.In the last adjustment, when the computed total current trace is judged to be reasonably well fitted in all five waveform features, computed time histories are presented, in Figure 4a-o as an example, as follows: for the NX2 operated at 11 kV, 2.6 Torr neon [18,37].

Comments on Computed Quantities
The computed total current trace typically agrees very well with the measured because of the fitting.The end of the radial phase is indicated in Figure 4a.Plasma currents are rarely measured.We had done a comparison of the computed plasma current with measured plasma current for the Stuttgart PF78 which shows good agreement of our computed to the measured plasma current [28].The computed plasma current in this case of the NX2 is shown in Figure 4b.The computed tube voltage is difficult to compare with measured tube voltages in terms of peak values, typically because of poor response time of voltage dividers.However the computed waveform shape in Figure 4c. is generally as expected.The computed axial trajectory and speed, agree with experimentally obtained time histories.Moreover, the behavior with pressure, running the code again at different pressures, agrees well with experimental results.The radial trajectories and speeds are difficult to measure.The computed trajectories Figure 4e agree with the scant experimental data available.The length of the radial structure is shown in Figure 4f.Computed speeds radial shock front and piston speeds and speed of the elongation of the structure are shown in Figure 4g.The computed inductance (Figure 4h) shows a steady increase of inductance in the axial phase, followed by a sharp increase (rising by more than a factor of 2 in a radial phase time interval about 1/10 the duration of the axial phase for the NX2).

SXR emission power
The inductive energy (0.5 LI 2 ) peaks at 70% of initial stored energy, and then drops to 30% during the radial phase, as the sharp drop of current more than offsets the effect of sharply increased inductance (Figure 4i).In Figure 4j is shown the work done by the magnetic piston, computed using force integrated over distance method.Also shown is the work dissipated by the dynamic resistance, computed using dynamic resistance power integrated over time.We see that the two quantities and profiles agree exactly.This validates the concept of half Ldot as a dynamic resistance, DR (see Section 6.1).The piston work deposited in the plasma increases steadily to some 12% at the end of the axial phase and then rises sharply to just below 30% in the radial phase.Dynamic resistance is shown in Figure 4k.The values of the DR in the axial phase, together with the bank surge impedance, are the quantities that determine I peak .The ion number density has a maximum value derived from shock-jump considerations, and an averaged uniform value derived from overall energy and mass balance considerations.The time profiles of these are shown in the Figure 4l.The electron number density (Figure 4m) has similar profiles to the ion density profile, but is modified by the effective charge numbers due to ionization stages reached by the ions.Plasma temperature too has a maximum value and an averaged uniform value derived in the same manner; are shown in Figure 4n.Computed neon soft x-ray power profile is shown in Figure 4o.The area of the curve is the soft x-ray yield in Joule.Pinch dimensions and lifetime may be estimated from Figure 4e,f.The model also computes the neutron yield, for operation in deuterium, using a phenomenological beam-target mechanism [25−27].The model does not compute a time history of the neutron emission, only a yield number Y n .
Thus as is demonstrated above, the model code when properly fitted is able to realistically model any plasma focus and act as a guide to diagnostics of plasma dynamics, trajectories, energy distribution and gross plasma properties.
Moreover, using such simulation, series of experiments have been systematically carried out to look for behavior patterns of the plasma focus.Insights uncovered by the series of experiments include: (i) pinch current limitation effect as static inductance is reduced; (ii) neutron and SXR scaling laws; (iii) a global scaling law for neutrons versus storage energy combining experimental and numerical experimental data; and (iv) insight into the nature and a fundamental cause of neutron saturation.These are significant achievements accomplished within a period of twenty months of intensive numerical experimentation.

Insight 1-Pinch Current Limitation Effect as Static Inductance Is Reduced towards Zero
In a recent paper [25], there was expectation that the large MJ plasma focus PF1000 in Warsaw could increase its discharge current, and its pinch current, and consequently neutron yield by a reduction of its external or static inductance L 0 .To investigate this point, experiments were carried out using the Lee Model code.Unexpectedly, the results indicated that whilst I peak indeed progressively increased with reduction in L 0 , no improvement may be achieved due to a pinch current limitation effect [23−24].Given a fixed C 0 powering a plasma focus, there exists an optimum L 0 for maximum I pinch .Reducing L 0 further will increase neither I pinch nor Y n .The numerical experiments leading to this unexpected result is described below.
A measured current trace of the PF1000 with C 0 = 1332 μF, operated at 27 kV, 3.5 torr deuterium, has been published [25], with cathode/anode radii b = 16 cm, a = 11.55 cm and anode length z 0 = 60 cm.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 [21,24,26] 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 5.We carried out numerical experiments for PF1000 using the machine and model parameters determined from Figure 5. Operating the PF1000 at 35 kV and 3.5 Torr, we varied the anode radius a with corresponding adjustment to b to maintain a constant c = b/a = 1.39 and in order to keep the peak axial speed at 10 cm/s.The anode length z 0 was also adjusted to maximize I pinch as L 0 was decreased from 100 nH progressively to 5 nH.
As expected, I peak increased progressively from 1.66 to 4.4 MA.As L 0 was reduced from 100 to 35 nH, I pinch also increased, from 0.96 to 1.05 MA.However, then unexpectedly, on further reduction from 35 to 5 nH, I pinch stopped increasing, instead decreasing slightly to 1.03 MA at 20 nH, to 1.0 MA at 10 nH, and to 0.97 MA at 5 nH.Y n also had a maximum value of 3.2 × 10 11 at 35 nH.
To explain this unexpected result, we examine the energy distribution in the system at the end of the axial phase (see Figure 5) just before the current drops from peak value I peak and then again near the bottom of the almost linear drop to the pinch phase indicated by the arrow pointing to 'end of radial phase'.The energy equation describing this current drop is written as follows: where L a is the inductance of the tube at full axial length z 0 , δ plasma is the energy imparted to the plasma as the current sheet moves to the pinch position and is the integral of 0.5(dL/dt)I 2 .We approximate this as 0.5L p I pinch 2 which is an underestimate for this case.δ cap is the energy flow into or out of the capacitor during this period of current drop.If the duration of the radial phase is short compared to the capacitor time constant, the capacitor is effectively decoupled and δ cap may be put as zero.From this consideration we obtain where we have taken f c = 0.7 and approximated f c 2 as 0.5.
Generally, as L 0 is reduced, I peak increases; a is necessarily increased leading [7] to a longer pinch length z p , hence a bigger L p .Lowering L o also results in a shorter rise time, hence a necessary decrease in z 0 , reducing L a .Thus, from Equation (2), lowering L 0 decreases the fraction I pinch /I peak .Secondly, this situation is compounded by another mechanism.As L 0 is reduced, the L-C interaction time of the capacitor bank reduces while the duration of the current drop increases (see Figure 6, discussed in the next section) due to an increasing a.This means that as L 0 is reduced, the capacitor bank is more and more coupled to the inductive energy transfer processes with the accompanying induced large voltages that arise from the radial compression.Looking again at the derivation of Equation (2) from Equation (1) a nonzero δ cap , in this case, of positive value, will act to decrease I pinch further.The lower the L 0 the more pronounced is this effect.
Summarizing this discussion, the pinch current limitation is not a simple effect, but is a combination of the two complex effects described above, namely, the interplay of the various inductances involved in the plasma focus processes abetted by the increasing coupling of C o to the inductive energetic processes, as L 0 is reduced.

Optimum L 0 for Maximum Pinch Current and Neutron Yield
From the pinch current limitation effect, it is clear that given a fixed C 0 powering a plasma focus, there exists an optimum L 0 for maximum I pinch .Reducing L 0 further will increase neither I pinch nor Y n .The results of the numerical experiments carried out are presented in Figure 6 and Table 1.
With large L 0 = 100 nH it is seen (Figure 6) that the rising current profile is flattened from what its waveform would be if unloaded; and peaks at around 12 μs (before its unloaded rise time, not shown, of 18 μs) as the current sheet goes into the radial phase.The current drop, less than 25% of peak value, is sharp compared with the current rise profile.At L 0 = 30 nH the rising current profile is less flattened, reaching a flat top at around 5 μs, staying practically flat for some 2 μs before the radial phase current drop to 50% of its peak value in a time which is still short compared with the rise time.With L 0 of 5 nH, the rise time is now very short, there is hardly any flat top; as soon as the peak is reached, the current waveform droops significantly.There is a small kink on the current waveform of both the L 0 = 5 nH, z 0 = 20 cm and the L 0 = 5 nH, z 0 = 40 cm.This kink corresponds to the start of the radial phase which, because of the large anode radius, starts with a relatively low radial speed, causing a momentary reduction in dynamic loading.Looking at the three types of traces it is seen that for L 0 = 100 nH to 30 nH, there is a wide range of z 0 that may be chosen so that the radial phase may start at peak or near peak current, although the longer values of z 0 tend to give better energy transfers into the radial phase.
The optimized situation for each value of L 0 is shown in Table 1.The table shows that as L 0 is reduced, I peak rises with each reduction in L 0 with no sign of any limitation.However, I pinch reaches a broad maximum of 1.05 MA around 40-30 nH.Neutron yield Y n also shows a similar broad maximum peaking at 3.2 × 10 11 neutrons.Figure 7 shows a graphical representation of this I pinch limitation effect.The curve going up to 4 MA at low L 0 is the I peak curve.Thus I peak shows no sign of limitation as L 0 is progressively reduced.However I pinch reaches a broad maximum.From Figure 7 there is a stark and important message: one must distinguish clearly between I peak and I pinch .In general one cannot take I peak to be representative of I pinch .We carried out several sets of experiments on the PF1000 for varying L 0 , each set with a different damping factor.In every case, an optimum inductance was found around 30-60 nH with I pinch decreasing as L 0 was reduced below the optimum value.The results showed that for PF1000, reducing L 0 from its present 20-30 nH will increase neither the observed I pinch nor the neutron yield, because of the pinch limitation effect.Indeed, the I pinch decreases very slightly on further reduction to very small values.We would add that we have used a set of model parameters which in our experience is the most reasonable to be used in these numerical experiments.Variations of the model parameters could occur but we are confident that these variations are not likely to occur with such a pattern as to negate the pinch current limitation effect.Nevertheless these variations should be actively monitored and any patterns in the variations should be investigated.

Computation of Neutron Yield-describing the Beam-target Mechanism
The neutron yield is computed using a phenomenological beam-target neutron generating mechanism described recently by Gribkov et al. [25] 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 [21][22][23][24]26] 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 [40] 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−150 kV; 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 [25], and for smaller lower-voltage machines the relevant energy could be lower at 30-60 keV [16].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.7 × 10 7 obtained by calibrating the yield [21][22][23]26], 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.It might be argued that an adjustment to the thermonuclear component could also be attempted in a similar way to the usage of the multiple to V max .However, the usage of the multiple to V max has some experimental basis due to ion energy measurements.Moreover the value of Vmax in each numerical experiment is calculated from the slug model leading to the slow compression phase whilst it is known experimentally that after the slow compression phase, instability effects set in which will increase the electric fields operating within the pinch.These are the basic arguments supporting the view that the operational beam energy has a value above Vmax.For the thermonuclear component a feasible model to adjust the yield upwards has yet to be suggested.

Scaling Laws for Neutrons from Numerical Experiments over a Range of Energies from 10 kJ 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 10 kJ to 25 MJ.
As shown earlier the PF1000 current trace has been used to fit the model parameters, with very good fitting achieved between the computed and measured current traces (Figure 5).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 [26,27].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.
The optimum pressure for this series of numerical experiments is 10 torr and the ratio c = b/a is retained at 1.39.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.The numerical experiments were carried out for C 0 ranging from 14 µF to 39,960 µF corresponding to energies from 8.5 kJ to 24.5 MJ [22].
For this series of experiments we find that the Y n scaling changes from Y n ~ E 0 2.0 at tens of kJ to Y n ~ E 0 0.84 at the highest energies (up to 25 MJ) investigated in this series.This is shown in Figure 8.This compares to an earlier study carried out on several machines with published current traces and Y n yield measurements, operating conditions and machine parameters including the Chilean PF400J, the UNU/ICTP PFF, the NX2 and Poseidon providing a slightly higher scaling laws: Y n ~ I pinch 4.7 and Y n ~ I peak 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 1 MA and voltages generally below the 35 kV [21].Summarising: Over wide ranges of energy, optimizing pressure, anode length and radius, the scaling laws for Y n [15] are listed here: Y n = 3.2 × 10 11 I pinch 4.5 ; Y n = 1.8 × 10 10 I peak 3.8 ; I peak (0.3 to 5.7), I pinch (0.2 to 2.4) in MA.
Y n ~ E 0 2.0 at tens of kJ to Y n ~ E 0 0.84 at MJ level (up to 25MJ).
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 sources.

Computation of Neon SXR Yield-the Equations Used in the Computation
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 modeled 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 modeled.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 [8,9,26−27], 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 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 modeling 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 [8,11,18] that a temperature around 300 eV is optimum for SXR production.Shan Bing's subsequent work [9] and our experience through numerical experiments suggest that around 2 × 10 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.

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 over the energy range 0.2 kJ to 1 MJ [30,34].In this case we apply it to a proposed modern fast plasma focus machine with optimized 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 : (i) the ratio c = b/a (kept at 1.5, which is practically optimum according to our preliminary numerical trials; (ii) the operating voltage V 0 (kept at 20 kV); (iii) static inductance L 0 (kept at 30 nH, which is already low enough to reach the I pinch limitation regime [23,24] over most of the range of E 0 we are covering) and; (iv) the ratio of stray resistance to surge impedance RESF (kept at 0.1, representing a higher performance modern capacitor bank).The model parameters [26,27] 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 so they represent the average values from the range of machines that we have studied.A typical example of a current trace for these parameters is shown in Figure 10.
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 tabulated in Table 2.We then plot Y sxr against I peak and I pinch and obtain SXR yield scales as Y sxr ~ I pinch 3.6 and Y sxr ~ I peak 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 behavior of I pinch in maintaining the scaling of Y sxr ~ I pinch 3.6 with less scatter than the Y sxr ~ I peak 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 [21−24,26−28].It may also be worthy to 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 [41].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 0.2 kJ to 1 MJ, on Mather-type devices, is within the range of 3.5 to 4 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 [42], 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.
Summarising: over wide ranges of energy, optimizing pressure, anode length and radius, the scaling laws for neon SXR are found to be: where I peak (0.1 to 2.4), I pinch (0.07 to1.3) in MA, Y sxr ~ E 0 1.6 (kJ range) to Y sxr ~ E 0 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 neon SXR sources.

Insight 4-Neutron Saturation
Besides being accurately descriptive and related to wide-ranging experimental reality, desirable characteristics of a model include predictive and extrapolative scaling.Moreover a useful model should be accessible, usable and user-friendly and should be capable of providing insights.Insight however cannot be a characteristic of the model in isolation, but is the interactive result of the model with the modeler or model user.
It was observed early in plasma focus research [43] that neutron yield Y n ~ E 0 2 where E 0 is the capacitor storage energy.Such scaling gave hopes of possible development as a fusion energy source.Devices were scaled up to higher E 0 .It was then observed that the scaling deteriorated, with Y n not increasing as much as suggested by the E 0 2 scaling.In fact some experiments were interpreted as evidence of a neutron saturation effect [44,45] as E 0 approached several hundreds of kJ.As recently as 2006 Krauz [46] and November 2007, Scholz [47] have questioned whether the neutron saturation was due to a fundamental cause or to avoidable machine effects such as incorrect formation of plasma current sheath arising from impurities or sheath instabilities [45].We should note here that the region of discussion (several hundreds of kJ approaching the MJ region) is in contrast to the much higher energy region discussed by Schmidt at which there might be expected to be a decrease in the role of beam target fusion processes [45,48].

The Global Neutron Scaling Law
Recent extensive numerical experiments [21,22] also showed that whereas at energies up to tens of kJ the Y n ~ E 0 2 scaling held, deterioration of this scaling became apparent above the low hundreds of kJ.
This deteriorating trend worsened and tended towards Yn ~ E 0 0.8 at tens of MJ.The results of these numerical experiments are summarized in Figure 1 with the solid line representing results from numerical experiments.Experimental results from 0.4 kJ to MJ, compiled from several available published sources [43,45−47,49−52], are also included as squares in the same figure.The combined experimental and numerical experimental results [31][32][33]38] (Figure 12) appear to have general agreement particularly with regards to the Y n ~ E 0 2 at energies up to 100 kJ, and the deterioration of the scaling from low hundreds of kJ to the 1 MJ level.The global data of Figure 12 suggests that the apparently observed neutron saturation effect is overall not in significant variance with the deterioration of the scaling shown by the numerical experiments.

The Cause of Neutron 'Saturation' is the Dynamic Resistance
A simple yet compelling analysis of the cause of this neutron saturation has been published [38].In Figure 1 left side is shown a schematic of the plasma dynamics in the axial phase of the Mather-type plasma focus.In that work the simplest representation was used, in which the current sheet is shown to go from the anode to the cathode perpendicularly.Observation shows that there is actually a canting of the current sheet [53] and also that only a fraction (typically 0.7) of the total current participates in driving the current sheet.These points are accounted for in the modeling by model parameters f m and f c .We now represent the plasma focus circuit as shown in Figure 13.We consider only the axial phase.By surveying published results of all Mather-type experiments we find that all deuterium plasma focus devices operate at practically the same speeds [7] and are characterized by a constancy of energy density (per unit mass) over the whole range of devices from the smallest sub-kJ to the largest MJ devices.The time varying tube inductance is L = (/2)ln(c) z, where c = b/a and  is the permeability of free space.The rate of change of inductance is dL/dt = 2 × 10 −7 ln(c) dz/dt in SI units.Typically on switching, as the capacitor discharges, the current rises towards its peak value, the current sheet is accelerated, quickly reaching nearly its peak speed and continues accelerating slightly towards its peak speed at the end of the axial phase.Thus for most of its axial distance the current sheet is travelling at a speed close to the end-axial speed.In deuterium the end-axial speed is observed to be about 10 cm/s over the whole range of devices.This fixes the rate of change of inductance dL/dt as 1.4 × 10 −2 H/s for all the devices, if we take the radius ratio c = b/a = 2.This value of dL/dt changes by at most a factor of 2, taking into account the variation of c from low values of 1.4 (generally for larger machines) to 4 (generally for smaller machines).This typical dL/dt may also be expressed as 14 m We need now to inquire into the nature of the change in the inductance L(t).
Consider instantaneous power P delivered to L(t) by a change in L(t) Induced voltage: Hence instantaneous power into L(t): Next, consider instantanteous power associated with the inductive energy (½ LI 2 ): We note that P L of Equation ( 7) is not the same as P of Equation (6).The difference = P -P L = (½)(dL/dt)I 2 is not associated with the inductive energy stored in L. We conclude that whenever L(t) changes with time, the instantaneous power delivered to L(t) has a component that is not inductive.Hence this component of power (½)(dL/dt)I 2 must be resistive in nature; and the quantity (½)(dL/dt) also denoted as half Ldot is identified as a resistance, due to the motion associated with dL/dt ; which we call the dynamic resistance DR [27,31,33,38].Note that this is a general result and is independent of the actual processes involved.In the case of the plasma focus axial phase, the motion of the current sheet imparts power to the shock wave structure with consequential shock heating, Joule heating, ionization, radiation etc.The total power imparted at any instant is just the amount (½)(dL/dt)I 2 , with this amount powering all consequential processes.We denote the dynamic resistance of the axial phase as DR 0 .
We have thus identified for the axial phase of the plasma focus a typical dynamic resistance of 7 m due to the motion of the current sheet at 10 cm/s.It should be noted here that similar ideas of the role of dL/dt as a resistance were discussed by Bernard et al. [45].In that work the effect of dL/dt was discussed only for the radial phase.In our opinion the more important phase for the purpose of neutron saturation is actually the axial phase for the Mather-type plasma focus.
We now resolve the problem into its most basic form as follows.We have a generator (the capacitor charged to 30 kV), with an impedance of Z 0 = (L 0 /C 0 ) 0.5 driving a load with a near constant resistance of 7 m.We also assign a value for stray resistance of 0.1Z 0 .This situation may be shown in Table 3 where L 0 is given a typical value of 30 nH.We also include in the last column the results from a circuit (L-C-R) computation, discharging the capacitor with initial voltage 30 kV into a fixed resistance load of 7 msimulating the effect of the DR 0 and a stray resistance of value 0.1Z 0 .Table 3. Discharge characteristics of equivalent PF circuit, illustrating the 'saturation' of I peak with increase of E 0 to very large values.The last column presents results using circuit (L-C-R) computation, with a fixed resistance load of 7 m, simulating the effect of the DR 0 and a stray resistance of value 0.1Z 0 .Plotting the peak current as a function of E 0 we obtain Figure 14, which shows the tendency of the peak current towards saturation as E 0 reaches large values; the deterioration of the curve becoming apparent at the several hundred kJ level.This is the case for I peak = V 0 /Z total and also for the L-C-R discharge with simulated value of the DR 0 .In both cases it is seen clearly that a capacitor bank of voltage V 0 discharging into a constant resistance such as DR 0 will have a peak current I peak approaching an asymptotic value of I peak = V 0 /DR 0 when the bank capacitance C 0 is increased to such large values that the value of Z 0 = (L 0 /C 0 ) 0.5 << DR 0 .Thus DR 0 causes current 'saturation'.In Section 4.2 we had shown the following relationships between Y n and I peak and I pinch as follows: Y n ~ I pinch 4.5 (8) Y n ~ I peak 3.8 (9) Hence saturation of I peak will lead to saturation of Y n .At this point we note that if we consider that only 0.7 of the total current takes part in driving the current sheet, as typically agreed upon from experimental observations, then there is a correction factor which reduces the axial dynamic resistance by some 40%.That would raise the asymptotic value of the current by some 40%; nevertheless there would still be 'saturation'.
Thus we have shown that current 'saturation' is inevitable as E 0 is increased to very large values by an increase in C 0 , simply due to the dominance of the axial phase dynamic resistance.This makes the total circuit impedance tend towards an asymptotic value which approaches the dynamic resistance at infinite values of E 0 .The 'saturation' of current inevitably leads to a 'saturation' of neutron yield.Thus the apparently observed neutron 'saturation' which is more accurately represented as a neutron scaling deterioration is inevitable because of the dynamic resistance.In line with current plasma focus terminology we will continue to refer to this scaling deterioration as 'saturation'.The above analysis applies to the Mather-type plasma focus.The Filippov-type plasma focus does not have a clearly defined axial phase.Instead it has a lift-off phase and an extended pre-pinch radial phase which determine the value of I peak .During these phases the inductance of the Filippov discharge is changing, and the changing L(t) will develop a dynamic resistance which will also have the same current 'saturation' effect as the Filippov bank capacitance becomes big enough.

Beyond Presently Observed Neutron Saturation Regimes
Moreover the 'saturation' as observed in presently available data is due also to the fact that all tabulated machines operate in a narrow range of voltages of 15−50 kV.Only the SPEED machines, most notably SPEED II [51,54] operated at low hundreds of kV.No extensive data have been published from the SPEED machines.Moreover SPEED II, using Marx technology, has a large bank surge impedance of 50 mwhich itself would limit the current.If we operate a range of such high voltage machines at a fixed high voltage, say 300 kV, with ever larger E 0 until the surge impedance becomes negligible due to the very large value of C 0 .then the 'saturation' effect would still be there, but the level of 'saturation' would be proportional to the voltage.In this way we can go far above presently observed levels of neutron 'saturation'; moving the research, as it were into presently beyond-saturation regimes.
Could the technology be extended to 1 MV?That would raise I peak to beyond 15 MA and I pinch to over 6 MA.Also multiple Blumleins at 1 MV, in parallel, could provide driver impedance of 100 m, matching the radial phase dynamic resistance and provide fast rise currents peaking at 10 MA with I pinch value of perhaps 5 MA.Bank energy would be several MJ.The push to higher currents may be combined with proven neutron yield enhancing methods such as doping deuterium with low % of krypton [55].Further increase in pinch current might be by fast current injection near the start of the radial phase.This could be achieved with charged particle beams or by circuit manipulation such as current-stepping [31,33,56−57].The Lee model is ideally suited for testing circuit manipulation schemes.

Conclusions
This paper has reviewed the extensive and systematic numerical experiments which have been used to uncover new insights into plasma focus fusion devices including the following.A plasma current limitation effect was unexpectedly found, as the static inductance of any focus device is reduced towards very small values.Scaling laws of neutron yield and soft X-ray yield as functions of storage energies, circuit peak current as well as plasma pinch current, were developed over wider range of parameters than attempted previously.A global scaling law for neutron yield as a function of storage energy was uncovered combining experimental and extensive numerical data, showing that scaling deterioration has been wrongly interpreted as neutron 'saturation'.However in keeping with conventional terminology, the effect of scaling deterioration will continue to be referred to as neutron 'saturation'.The cause of neutron 'saturation' as device storage energy is increased was found to be the axial phase 'dynamic resistance'.With the fundamental cause discovered, it is suggested that beyond 'present saturation' regimes may be reached by going to higher voltages, and using plasma current enhancement techniques such as current-steps.It is expected that numerical experiments will continue to play a major role complementing and even guiding laboratory measurements and practices.The ground-breaking insights thus gained will completely open up the directions of plasma focus fusion research.

Figure 2 .
Figure 2. Schematic of radius versus time trajectories to illustrate the radial inward shock phase when r s moves radially inwards, the reflected shock (RS) phase when the reflected shock moves radially outwards, until it hits the incoming piston r p leading to the start of the pinch phase (t f ) and finally the expanded column phase.

Figure 3 .
Figure 3.The 5-point fitting of computed current trace to the measured (or the reference) current trace.Point 1 is the current rise slope.Point 2 is the topping profile.Point 3 is the peak value of the current.Point 4 is the slope of the current dip.Point 5 is the bottom of the current dip.Fitting is done up to point 5 only.Further agreement or divergence of the computed trace with/from the measured trace is only incidental and not considered to be important.

Figure 4 .
Figure 4. Computed plasma and electrodynamic properties of the NX2.

Figure 5 .
Figure 5. Fitting 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.The measured current trace was for the PF1000 at 27 kV, storage capacity of 1332 F and fitted static inductance of 33.5 H.

Figure 6 .
Figure 6.PF1000 current waveforms computed at 35 kV, 3.5 Torr D 2 for a range of L 0 showing the changes in waveforms as L 0 varies.

Figure 8 .
Figure 8. Y n plotted as a function of E 0 in log-log scale, showing Y n scaling changes from Y n ~ E 0 2.0 at tens of kJ to Y n ~ E 0 0.84 at the highest energies (up to 25 MJ).

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

Figure 11 .
Figure 11.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 f m , f c , f mr , f cr at 0.06, 0.7, 0.16 and 0.7 respectively.The white data points are for specific machines which have different values for the parameters c, L 0 and V 0 .

Figure 12 .
Figure 12.The global scaling law, combining experimental and numerical data.The global data illustrates Y n scaling deterioration observed in numerical experiments from 0.4 kJ to 25 MJ (solid line) using the Lee model code, compared to measurements compiled from publications (squares) of various machines from 0.4 kJ to 1 MJ.

Figure 13 .
Figure 13.Plasma focus circuit schematic.The capacitor bank with static inductance L 0 and stray resistance r 0 is switched into the plasma focus tube where a fraction f c of the circuit current I(t) effectively drives the plasma creating a time-varying inductance L(t) in the focus tube.

Figure 14 .
Figure 14.I peak versus E 0 on log-log scale, illustrating I peak 'saturation' at large E 0 .

Table 1 .
Effect on currents and ratio of currents as L 0 is reduced-PF1000 at 35kV, 3.5 Torr Deuterium.