An ion species model for positive ion sources: I. Description of the model

Before plasma transport is examined, the primary flux and energy flow needed to sustain the plasma is required. The plasma modelling is divided into three parts, the primary electron input, the plasma itself and finally the plasma sheath at the plasma grid, which normally floats at a negative potential due to the influx of primaries. In most sources there is zero current to this grid as the ion and total electron currents cancel. The code however retains the possibility of applying an external potential bias and positive currents are defined as net ion flow to the grid. A schematic of a generic source and the potentials discussed below is shown in figure 1.

Zoom In Zoom Out Reset image size

The plasma itself is biased by a potential, φ, positive with respect to the anode (ion source outer wall). The plasma electrons and the degraded primaries have two loss routes: to the anode and across the filter and each electron removes an energy eT₁(T₁ is the driver plasma temperature) as it escapes. The model assumes a single loss route to the wall and this argument is justified by the existence of a potential difference across the filter region which causes the hot thermal electrons to lose almost all their energy, eT₁. This avoids having to implement a spatial variation of energy removal as the electron temperature decreases across the filter within the source. There is an energy balance for the main plasma where the driver thermal electron temperature is determined by the input energy to the thermal plasma derived from the energy transfer from the fast primary electrons. Current balance at the anode determines the plasma potential, φ which determines the energy removed by the escaping plasma ions.

The ions escape with an energy eφ if they are collected on the anode and a much larger value, e(φ + V_g) if they are collected on the plasma grid. Note that V_g is defined as a positive value when plasma grid is negatively biased with respect to theanode. The net current to the plasma grid is a function of the combined energy distribution of the thermal and primary electrons. The potential energy, V_g + φ, lost by both groups of electrons in crossing the sheath in front of the plasma grid is transferred to the ions so they can escape with an equal energy. For this reason, the plasma grid net current is input data to the model and is used to determine the grid potential, V_g.

The primary electron power is I_e(V − V_fil) where I_e is the primary emission current (close but not equal to the arc current) and is an input variable to the code, V is the arc voltage and V_fil is the drop in voltage across the filament. For simplicity, we ignore the fact that the primaries are slightly more energetic by the plasma potential, eφ. This is justified as φ is typically 2 to 5 V and is very small compared with the arc voltage. This energy is either transferred to the plasma electrons or is lost through inelastic collisions using the value developed by Hiskes and Karo [12] (which includes ionization). The primaries then escape with just the thermal temperature, T₁, if degraded by inelastic collisions but a few (a fraction, k_f) escape with their full energy if they find a leakage path from the plasma. The two leakage areas are either the anode or the plasma grid but the exact leak rates are not known at the start of the calculation and are determined by iteration.

There is a final current balance at the grid. The primaries and plasma electrons are retarded by the grid potential, V_g, so that the algebraic sum of the electron and ion currents is equal to the grid current. The ion flux in turn is defined by the ionization rate less any losses to anode or filament and matches the ion flux surviving the transit across the plasma. However, again, the ion leaks are not known at the start of the calculation and are found by iteration. There can be a situation where an external amplifier or power supply is used to control the grid bias. In this instance, an extra power input of I_gV_g is added, where I_g can be positive or negative (the latter occurs when V_g is less than the floating potential where I_g is zero). This extra power is added to the ion energy flux at the grid. This is described in the section on the plasma grid.

In order to solve these equations some boundary conditions must be set. The initial primary temperature (eV) is mainly determined by the arc voltage, V, less the voltage drop across the filament, V_fil and is:

$$\begin{equation} T_{\rm f0}=2(V-V_{\rm fil})/3. \end{equation} \tag{ 3 }$$

The plasma potential is ignored here as it is small compared with the filament drop. The 2/3 factor arises from the kinetic theory of gases relating temperature to total energy and assumes that the ballistic nature of the primary electrons as they are emitted from the filament is modified into a distribution close to a Maxwellian; experiment indicates that this is applicable at pressures as low as ∼0.5 Pa [17]. The primary electron flux, n_f0v_f, is found from the emission current, I_e, (this is not the cathode or arc current) less any electron emission current collected by the anode. The initial primary density is given by an equation similar to equation (11) in Holmes [9].

$$\begin{eqnarray} &&\fl\left( {I_{\rm e}-k_{\rm f}I_{\rm fg}} \right)/e=k_{\rm f}A_{\rm a}n_{\rm f0}v_{\rm f}/4\nonumber\\ &&+\left( {n_{\rm f0}NR_{\rm in}/T_{\rm f}+2n_1k\lambda T_{{\rm f}}^{-3/2}} \right)A_{{\rm g}} D_{{\rm s}} , \end{eqnarray} \tag{ 4 }$$

where I_fg is the fast electron current to the plasma grid of area A_g. Note that I_fg is a positive quantity. The second term in the bracket is the Coulomb drag loss, which heats the thermal density, n₁. The term A_a is the anode cusp loss area and R_in is the inelastic energy loss rate from Hiskes and Karo [12]. The parameter, k_f, represents the effectiveness of the lost primaries in imparting energy through inelastic collisions (i.e. if k_f = 1 then these electrons do not have any inelastic processes). This rate, R_in, is:

$$\begin{equation*} R_{\rm in}=2.4\times 10^{-12}\exp \left( {-28/T_{\rm f0}} \right)\quad \hbox{eV\,m$^{3}$\,s$^{-1}$}. \end{equation*}$$

The rate, R_in, can be converted to a simple collision rate by dividing by the primary temperature. A simple expression for k_f can be derived by comparing the characteristic time for primary electron loss to the walls, τ_w with the inelastic loss collision time. If the former is much shorter than the latter, then k_f approaches unity while the reverse leads to k_f nearing zero. A suitable function would employ the exponential form, so that:

$$\begin{equation*} k_{\rm f}=\exp \left( {-\nu \tau_ {\rm w}} \right)=\exp \left( {-\frac{NR_{\rm in}D_{\rm s}}{v_{\rm f0}T_{\rm f0}}} \right). \end{equation*}$$

Hershkowitz et al [18] have argued that the plasma loss to the anode is equal to the total cusp length, C, multiplied by twice the hybrid Larmor diameter. This argument is extended to the energetic primary electrons, but as they are decoupled from the plasma electrons by their much higher energy, the primary electron loss area to the anode is twice the primary Larmor diameter multiplied by the total cusp length, C. Thus:

$$\begin{equation*} A_{\rm a}=\frac{4Cm}{eB_{\rm c}}\left( {\frac{8eT_{\rm f0}}{\pi m}} \right)^{1/2}. \end{equation*}$$

Unfortunately the solution of (1) and (2) is very prone to instability, particularly as F_f and Q_f are both functions of the distance. A more practical method is to seek a solution that is consistent with (1) and (2) and retains a realistic attenuation of the current density and energy flux. Firstly (1) and (2) are rewritten in a more compact form where the subscripts have been dropped:

$$\begin{equation} f=F_{{\rm f}} =-aTn'+0.5bnT' \end{equation} \tag{ 5 }$$

$$\begin{equation} q=-Q_{{\rm f}} /T=aTn'+1.92bnT'. \end{equation} \tag{ 6 }$$

The terms a and b are defined as:

$$\begin{equation*} a=\frac{e}{m\upsilon _{\rm f1}\left( {1+\left( {\omega /\upsilon_ {\rm f1}} \right)^2} \right)} \quad b=\frac{e}{m\upsilon_{{\rm f2}} \left( {1+\left( {\omega /\upsilon_{{\rm f}2}} \right)^2} \right)}. \end{equation*}$$

Note that f is a positive and q is a negative quantity and that n', T' are also negative. The term a is positive but b is negative if Coulomb collisions are small due to the high value of the primary temperature. A solution is needed where n and T decrease in an assumed form of exponential decay with distance from the filament at z = 0. One possibility is:

Let

$$\begin{equation} n=n_0\exp \left( {-\int {\frac{\alpha}{a}\,{\rm d}z}} \right) \end{equation} \tag{ 7 }$$

and

$$\begin{equation} T=T_0\exp \left( {-\int {\frac{\beta}{b}\,{\rm d}z}} \right). \end{equation} \tag{ 8 }$$

After substituting (7) and (8) into (1) and (2) and applying the boundary conditions, the solution for the density and temperature is:

$$\begin{equation} n_{{\rm f}} =n_{{\rm f0}} \exp \left( {-\int {\frac{0.234\upsilon_ {\rm f1}\;\left( {1+\left( {\omega /\upsilon_ {\rm f1}} \right)^2} \right)}{\sqrt {eT_{\rm f}/m}}\,{\rm d}z}} \right) \end{equation} \tag{ 9 }$$

and

$$\begin{equation} T_{\rm f}=T_{\rm f0}\exp \left( {\int {\frac{0.33\upsilon_{{\rm f2}} \left( {1+\left( {\omega /\upsilon_{{\rm f}2}} \right)^2} \right)}{\sqrt {eT_{\rm f}/m}}\,{\rm d}z}} \right). \end{equation} \tag{ 10 }$$

Note that υ_f2 is negative in (10), so both n_f and T_f decrease with distance from the start line (z = 0) as expected. The total ionization current, I_i, is:

$$\begin{equation} I_{\rm i}=\int_0^{{\rm d}s} {n_{\rm f}(N_2+x_0N_1)q_{\rm i}(T_{\rm f})\,{\rm d}z-I_{\rm rec}} , \end{equation} \tag{ 11 }$$

where N₂ is the molecular gas density, N₁ is the atomic gas density, q_i(T_f) is the molecular ionization rate (see section 6), x₀ is the multiplier for atomic ionization (x₀ ≈ 0.66 ± 0.04) and I_rec is the ionization current lost by recombination of ${\rm H}_{3}^{+}$ with H to reform H₂ molecules (see section 8.1). Note that the atomic ionization rate is almost exactly 0.66 of the molecular rate at all energies slightly above the threshold, so for convenience, only the molecular rate is used. This I_rec current is not known in the first cycle but can be derived from the build-up of the ion species described in section 8.1.

The primary effectiveness, k_f is initially set to zero and hence a value for n_f0 can be derived from (4) and the arc voltage gives a value for T_f0. This latter value does not change in later cycles. As T_f appears in the argument of the exponent in (10), care is required for step-wise integration because of the changing value of ω, the cyclotron frequency. Differentiating (10) gives:

$$\begin{equation*} T'_{{\rm f}} =T_{\rm f}\frac{0.3298\upsilon _{\rm f2}\left( {1+\left( {\omega /\upsilon_ {\rm f2}} \right)^2} \right)}{\sqrt {eT_{\rm f}/m}}. \end{equation*}$$

The derivative of T_f varies with the local square root of T_f, which can be integrated (where ω is a local constant) to get:

$$\begin{equation*} T_{\rm fn}=\left[ {\frac{0.3298\upsilon_ {\rm f2}\left( {1+\left( {\omega /\upsilon_{\rm f2}} \right)^2} \right)}{2\sqrt {e/m}}\delta +\sqrt {T_{{\rm fn}-1}}} \right]^{2}, \end{equation*}$$

where δ is the local increment along the z axis from the (n − 1) position. Once the new value of T_f is obtained, an average root value is made over the step and used in (9) to obtain the change in n_f. At the same time the local ionization current, I_i, is obtained by multiplying by the ionization rate and the density of atomic and molecular gas.

Once a solution is found for the plasma density in the driver region, the coulomb drag element can be included as a correction in subsequent cycles and a new value of the coefficient, k_f, can be assigned.

The next step is to derive the atomic gas fraction. For simplicity no negative ion rates are included as these are not expected to impact significantly on the atomic gas fraction. For hydrogen, the ionization rate for atomic hydrogen is 0.66 of the molecular rate at all energies within an error of ±5% (Chan [14]). The atomic fraction is assumed to be constant everywhere inside the source, so the calculation is global for the entire source chamber. This is particularly true for atomic loss on the walls, which is obviously a global process. The complicating factor is the contributions made to protons and atomic hydrogen by ${\rm H}_{3}^{+}$ dissociation and recombination. Unlike the other processes these are local and mainly take place away from the filament region.

To circumvent this problem, these processes will be treated as perturbations that are initially zero and introduced as additional source or sink terms in later cycles. The processes that contribute are from Chan [14], improved from Green [19] except the ionization rate which is from Freeman and Jones [20].

The rates are introduced as algebraic empirical expressions, shown in table 2, that fit the experimental data over the temperature range up to about 100 eV for ionization and other energetic collisions. For low energy collision processes it is valid up to about 20 eV. This avoids problems with polynomial fits which diverge if the solution process goes out of range. The only exception is the xs2 rate which is fitted by two different expressions above or below T_f = 22 eV.

Table 2. Analytic forms of reaction rates used in the model.

Process	label	Expression in m3 s−1	Reference
H2 + e → 2 H + e	xs2	xs2 = 10−15exp(3.2)exp(−8/T)F	[14]
		if T > 22, then F = 1
		else F = exp(−0.068(T − 22))
${\rm H}_{2}^{+}+\rm e \to 2\,H$	xs4	xs4 = 10−15exp(3.5)(1 − exp − (5.2/T))	[14]
${\rm H}_{2}^{+}+\rm e \to H^{+} + H + e$	xs5	xs5 = 10−15exp(4.5)exp(−2.5/T)	[14]
${\rm H}_{2}^{+} +\rm H_{2}\to {\rm H}_{3}^{+} + H$	xs6	2.1 × 10−15	[14]
${\rm H}_{3}^{+} +\rm e \to H^{+} + H_{2} + e$	xs7	xs7 = 10−15exp(6.5)exp(−13/T)	[14]
${\rm H}_{3}^{+} +\rm e \to H_{2} + H$	xs8	xs8 = 10−15exp(3.8)(1 − exp(−6/T))	[14]
${\rm H}_{2} +\rm e \to {\rm H}_{2}^{+} +2e$	qi	qi = 10−15exp(3.98)exp(−22.84/T)	[19]

Chan et al [21] has derived a particle balance equation that can be used to derive the atomic density, N₁. In the first order any process that contains ${\rm H}_{3}^{+}$ is ignored to obtain:

$$\begin{eqnarray*}&&\fl\frac{N_1\gamma}{t_1}+0.66N_1n_{\rm f0}^ \ast q_{{\rm i}} =2N_2\left( {n_{\rm e}^\ast xs2_{\rm e}+n_{\rm f}^\ast xs2_{\rm f}} \right)\nonumber\\ &&\fl\quad+n_{{\rm i}2}^ \ast n_{\rm e}^\ast \left( {xs5+2\times xs4} \right) +N_2n_{{\rm i2}} ^\ast xs6+n_{{\rm i}3}^{\ast} n_{\rm e}^\ast xs8. \end{eqnarray*}$$

Here N₂ is the molecular gas density, n₂ is the ${\rm H}_{2}^{+}$ ion density and n_e is the plasma density. The term γ is the accommodation coefficient on the walls (Chan [21] selected a value of 0.3 for γ). The term, t₁, is the transit time of an atom in the source (= 4D/v_atom). The asterisk indicates that an average is taken of the densities across the source depth (i.e. initial/2+final/2) to allow for the spatial distribution.

In addition we have the condition that the total mass flow, Q_g, of gas through the source is constant, so:

$$\begin{eqnarray} &&\fl Q_{\rm g}=2m_{\rm p}\,N_0\kappa T_{{\rm w}}^{1/2}=2m_{\rm p}\,N_2\kappa T_{\rm g}^{1/2}+\,m_{\rm p}\,N_1\sqrt 2 \kappa T_{\rm g}^{1/2}.\nonumber\\ \end{eqnarray} \tag{ 17 }$$

Here the initial gas input is at the wall temperature, denoted by N₀ and T_w but emerges into the accelerator as hot gas of temperature, T_g, with densities N₁ and N₂. The term κ is the fixed part of the gas conductance and includes the number and area of the extraction apertures.

If N₁/N₂ = G then rationalization gives:

$$\begin{equation} \sqrt {\frac{T_{\rm w}}{T_{\rm g}}} N_0=N_2\left( {1+G/\sqrt 2} \right)=N_{{\rm c}} . \end{equation} \tag{ 18 }$$

Here N_c is the gas density in the presence of the discharge (i.e. the gas is hot) but without dissociation. Note that N₁/2 + N₂ is not equal to N_c because of the increased mobility of the atomic gas.

Each of the ion densities is related to the plasma density, n_e, by a coefficient, c_n1, c_n2 or c_n3:

$$\begin{eqnarray} &&n_{{\rm i1}} =n_{\rm e}c_{{\rm n1}}\\ &&n_{{\rm i}2} =n_{\rm e}c_{\rm n2}\\ &&n_{{\rm i}3} =n_{{\rm e}} c_{{\rm n3}} . \end{eqnarray}$$

The values of the molecular ion fractions, c_n2 and c_n3, are not known initially so a preliminary guess of 0.333 is assumed but this is updated after each major cycle of convergence.

The values of n_e, T_f0 and T_e are found from the solution of the driver region giving all rates. The atomic gas velocity is derived from an assumed initial gas temperature of 500 K. The initial value of n_f0 is known and the undisturbed gas density, N₀, is obtained from the initial gas temperature. Simplifying the particle balance equations developed by Chan et al [21] using the coefficients below:

$$\begin{eqnarray*}&&\fl {R}'=\frac{\gamma}{t_1}+0.66n_{{\rm f0}}^{{\rm \ast}} q_{{\rm i}} \quad {H}'=2\left( {n_{{\rm e}}^{{\rm \ast}} xs2_{\rm e}+n_{{\rm f0}}^{{\rm \ast}} xs2_{\rm f}} \right) \\ &&\fl {B}'=c_{\rm n2}n_{{\rm e}}^{{\rm \ast \, 2}}\left( {xs5+2xs4} \right) \quad {C}'=c_{\rm n2}n_{{\rm e}}^{{\rm \ast}} xs6 \\ &&\fl {D}'=c_{\rm n3}n_{{\rm e}}^{{\rm \ast \,2}}xs8. \end{eqnarray*}$$

The balance equation for molecular gas is then:

$$\begin{equation*} GN_2{R}'={H}'N_2+{B}'+{C}'N_2+{D}'. \end{equation*}$$

Substituting N₂ with N_c using (18) gives after rationalization:

$$\begin{equation} G=\frac{{H}'+{C}'+\left( {{B}'+{D}'} \right)/N_{{\rm c}}}{{R}'-\left( {{B}'+{D}'} \right)/(\sqrt 2 N_{{\rm c}} )}. \end{equation} \tag{ 19 }$$

Once G is known, N₂ can be found and hence:

$$\begin{equation} N_1=N_{{\rm c}} /(G^{-1}+1/2^{0.5}). \end{equation} \tag{ 20 }$$

There are several channels that contribute to the gas temperature such as electron–molecule collisions, Coulomb drag, ion–molecular collisions arising from the electric field and hot molecules arising from backscattering of ions accelerated across the plasma sheath. At high arc currents, the Coulomb drag term dominates due to its quadratic scaling with plasma density. At low arc currents, the gas temperature tends towards the source wall temperature as the plasma density is low so the impact of the other heating channels is weak and can be neglected for simplicity. For the case where the ions are predominantly heated by energy transfer from the electrons by Coulomb collisions and using the driver region conditions where the plasma density is highest, the energy balance is:

$$\begin{eqnarray*}&&\fl eA_{{\rm w}} \gamma_{{\rm h}} \left( {T_{\rm g}-T_{\rm w}} \right)Nv_{\rm gas}/4=Q_{\rm i}\\ &&=k_{\varepsilon} eA_{{\rm g}} D_{\rm s}\lambda n_1^2T_1^{-1/2}/a_{{\rm mass}} . \end{eqnarray*}$$

Here T_w is the wall temperature in eV, T_g is the hot gas temperature, k_ε = 3.2 × 10⁻¹⁵ (eV)^3/2 m³ s⁻¹ is the Coulomb rate coefficient for energy equilibration between ions and electrons [15], λ is the Coulomb logarithm (= 10), A_w is the total gas cooling surface area, A_gD_s is the grid area times plasma depth. The term γ_h is the thermal accommodation coefficient and is probably similar to the atomic recombination coefficient and v_gas is the gas velocity. However not all the ions have charge exchanging collisions before hitting the anode so an additional reduction of 1 − exp(−NσD_s) is included in the numerator on the rhs. The product Nv_gas is conserved as the gas flow is constant and this simplifies to:

$$\begin{eqnarray} &&\fl T_{\rm g}=T_{\rm w}+\frac{4A_{{\rm g}} d_{\rm s}k_{\varepsilon} n_1^2T_1^{-1/2}\left[ {1-\exp \left( {-N\sigma D_{{\rm s}}} \right)} \right]}{A_{{\rm w}} a_{{\rm i}} \gamma_{{\rm h}} N_0\left( {\frac{2.55eT_{\rm w}}{2m_{{\rm p}}}} \right)^{1/2}}.\nonumber\\ \end{eqnarray} \tag{ 21 }$$

The gas density scales as (T_w/T_g)^1/2 as the gas mass flow is constant as shown in (18). This is used to update the value of k_f at the end of a cycle.

In the above analysis two wall coefficients have been used; γ, the atomic hydrogen recombination rate and γ_h, the thermal accommodation rate. In the work by Chan et al [21], only the first was used and was treated as a fixed constant. The values are not known with any certitude and can only be determined by comparing the code results with experimental results. However Chan, et al argue that γ should be about 0.2 for deuterium discharges although experimental data from Wood and Wise [22] suggest it should a little smaller. In the results presented here the values given in table 3 are used.

The values of γ are closer to the results of Wood and Wise [22] but it should be emphasized that the presence of the plasma and the evaporated tungsten on the copper walls makes it difficult to have a definitive value. A sensitivity study varying these parameters over the range 0.1 to 0.4 showed a modest influence on the species ratios for both with the stronger effect being a decrease in deuteron yield with increasing γ but an increase in yield with γ_h.

The directed initial current density of ions and electrons towards the plasma grid are zero at the plane of origin at the back of the source as no ionization has yet taken place. By analogy, the directed electron energy flux is also zero. The ionic grid current is defined as a fraction, f_i, of the total ion production and the electron current is a further fraction, f_e of this current. The former has an initial value close to unity while the latter is close to zero.

In the equations below the diamond bracket indicates a simple average of the flux product taken over the previous step (the first step assumes that values are zero). The velocities in the denominator are the mass weighted thermal ion velocities based on the gas temperature which is also the effective ion temperature due to the high ion–neutral collisionality. The term x₀ is the ratio of ionization rates for atomic and molecular gas and is 0.66. After each step the total fluxes of all ion species are re-evaluated by integrating the following expressions from the source backplate to the plasma grid:

$$\begin{eqnarray} &&\fl \frac{{\rm d}F_1}{{\rm d}z}=n_{\rm f}N_{1} x_0q_{\rm i}+\frac{\left\langle {F_2n_{\rm e}} \right\rangle xs5}{v_2}\nonumber\\ &&\quad\quad+\,\frac{\left\{{\left\langle {F_3n_{\rm e}} \right\rangle xs7_{\rm e}+\left\langle {F_3n_{\rm f}} \right\rangle xs7_{\rm f}} \right\}}{v_3} \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} &&\fl \frac{dF_2}{dz}=n_{\rm f}N_{2} q_{\rm i}-\frac{\left\langle {F_2n_{\rm e}} \right\rangle xs5}{v_2}-N_2\left\langle {F_2} \right\rangle xs6 \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} &&\fl \frac{{\rm d}F_3}{{\rm d}z}=N_2\left\langle {F_2} \right\rangle xs6-\frac{\left\{{\left\langle {F_3n_{\rm e}} \right\rangle xs7_{\rm e}+\left\langle {F_3n_{\rm f}} \right\rangle xs7_{\rm f}} \right\}}{v_3}\nonumber\\ &&\quad\quad-\,\frac{\left\langle {F_3n_{\rm e}} \right\rangle xs8}{2v_3}. \end{eqnarray} \tag{ 29 }$$

After integration from the backplate to the plasma grid, the total directed ion flux, F_ig, is composed of the three individual fluxes, F_1g, F_2g and F_3g:

$$\begin{equation} F_{{\rm ig}} =\left( {F_{1{\rm g}} +F_{2{\rm g}} +F_{{\rm 3g}}} \right) \end{equation} \tag{ 30 }$$

and:

$$\begin{equation*} C_{1} =F_{1{\rm g}} /F_{{\rm ig}} \quad C_{2} =F_{2{\rm g}} /F_{{\rm ig}} \quad C_{3} =F_{{\rm 3g}} /F_{{\rm ig}}. \end{equation*}$$

This integration is made after the electron spatial distribution has been found using (24) to (26) which is based on the ion current density fractions, C₁, C₂ and C₃, derived in an earlier cycle. The new integration allows an update of these fractions to be made. The values of the plasma density fractions, c_n1, c_n2 and c_n3 are derived by dividing the current density fractions by the normalized velocities appropriate to that ion.

The fraction of the total ion current that goes to the plasma grid is only part of the total ion production, typically ∼60%, with the rest going to the anode or the filament structure. The addition of a magnetic filter increases the latter currents significantly and reduces the former current. These new fractions are used in the entire transport calculations (26) to define that average ionic transport coefficient for the cycle that follows the definition of the C₁, C₂ and C₃ values.

In the case of electron flow, it is assumed that the electron flux, F_e, is a fraction, f_e, of the total ion flux that is collected by the plasma grid. Hence:

$$\begin{equation*} F_{\rm e}=f_{\rm e}I_{{\rm ig}} =f_{{\rm e}} f_{{\rm i}} e\left( {F_{{\rm 1g}} +F_{2{\rm g}} +F_{3{\rm g}}} \right), \end{equation*}$$

where

$$\begin{equation*} I_{{\rm g}} =f_{{\rm i}} I_{{\rm i}} \left( {1-f_{{\rm e}}} \right). \end{equation*}$$

The electron current to the plasma grid is a low fraction of the total electron production in almost all situations. The values of f_e and f_i are found by iteration over many cycles so that it matches the assumptions about the plasma grid current, I_g, that form part of the input data. This is discussed further in section 8.2. However it is necessary to guess the values initially but these are updated in later cycles.

If all the three ion flux gradients (27) to (29) are added, all the terms cancel apart from ionization (the q_i terms) and the xs8 term in the gradient of F₃. This last term is the recombination of ${\rm H}_{3}^{+}$ and H₂ and represents a loss of ionization. Thus:

$$\begin{equation} I_{\rm rec}=A_{{\rm g}} \int_0^{{\rm d}s} {\frac{e\left\langle {F_3n_{\rm e}} \right\rangle xs8}{2v_3}} \,{\rm d}z. \end{equation} \tag{ 31 }$$

This current is included in the total ion current balance but is usually quite small.

The random ion current density controlled by the ion sound speed that reaches the plasma grid of area A_g must equal the directed ion current density that is transported via (26). This can be expressed by the equality:

$$\begin{equation} ef_{{\rm i}} F_{{\rm ig}} =I_{{\rm ig}} /A_{{\rm g}} =en_{{\rm ig}} \psi_ {\rm i}T_{\rm eg}^{1/2}. \end{equation} \tag{ 32 }$$

Here n_ig is the ion density at the grid following a solution of the equations in the transport section 8 and T_eg is the matching temperature. Equation (32) creates a new value of the ion fraction, f_i based on the value of n_ig and T_eg and the value of F_ig derived in section 8.1. The value of f_i, used in the next cycle, is created by mixing the original estimate with the new value derived from (32) above. Typically between 5% and 10% of the new solution is used, giving an under-relaxation parameter of around 0.1.

The other quantities to be derived are the plasma grid bias voltage, V_g and the value of f_e. The former is found from current continuity:

$$\begin{eqnarray} &&\fl A_{\rm g}e\Bigg[ \psi_{{\rm e}} n_{\rm fg}T_{\rm fg}^{1/2}\exp \left( {-\frac{V_{\rm g}}{T_{\rm fg}}} \right)\nonumber\\ &&+\,\psi_{{\rm n}} n_{\rm eg}T_{\rm eg}^{1/2}\exp \left( {-\frac{V_{\rm g}}{T_{\rm eg}}} \right) \Bigg]=I_{{\rm eg}} +I_{{\rm fg}} =I_{\rm ig}-I_{{\rm g}} .\nonumber\\ \end{eqnarray} \tag{ 33 }$$

The negative ions are assumed to behave as thermal electrons but having a mass, M_n, equal to the atomic ion mass at the sheath, hence the composite value of ψ_n for electrons and negative ions can be derived as:

$$\begin{equation*} \psi_{{\rm n}} =\frac{1}{4}\left( {\frac{8e}{\pi m}} \right)^{1/2}\left( {1-q} \right)+\frac{1}{4}\left( {\frac{8e}{\pi M_{{\rm n}}}} \right)^{1/2}\left( q \right). \end{equation*}$$

The actual grid current is I_g, defined as positive if flowing conventionally from the grid (i.e. very negative bias of the plasma grid creates a positive value for I_g). This last equation is highly non-linear and must be solved using a step-wise procedure. The values of n_fg and T_fg are derived from the last primary electron distribution. If the filter is moderately strong there is no primary current. There is no explicit negative ion current at the plasma grid used in (33) as the negative ion velocity is small relative to that of the ions, but their presence is allowed for through the fact that the electron density is smaller than that of the positive ions.

The value of V_g is used to derive an update for the electron fluxes via the expression:

$$\begin{equation} f_{{\rm e}} I_{{\rm ig}} /A_{\rm g}=e\psi_ {\rm e}n_{\rm eg}T_{\rm eg}^{1/2}\exp \left( {-\frac{V_{\rm g}}{T_{\rm eg}}} \right) \end{equation} \tag{ 34 }$$

$$\begin{equation} I_{\rm fg}=eA_{\rm g}\psi_ {\rm e}n_{\rm fg}T_{\rm fg}^{1/2}\exp \left( {-\frac{V_{\rm g}}{T_{\rm fg}}} \right). \end{equation} \tag{ 35 }$$

Again the actual negative ion current density is ignored as it will be small compared with that of the electrons. Similarly to the convergence of the value of f_i, severe under-relaxation is needed for f_e. These new values of f_i, f_e and I_fg are returned to the primary modelling described in sections 4 and 4.1 together with the recombination current given by (31). Typically it takes about one thousand cycles to converge to better than 1% although the convergence rate depends strongly on gas pressure.

This is a cubical magnetic multipole source (of 10 inch dimensions, hence the name) with a filter formed by small bar magnets inside water cooled tubes placed inside the plasma approximately 95 mm from the plasma grid, creating a dipole filter. By varying the magnets, Pincosy [7] created several magnetic filter geometries. The ion species were measured by a magnetic momentum analyser and the plasma properties were derived from Langmuir probes. Figure 4 shows the source species for a deuterium discharge as a function of ion current density at the extraction plane without a magnetic filter and for comparison the output from the model of this source.

Zoom In Zoom Out Reset image size

The addition of a magnetic filter to give integrated flux of 1.05 × 10⁻⁴ or 2.1 × 10⁻⁴ T m increases the deuteron fraction as shown in figures 5 and 6, again with the results of the model for comparison. Figures 4 to 6 indicate qualitative agreement between experiment and the model, which is able to reproduce the species variations observed in the discharge. There is a tendency to over-estimate the improvement in the deuteron yield with increasing filter field but the main trend of the improvement in deuteron fraction is clear. The deuteron fraction increases from zero discharge current density and the rate of rise and its limiting value at high arc currents follow the experimental data closely.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The filter field increases the plasma density in the driver region where the filaments are situated and inhibits the plasma flow to the extraction plane, hence the truncation of the data at lower current density in figure 6. The primary electrons cross the filter field by collisions, losing energy in the process, which serves to reduce the plasma temperature. These effects are illustrated in figure 7 for a 1.05 × 10⁻⁴ T m filter field in a 35 kW deuterium discharge. The data are derived from Langmuir probe measurements. The electron temperature at the extraction plane is still quite high in this source and rises across the filter plane. However the model cannot follow the rise observed experimentally and gives a result 1 eV less than experimentally observed in the driver region possibly due to the difficulty of including the physical presence of the rods within the model.

Zoom In Zoom Out Reset image size

It is shown in the second paper that the main contributing factor to deuteron fraction is the large rise in current density in the driver region (about 4 fold) which increases molecular gas dissociation by electron collisions and hence direct deuteron production by atomic ionization. Molecular ion dissociation (which was thought to be the main reason for the increase [7]) is still present but is a lesser contributor

The PINI source used on the JET tokamak as a neutral beam injector for heating the toroidal plasma has interior dimensions 0.55 × 0.31 × 0.21 m³. The usual discharge gas is deuterium and there are extensive measurements of the beam fractions using a Doppler shift diagnostic which measures the intensity of full energy, half energy and one third energy D⁰ atom velocities. The latter two particles originate from break-up of the ${\rm D}_{2}^{+}$ and ${\rm D}_{3}^{+}$ ions after ion-gas molecule collisions in the neutralizer of the injector. The intensities of these two D⁰ fractions can be traced back to the original ion species in the discharge [26]. In all of the PINI data the plasma grid is floating, so the grid current is zero.

Without a filter field the maximum deuteron fraction approaches 75% with reasonable quantitative agreement of the model as can be seen in figure 8. With hydrogen, see paper II, the proton fraction is lower, reaching a maximum value of only 60%.

Zoom In Zoom Out Reset image size

The supercusp is formed by arranging a net north pole on the backplane of the source (behind the filaments) and a ring of net south poles near the extraction plane so that the field has a form like a 'circus tent' with the filaments on the outside of the 'tent' (see figure 4). This minimizes the uniformity problems encountered by dipole filter used by Pincosy et al [7] arising from J × B drift orthogonal to the filter field and the source central axis. It was a successful approach [5, 6] as the non-uniformity was less than 5% of the plasma current density over the entire aperture array of 450 × 200 mm². Despite the filter field being at an angle of about 40⁰ to the source axis, the model works well in this situation, providing the field is reduced by the cosine of this angle, as is shown in figure 9. In this instance, the deuteron fraction reaches 91%.

Zoom In Zoom Out Reset image size

Another test of the model that can be applied to the PINI source is calculation of the discharge efficiency; i.e. the current density available at the extraction plane for a given arc current as measured at the arc power supply (this is larger than the emission current used in the model due to the back ion bombardment of the cathode). This is shown in figure 10 below for a deuterium discharge with no filter field and with a supercusp filter:

Zoom In Zoom Out Reset image size

The efficiency is noticeably lower for the supercusp source than the filterless source, a result that is also noted for the 10 × 10 source [7]. The main cause of the lower efficiency is the result of the lower plasma transport across the magnetic filter. As almost all of the ionization by fast primary electrons occurs in the driver region, near the filaments, the extraction plane becomes isolated from the ionization region and the plasma density and current density falls. The driver region, in contrast, shows a rapid rise in plasma density as the loss areas in the driver region (the filament structure and anode cusp lines) are small relative to the plasma grid area. This is seen in figure 7 for the 10 × 10 source, where the driver current density is much higher than at the extraction plane. The cusp losses, which are normally trivial, become significant at higher filter fields. The arc efficiency of very strong filters is even lower and will be discussed in paper II.

Several variants of the Jet PINI source have been developed including a version with a much greater depth which tests the effects of increasing the distance over which the molecular ions can be broken up both with and without a filter field and a version used to create negative ions [8] which has a much stronger supercusp filter (see paper II).

The long PINI, i.e. the source of greater depth, which is 270 mm deep compared with the standard design of 200 mm depth, has the same filter field in the same position relative to the backplane of the source, so that the filter field is now further from the extraction plane. This does increase the deuteron yield as shown in figure 11.

Zoom In Zoom Out Reset image size

However the increase in deuteron yield from 90% to 93% is at the expense of a decrease in discharge efficiency from 1.55 A m⁻² A⁻¹ of arc current to about 0.95 A m⁻² A⁻¹ of arc current by comparison with figure 9. For this reason, the long PINI was not considered a success and a similar result was reported by Pincosy [7].