# A 1.5D fluid-Monte Carlo model of a hydrogen helicon plasma

R Agnello, G Fubiani, I Furno, Ph Guittienne, A Howling, R Jacquier, F Taccogna

► To cite this version:

R Agnello, G Fubiani, I Furno, Ph Guittienne, A Howling, et al. A 1.5D fluid-Monte Carlo model of a hydrogen helicon plasma. Plasma Physics and Controlled Fusion, 2022, 64 (5), pp.055012. 10.1088/1361-6587/ac5ca2. hal-04828580

HAL Id: hal-04828580

https://ut3-toulouseinp.hal.science/hal-04828580v1

Submitted on 10 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# A 1.5D fluid—Monte Carlo model of a hydrogen helicon plasma

To cite this article: R Agnello et al 2022 *Plasma Phys. Control. Fusion* **64** 055012

View the [article online](#) for updates and enhancements.

## You may also like

- Diagnosis of helicon plasma by local OES
Ruilin Cui, Ruoyu Han, Kaiyi Yang et al.

- Experimental characterization and modelling of the resistive wall mode response in a reversed field pinch
E A Saad and P R Brunsell

- Interpreting radial correlation Doppler reflectometry using gyrokinetic simulations
J Ruiz Ruiz, F I Parra, V H Hall-Chen et al.
# A 1.5D fluid—Monte Carlo model of a hydrogen helicon plasma

R Agnello<fnref n="1" />,*, G Fubiani<fnref n="2" /> ID, I Furno<fnref n="1" /> ID, Ph Guittienne<fnref n="1" />, A Howling<fnref n="1" /> ID, R Jacquier<fnref n="1" /> and F Taccogna<fnref n="3" /> ID

E-mail: riccardo.agnello@epfl.ch

Received 5 October 2021, revised 23 February 2022

Accepted for publication 10 March 2022

Published 11 April 2022

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

Helicon plasma sources operating with hydrogen or deuterium might be attractive for fusion applications due to their higher power efficiency compared to inductive radiofrequency plasma sources. In recent years, the resonant antenna ion device (RAID) has been investigating the physics of helicon plasmas and the possibility of employing them to produce negative ions for heating neutral beam injectors (HNBs). Herein, we present a fluid Monte Carlo (MC) model that describes plasma species transport in a typical helicon hydrogen plasma discharge. This work is motivated by an interest in better understanding the basic physics of helicon plasma devices operating in hydrogen and, in particular, the volume production of negative ions. This model is based on the synergy between two separate self-consistent approaches: a plasma fluid model that calculates ion transport and an MC model that determines the neutral and rovibrational density profiles of H₂. By introducing the electron density and the temperature profiles measured by Langmuir probes as model constraints, the densities of ion species (H⁺, H₂⁺, H₃⁺, H⁻) are computed in a 1.5D (dimensional) geometry. The estimate of the negative ion density profile represents a useful benchmark that is comparable with dedicated diagnostics, such as cavity ring-down spectroscopy and Langmuir probe laser photodetachment. Neutral gas particles (atoms and molecules) are calculated assuming a fixed plasma background. This gas-plasma decoupling is necessary due to the different timescales of plasma (microseconds) and gas kinetics (milliseconds).

Keywords: high-power helicon plasmas, volume production of negative ions, fluid plasma model, Monte Carlo model

(Some figures may appear in colour only in the online journal)

<figure><img src="image_3.png" /></figure>

Original Content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ École Polytechnique Fédérale de Lausanne (EPFL), Swiss Plasma Center (SPC), CH-1015 Lausanne, Switzerland</footnote>

<footnote n="2">² LAPLACE, Université de Toulouse, CNRS, INPR, UPS, Toulouse, France</footnote>

<footnote n="3">³ Istituto di Metodologie Inorganiche e dei Plasmi-CNR, Via Amendola 122/D, 70126 Bari, Italy</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
# 1. Introduction

The fusion reactor ITER will be equipped with two heating neutral beam (HNB) injectors to heat the plasma to the required conditions for controlled nuclear fusion [1]. To produce high-energy (up to 1 MeV) atomic beams, negative ions ($H^-$ or $D^-$) are accelerated and then neutralized by charge exchange with a background gas along the beam propagation direction. Negative hydrogen and deuterium ions for HNB applications are routinely produced in large plasma sources, such as SPIDER [2] and ELISE [3], by means of radiofrequency inductive plasma sources. Our understanding of the physics mechanisms describing the transport and extraction of negative hydrogen and deuterium ions in these sources is supported by numerical simulations [4, 5]. Low-pressure hydrogen and deuterium plasmas are composed of different ionic and gas species, including positive and negative ions, ground-state atoms and molecules, atomic and molecular electronic excited states, molecular rovibrational states, and radicals.

In recent years, the possibility of producing negative ions ($H^-$ or $D^-$) for HNB injectors by means of helicon plasma sources has been investigated at the resonant antenna ion device (RAID) experiment [6, 7]. In parallel, other helicon plasma sources operating worldwide have studied the volume production of negative ions in helicon-sustained plasmas [8–11]. In RAID, plasmas are sustained by a resonant birdcage antenna used as a helicon plasma source [12] for up to 10 kW long-pulse operation. Numerical simulations are under way to investigate the physics of helicon waves [13] and plasma transport [14] in RAID for typical H/D plasma discharges. As part of our numerical simulation efforts and cross-check versus experimental data, we present herein a 1.5D fluid–Monte Carlo (MC) model aimed at the description of the transport of plasma and neutral species in a typical hydrogen helicon plasma produced by RAID. As a ‘1.5D’ model, we intend that it should describe the dynamics along the radial direction while retaining some axial physics. In particular, ion losses along the direction of the magnetic field are considered. The 1.5D geometry for RAID is justified by the cylindrical symmetry in the density and temperature of helicon plasmas and their uniformity along the plasma column’s axis.

Fluid models can be employed to understand the physical principles of certain magnetized plasmas; in particular, they can be very successful in low pressure, low beta plasmas, which are relatively more stable. Advanced 2D and 3D models, required by the lack of particular symmetries, are currently used to describe plasmas produced by inductive plasma sources for HNB injectors [15]. Fluid models have been employed to describe highly collisional plasmas, such as those on the edge of fusion plasmas, small toroidal experiments (TORPEX), and steady-state linear plasma devices [16]. Unlike other helicon linear reactors in which the plasma is bounded by long and narrow dielectric tubes, the RAID vessel walls are far from the edge of the plasma column and are electrically conductive. This suggests that the boundary

conditions are mainly dominated by the endplates, on which important plasma losses and recombination processes occur. A comprehensive description of a hydrogen helicon discharge would require coupling the helicon wave propagation and power deposition mechanisms to a transport/chemistry model. This is beyond the scope of the present work, which aims to investigate the roles of plasma reactions and species profiles, which cannot easily be separately measured.

The paper is structured as follows: in section 2, the RAID device and its parameters are presented. In section 3, the equations of the fluid model are discussed, as well as the plasma chemistry and the model inputs. In section 4, the MC model for neutral particles is described. In section 5, the output data of the model are presented, including plasma species profiles, source profiles, generation/loss mechanisms for $H^-$, and a comparison with experimental measurements. In section 6, we discuss the conclusions and the possible applications of this model.

# 2. Experimental setup and parameters

Plasma discharges in the RAID device are excited by a radiofrequency (RF) antenna in a birdcage geometry that is able to sustain the propagation of helicon waves [12]. The antenna is wrapped around a ceramic cylinder with an external diameter of 11.5 cm, an internal diameter of 9.5 cm, and a length of 38 cm mounted on one end of a stainless steel cylindrical vessel 1.5 m long and 40 cm in diameter; it can deliver up to 10 kW to the plasma. A set of six coils generates a DC magnetic field of up to 800 mT along the axis. Further details of RAID can be found in other works [6, 7, 17].

The fluid model can be applied provided that the ion-neutral collision mean free path (MFP) is smaller than the typical size of the device. Table 1 summarizes typical RAID hydrogen plasma and device parameters. As in the case of hot fusion plasmas, the magnetic field can play an important role in plasma transport in low-pressure plasma sources [5]. The typical values of the RAID magnetic field (0.01–0.1 T) are such that the ion and electron Larmor radii are in the submillimeter range, therefore the plasma is strongly magnetized and charged particles are constrained to flow along the magnetic field lines. Cross-field diffusion is, however, possible due to particle collisions that produce displacements along the radial direction, whose size is of the order of the Larmor radius. Along the axial direction, charged particles freely stream and undergo collisions with the background gas ($H_2$) with an MFP of the order of 0.1 m, shorter than the plasma axial length (1.8 m). The degree of anisotropy in the ion transport is such that the typical ratio between the axial and perpendicular diffusion coefficients is [18]:

$$D_{\parallel}/D_{\perp} = 1 + \left(\frac{\omega_c}{\nu}\right)^2 \approx 10^2 - 10^3, \quad (1)$$

where $\omega_c$ is the ion cyclotron frequency and $\nu$ the typical ion–gas collision frequency.
<table><caption>**Table 1.** Typical parameters for a hydrogen plasma discharge in RAID. The size of the ion Larmor radius, compared to the size of the vessel, makes it possible to describe the ion transport using a fluid model.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Vessel radius</td><td>0.2 m</td></tr><tr><td>Gas pressure</td><td>0.15–1.5 Pa</td></tr><tr><td>Magnetic field</td><td>0.01–0.1 T</td></tr><tr><td>Plasma axial length</td><td>1.8 m</td></tr><tr><td>Plasma density</td><td>≤2 × 10<sup>18</sup> m<sup>-3</sup></td></tr><tr><td>Electron temperature</td><td>≤8 eV</td></tr><tr><td>Ion temperature</td><td>≤0.1 eV</td></tr><tr><td>Ion cyclotron frequency</td><td>3.05 × 10<sup>5</sup> s<sup>-1</sup> (B = 20 mT)</td></tr><tr><td>Electron cyclotron frequency</td><td>5.6 × 10<sup>8</sup> s<sup>-1</sup> (B = 20 mT)</td></tr><tr><td>H<sub>2</sub><sup>+</sup>–H<sub>2</sub> collision frequency</td><td>2.36 × 10<sup>4</sup> s<sup>-1</sup></td></tr><tr><td>Electron–H<sub>2</sub> collision frequency</td><td>1.19 × 10<sup>5</sup> s<sup>-1</sup></td></tr><tr><td>Electron–H<sub>2</sub><sup>+</sup> collision frequency</td><td>2.5 × 10<sup>6</sup> s<sup>-1</sup></td></tr><tr><td>Electron–electron collision frequency</td><td>6.25 × 10<sup>6</sup> s<sup>-1</sup></td></tr><tr><td>H<sub>2</sub><sup>+</sup>–H<sub>2</sub><sup>+</sup> collision frequency</td><td>2.61 × 10<sup>7</sup> s<sup>-1</sup></td></tr><tr><td>H<sub>2</sub><sup>+</sup>–H<sub>2</sub> mean free path</td><td>≈0.1 m</td></tr><tr><td>electron–H<sub>2</sub><sup>+</sup> mean free path</td><td>≈5 mm</td></tr><tr><td>Ion Larmor radius</td><td>≤1 mm</td></tr><tr><td>Electron Larmor radius</td><td>≤0.25 mm</td></tr><tr><td>Debye length</td><td>≈10<sup>-5</sup> m</td></tr></tbody></table>

The RAID vessel is supplied at an H₂ flow rate of 10 sccm and a constant pressure of 0.3 Pa, corresponding to a molecular density of 7.25 × 10¹⁹ m⁻³ at room temperature (300 K). Given the typical plasma density generated in the center of a hydrogen plasma column (2.5 × 10¹⁸ m⁻³), we make the assumption that plasma species are diluted in an H₂ background. Ion transport can then be described by classical drift-diffusion equations in the presence of a uniform magnetic field. At this low RF power (≤5 kW) and gas pressure (p ≤ 1 Pa), neutral depletion [19], consisting of the displacement of neutrals to the walls of the device at high plasma densities (≥10¹⁹ m⁻³), can be neglected. As far as Coulomb collisions are concerned, in the case of like particles (electron–electron and ion–ion collisions), these do not cause diffusion [20]. Only electron–ion collisions might play a role. The electron–ion collision frequency is about one order of magnitude larger than the electron–neutral collision frequency; however, due to the large mass difference, the relative momentum exchange is of the order of Δp/p ≈ √(m_e/m_i)√(T_e/T_i) ≈ 0.15; therefore, Coulomb collisions should not greatly affect ion diffusion with respect to the background neutral gas. Plasma instability and turbulence might occur in a mid/high power helicon plasma device, such as RAID, which would require us to take account of nonlinear terms, dramatically complicating the model. However, as we will see later on, RAID plasmas are well described by classical diffusion in a magnetic field, so we do not take nonlinear phenomena into account in our first-order modeling.

# 3. Description of the model

## 3.1. Basic equations

The equations composing the fluid description are the continuity, momentum, and Poisson equations, which are given below. Each ion species is denoted by the subscript *s*, where *s* can be H⁺, H₂⁺, H₃⁺, or H⁻. The continuity and the ion momentum conservation equations are:

$$ \frac{\partial n_s}{\partial t} + \nabla \cdot \vec{\Gamma}_s = S_s - L_s, \quad (2) $$

$$ m_s n_s \left[ \frac{\partial \vec{u}_s}{\partial t} + (\vec{u}_s \cdot \nabla) \vec{u}_s \right] = q_s n_s (\vec{E} + \vec{u}_s \times \vec{B}) - \nabla \cdot \vec{P}_s \\ - m_s n_s \nu_s \vec{u}_s, \quad (3) $$

where $n_s$ is the number density, $m_s$ is the mass, $\vec{u}_s$ is the fluid velocity of the species, $\vec{E}$ is the electric field, $\vec{B}$ is the magnetic field, $\vec{P}_s$ is the partial pressure tensor (due to ions only), $q_s$ is the electron charge magnitude, $\vec{\Gamma}_s = n_s \vec{u}_s$ is the flux of species, $\nu_s$ the momentum transfer frequency, $S_s$ and $L_s$ are the source and loss terms.

The electric field and the plasma potential *V* are calculated using the Poisson equation:

$$ \nabla^2 V = -\frac{e}{\epsilon_0} (n_{H^+} + n_{H_2^+} + n_{H_3^+} - n_{H^-} - n_e), \quad (4) $$

with

$$ \vec{E} = -\nabla V, \quad (5) $$

where $n_e$ is the electron density and $\epsilon_0$ is the vacuum dielectric constant. Plasma quasi-neutrality is guaranteed by the electric field generated by the charge separation. Due to the cylindrical symmetry of the RAID plasma, the magnetic DC field as well as the electric field generated by the plasma potential have no azimuthal components. However, they produce an $\vec{E} \times \vec{B}$ and a diamagnetic drift along the azimuthal direction of the charged fluids with associated tangential fluxes. If we consider axial invariance, the total flux can be expressed as the sum of a radial and an azimuthal component, namely $\vec{\Gamma} = \vec{\Gamma}_r + \vec{\Gamma}_\phi$. Equations (2)–(4) are then solved with the PDE (partial differential equation) module of the finite element calculator COMSOL multiphysics [21].

Assuming cylindrical symmetry of the plasma discharge, the continuity equation can be written as:

$$ \frac{1}{r} \frac{\partial(r\Gamma_r)}{\partial r} = S(r) - 2 \frac{\Gamma_z(r)}{L_z}, \quad (6) $$

where $S(r)$ is the source term of any species along the radial direction (it takes into account the sum of sources and losses in volume) and $2\frac{\Gamma_z(r)}{L_z}$ is the term describing the particle losses along the axial direction *z* due to the free streaming along the magnetic field lines up to both endplates. Due to this term, we define the model to have a 1.5D geometry. $\Gamma_z(r)$ is the flux of
particles at the discharge ends at radius r. For positive ions, $\Gamma_z(r)$ reads:

$$ \Gamma_z(r) = \alpha n_s u_{\text{Bohm}}(r), \quad (7) $$

where $u_{\text{Bohm}}$ is the Bohm velocity $u_{\text{Bohm},s} = \sqrt{q_s T_e / m_s}$. Here, $\alpha$ is a factor that depends on the presheath voltage drop and the conditions in the collisional presheath, in particular, on the plasma flow velocity and the electron and ion temperatures. The presheath drop dictates the density at the sheath.

From equation (3) one can derive the expression for the radial flux; for positive ions, this reads:

$$ \Gamma_r = \frac{q_s \nu_s n_s}{m_s (\nu_s^2 + \omega_{c,s}^2)} \left( E_r - \frac{T_s}{n_s} \frac{dn_s}{dr} \right), \quad (8) $$

where $\omega_{c,s} = q_s B / m_s$ is the ion cyclotron frequency.

## 3.2. Assumptions

It is useful to summarize the main simplifying assumptions in this paragraph.

We neglect the acceleration term, $\partial \vec{u}_s / \partial t$, since we assume we are in the steady-state regime. For the ion inertia term, $(\vec{u}_s \cdot \nabla) \vec{u}_s$, we assume that diffusive timescales are much longer than bulk flow timescales [18]. The inclusion of the inertial term would lead to the separation of the momentum equation in each dimension for each species, increasing the complexity of the model [22].

The ion distribution function is considered isotropic and, also, we neglect the gradient of the ion temperature (we assume a 0.1 eV uniform ion temperature and the same for all ion species); therefore, $\nabla \cdot \vec{P}_s = k_b T_s \nabla n_s$, where $k_b$ is Boltzmann's constant.

The magnetic field $\vec{B}$ only takes into account the DC field generated by the external coils (200 mT), since the magnetic field generated by the plasma currents is negligible (the magnetic field of the helicon wave is a few Gauss [13]).

The magnetic pressure is significantly larger than the plasma pressure, rendering all plasma diamagnetic effects negligible. The ratio between the plasma kinetic pressure and the plasma magnetic pressure is given by:

$$ \frac{p_{\text{magnetic}}}{p_{\text{plasma}}} = \frac{B^2 / 2\mu_0}{k_b n_e T_e} \approx 100. \quad (9) $$

## 3.3. Boundary conditions

The vessel walls and endplates play an important role in the plasma equilibrium. They govern the rate at which particles are lost from the plasma bulk as well as the surface chemistry. This is particularly relevant in the case of a linear plasma device, such as RAID, in which magnetic field lines intersect the boundaries of the device. In the literature, some authors have investigated the effects of changing the positions and the electric properties of the endplates [23], which have an effect

on the axial boundary conditions of the helicon wave and, in turn, on the discharge properties.

A rigorous description of wall losses would require a self-consistent treatment of the plasma sheath. In the sheath region (a few mm in front of the wall), quasi-neutrality is lost and therefore strong electric fields appear. From the modeling point of view, this would imply strong charged-particle acceleration and problems of numerical convergence. To accurately resolve the dynamics of ion free fall across the plasma sheath, one should locally redefine the meshing with at least a submillimetric resolution. The coupling of the bulk plasma region with the plasma sheath region would require us to considerably increase the complexity of the transport model without adding any particularly relevant element to the overall transport dynamics and chemistry. Therefore, the ion losses at the boundaries are estimated according to the Bohm expression, which reads as follows for positive ions:

$$ \Gamma_{\text{wall}} = -\sqrt{(qT_e/m_s)} n_s e^{-\frac{1}{2}}. \quad (10) $$

Negative ions (such as $H^-$), due to their negative charge, are confined in the plasma bulk in electropositive plasmas [24]. The typical temperature of ions in low-temperature plasma devices (~0.1 eV) is indeed insufficient to overcome the plasma potential well. $H^-$ ions are therefore produced and destroyed in the plasma volume. In a volume negative ion source, such as RAID, it is therefore crucial to understand the interplay between the production and destruction channels of negative ions (discussed in section 5.2), and how they can be controlled to optimize negative ion yield.

## 3.4. Reaction rates

In this model, we consider volume reactions between electrons and neutral/ion species, between atomic and ion species, and between ion species in the plasma bulk. We do not consider three-body reactions, since they are negligible in this low-pressure regime. The reaction rate R is given by the product of the densities $n_1$ and $n_2$ of the reactants and is regulated by a rate coefficient $k(T_1, T_2)$ that depends on the reactants' temperatures, $T_1$ and $T_2$:

$$ R = k(T_1, T_2) n_1 n_2 \quad (\text{m}^{-3} \text{s}^{-1}). \quad (11) $$

To describe the chemistry of the hydrogen plasma discharge in RAID, we employed the reaction set in table 2, which shows the complexity of the plasma chemistry in this model.

For reaction 3, we employed the dissociative attachment reaction rate provided in [25], whereas for the others, we employed the reaction rates provided in [22]. The channels for the volume production of $H^-$ are given by reactions 3, 7, and 22. We will see that reaction 3, the dissociative attachment, is by far the most effective channel for $H^-$ production in the plasma. The net source term for each ion species is reported in the set of equations (12):
<table><caption>**Table 2.** Plasma reactions and reaction rates used for this fluid transport model.</caption><thead><tr><th></th><th>Reaction</th><th>Reaction rate (k<sub>n</sub>) (m<sup>3</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1</td><td>e + H<sub>2</sub> → H<sub>2</sub><sup>+</sup> + 2e</td><td>7.27 × 10<sup>-15</sup> T<sub>e</sub><sup>0.549</sup> e<sup>-15.5/(T<sub>e</sub>-0.001)</sup> + 1.37 × 10<sup>-14</sup> T<sub>e</sub><sup>-0.557</sup> e<sup>-20.3/T<sub>e</sub></sup></td></tr><tr><td>2</td><td>e + H<sub>2</sub> → 2H + e</td><td>1.54 × 10<sup>-14</sup> T<sub>e</sub><sup>0.06</sup> e<sup>8.63/T<sub>e</sub></sup> + 1.11 × 10<sup>-13</sup> T<sub>e</sub><sup>-0.813</sup> e<sup>-13.4/T<sub>e</sub></sup></td></tr><tr><td>3</td><td>e + H<sub>2</sub>(ν = 0 - 14) → H + H<sup>-</sup></td><td>1.972σ<sub>ν</sub><sup>(0)</sup> <sup>T<sub>e</sub><sup>1/2</sup></sup>/<sub>1+T<sub>e</sub>/E<sub>0</sub></sub> e<sup>-<sup>ΔE<sub>ih,ν</sub></sup>/<sub>T<sub>e</sub></sub> [<sup>|E<sub>ih,ν</sub>|</sup>/<sub>T<sub>e</sub></sub> + <sup>1</sup>/<sub>1+T<sub>e</sub>/E<sub>0</sub></sub>] × (10<sup>-14</sup>)</td></tr><tr><td>4</td><td>e + H → H<sup>+</sup> + 2e</td><td>k<sub>4</sub> = 9.74 × 10<sup>-15</sup> T<sub>e</sub><sup>-0.174</sup> e<sup>-14.3/(T<sub>e</sub>-0.001)</sup> + 6.35 × 10<sup>-15</sup> T<sub>e</sub><sup>0.433</sup> e<sup>-16.4/T<sub>e</sub></sup></td></tr><tr><td>5</td><td>H + H<sub>2</sub><sup>+</sup> → H<sub>2</sub> + H<sup>+</sup></td><td>k<sub>4</sub> = 1.54 × 10<sup>-14</sup> T<sub>e</sub><sup>0.859</sup> e<sup>-4.61/(T<sub>1</sub>-0.786)</sup> + 1.64 × 10<sup>-15</sup> T<sub>i</sub><sup>-0.353</sup> e<sup>-0.258/T<sub>i</sub></sup></td></tr><tr><td>6</td><td>e + H<sub>2</sub><sup>+</sup> → H + H<sup>+</sup> + e</td><td>3.50 × 10<sup>-13</sup> T<sub>e</sub><sup>-1.25</sup> e<sup>-3.19/(T<sub>e</sub>-0.008)</sup> + 1.77 × 10<sup>-13</sup> T<sub>e</sub><sup>-0.0924</sup> e<sup>-3.04/T<sub>e</sub></sup></td></tr><tr><td>7</td><td>e + H<sub>3</sub><sup>+</sup> → H<sub>2</sub><sup>+</sup> + H<sup>-</sup></td><td>1.93 × 10<sup>-15</sup> T<sub>e</sub><sup>-1.07</sup> e<sup>-6.26/(T<sub>e</sub>+0.131)</sup> + 5.35 × 10<sup>-17</sup> T<sub>e</sub><sup>-0.371</sup> e<sup>-2.07/T<sub>e</sub></sup></td></tr><tr><td>8</td><td>H<sub>2</sub> + H<sub>2</sub><sup>+</sup> → H + H<sub>3</sub><sup>+</sup></td><td>6.29 × 10<sup>-15</sup> T<sub>i</sub><sup>-1.46</sup> e<sup>-2.22/(T<sub>i</sub>+0.356)</sup> + 2.71 × 10<sup>-16</sup> T<sub>e</sub><sup>-1.30</sup> e<sup>-0.317/T<sub>i</sub></sup></td></tr><tr><td>9</td><td>H<sup>+</sup> + H<sup>-</sup> → 2H</td><td>4.46 × 10<sup>-14</sup> T<sub>i</sub><sup>-0.281</sup> + 1.26 × 10<sup>-14</sup> e<sup>-1.96/T<sub>i</sub></sup></td></tr><tr><td>10</td><td>H<sub>2</sub><sup>+</sup> + H<sup>-</sup> → H<sub>2</sub> + H</td><td>2.23 × 10<sup>-14</sup> T<sub>i</sub><sup>0.425</sup> + 8.93 × 10<sup>-14</sup> T<sub>i</sub><sup>-0.261</sup></td></tr><tr><td>11</td><td>H<sub>3</sub><sup>+</sup> + H<sup>-</sup> → 0.5H<sub>2</sub> + 3H</td><td>1.70 × 10<sup>-14</sup> T<sub>i</sub><sup>0.313</sup> + 5.75 × 10<sup>-14</sup> T<sub>i</sub><sup>-0.288</sup></td></tr><tr><td>12</td><td>e + H<sup>-</sup> → H + 2e</td><td>4.58 × 10<sup>-13</sup> T<sub>e</sub><sup>0.287</sup> e<sup>-4.41/(T<sub>e</sub>+0.117)</sup> + 2.71 × 10<sup>-14</sup> T<sub>e</sub><sup>0.62</sup> e<sup>-1.82/T<sub>e</sub></sup></td></tr><tr><td>13</td><td>e + H<sub>3</sub><sup>+</sup> → 2H + H<sup>+</sup> + e</td><td>2.69 × 10<sup>-13</sup> T<sub>e</sub><sup>-0.245</sup> e<sup>-15.6/(T<sub>e</sub>+0.003)</sup> + 1.01 × 10<sup>-12</sup> T<sub>e</sub><sup>-0.464</sup> e<sup>-26.8/T<sub>e</sub></sup></td></tr><tr><td>14</td><td>H<sub>2</sub> + H<sup>+</sup> → H<sub>3</sub><sup>+</sup> + hν</td><td>1.63 × 10<sup>21</sup></td></tr><tr><td>15</td><td>e + H<sub>2</sub><sup>+</sup> → 2H</td><td>2.29 × 10<sup>-14</sup> T<sub>e</sub><sup>-0.571</sup> + 3.31 × 10<sup>-15</sup> T<sub>e</sub><sup>-0.152T<sub>e</sub></sup></td></tr><tr><td>16</td><td>e + H<sub>3</sub><sup>+</sup> → 3H</td><td>3.36 × 10<sup>-15</sup> T<sub>e</sub><sup>-0.716</sup> + 3.73 × 10<sup>-14</sup> T<sub>e</sub><sup>-0.67</sup> e<sup>6.40/T<sub>e</sub></sup></td></tr><tr><td>17</td><td>e + H<sub>3</sub><sup>+</sup> → H + H<sub>2</sub></td><td>2.03 × 10<sup>-15</sup> (T<sub>e</sub><sup>-0.189</sup> + 0.040T<sub>e</sub><sup>-1.49</sup>) + 5.57 × 10<sup>-14</sup> T<sub>e</sub><sup>1.23</sup> e<sup>-6.21/T<sub>e</sub></sup></td></tr><tr><td>18</td><td>H + H<sup>-</sup> → H<sub>2</sub> + e</td><td>2.16 × 10<sup>-13</sup> T<sub>i</sub><sup>-1.89</sup> e<sup>-12.7/(T<sub>i</sub>+1.17)</sup> + 1.30 × 10<sup>-15</sup> T<sub>i</sub><sup>-0.418</sup> e<sup>-0.192/T<sub>i</sub></sup></td></tr><tr><td>19</td><td>H<sub>2</sub> + H<sup>-</sup> → H + H<sub>2</sub> + e</td><td>1.62 × 10<sup>-16</sup> T<sub>i</sub><sup>0.417</sup> e<sup>-6.47/(T<sub>i</sub>+0.132)</sup> + 5.70 × 10<sup>-16</sup> T<sub>i</sub><sup>0.550</sup> e<sup>-2.19/T<sub>i</sub></sup></td></tr><tr><td>20</td><td>H + H<sup>-</sup> → 2H + e</td><td>3.81 × 10<sup>-15</sup> T<sub>i</sub><sup>0.280</sup> e<sup>-3.76/(T<sub>i</sub>+0.626)</sup> + 4.55 × 10<sup>-16</sup> T<sub>i</sub><sup>0.603</sup> e<sup>-0.375/T<sub>i</sub></sup></td></tr><tr><td>21</td><td>H<sub>2</sub> + H<sup>+</sup> → H + H<sub>2</sub><sup>+</sup></td><td>5.54 × 10<sup>-16</sup> T<sub>i</sub><sup>-0.453</sup> e<sup>-3.26/(T<sub>i</sub>-0.001)</sup> + 5.98 × 10<sup>-18</sup> T<sub>i</sub><sup>-2.88</sup> e<sup>-0.310/T<sub>i</sub></sup></td></tr><tr><td>22</td><td>e + H → H<sup>-</sup> + hν</td><td>5.75 × 10<sup>-20</sup> (T<sub>e</sub><sup>-0.0285</sup> - 0.94T<sub>e</sub><sup>-0.05</sup>) + 6.54 × 10<sup>-19</sup> T<sub>e</sub><sup>-5.18</sup> e<sup>-72.4/T<sub>e</sub></sup></td></tr></tbody></table>

$$
\left\{
\begin{aligned}
\nabla \cdot \Gamma_{H^+} &= k_4 n_e n_H + k_5 n_H n_{H^+} + k_6 n_e n_{H_2^+} - k_9 n_{H^+} n_{H^-} + k_{13} n_e n_{H_3^+} - k_{14} n_{H_2} n_H - k_{21} n_{H_2} n_{H^+}, \\
\nabla \cdot \Gamma_{H_2^+} &= k_{11} n_e n_{H_2} - k_5 n_H n_{H_2^+} - k_6 n_e n_{H_2^+} + k_7 n_e n_{H_3^+} - k_8 n_{H_2} n_{H_2^+} - k_{10} n_{H_2^+} n_{H^-} - k_{15} n_e n_{H_2^+} + k_{21} n_{H_2} n_{H^+}, \\
\nabla \cdot \Gamma_{H_3^+} &= -k_7 n_e n_{H_3^+} + k_8 n_{H_2} n_{H_2^+} - k_{11} n_{H_3^+} n_{H^-} - k_{13} n_e n_{H_3^+} + k_{14} n_{H_2} n_{H^+} - k_{16} n_e n_{H_3^+} - k_{17} n_e n_{H_3^+}, \\
\nabla \cdot \Gamma_{H^-} &= n_e \sum_{\nu=0}^{14} k_3(\nu) n_{H_2(\nu)} + k_7 n_e n_{H_3^+} - k_9 n_{H^+} n_{H^-} - k_{10} n_{H_2^+} n_{H^-} - k_{11} n_{H_3^+} n_{H^-} - k_{12} n_e n_{H^-} + \\
&\quad -k_{18} n_H n_{H^-} - k_{19} n_{H_2} n_{H^-} - k_{20} n_H n_{H^-} + k_{22} n_e n_H.
\end{aligned}
\right.
\quad (12)
$$

The generation rate of H⁻ by dissociative attachment from ro-vibrationally excited H₂, which is the term $\sum_{\nu=0}^{14} k_3(\nu) n_{H_2(\nu)}$, is calculated by taking into account the density of H₂ molecules in each vibrational state $\nu$ computed by the MC model (see section 4). The most important H⁻ destruction processes are mutual neutralization (MN) with positive ions (R9–11), electron detachment (ED) (R12) and associative detachment (AD) with neutrals (R18–20).

### 3.5. Structure of the model

A complete fluid model describing all species would require the inclusion of a balance equation for electrons and an equation for the power deposition. In RAID, the electron temperature and density profiles can be measured by Langmuir probes and microwave interferometry across the entire volume of the device [4]. These data are introduced as a physical constraint on the basis of which the transport and chemistry of the
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Electron density and temperature radial profiles used as inputs in the model, obtained by fitting experimental data [26]. These profiles are representative for a hydrogen plasma discharge with 3 kW RF power and 0.3 Pa gas pressure at 1 m from the center of the antenna.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Flow chart describing the coupling between the fluid and the Monte Carlo model.</figcaption></figure>

other ion species and atomic hydrogen can be computed. Typical electron density and temperature profiles at 1 m from the center of the antenna, 3 kW RF power, and 0.3 Pa of gas pressure are shown in figure 1; these are taken as the input profiles for the modeling. Figure 2 provides a diagrammatic description of the structure of the model.

### 3.6. Dissociative attachment reactions

In plasmas, molecular hydrogen can be ro-vibrationally excited by electron impact or by decay from the excited molecule, sustaining 15 vibrational states in the electronic ground state. Experimental and theoretical studies have revealed that dissociative attachment to these vibrationally excited states is the main mechanism responsible for volume production of negative ions [27]. A detailed description of the overall hydrogen plasma transport would require the inclusion of reactions involving the collisions between H₂ in the various vibrational states and other plasma species and the respective

transport equations. This would increase the number of possible reactions to a few hundreds, increasing the complexity of the model without substantially changing the transport of ion species. Indeed, since for the RAID plasma regime we expect a small amount of excited H₂ (of the order of 1% compared to background H₂ [4]), the overall transport and chemistry of ions are expected to be weakly influenced by the presence of excited H₂. To approximate the contribution of the dissociative attachment to H⁻ production, the MAR (molecular activated recombination) rate coefficient could be used as a reasonable approximation of the dissociative attachment reaction rate $k_{DA}$ [28]. In this study, however, we use an MC model to estimate the radial density profiles of each ro-vibrational state of H₂ which are then inserted as inputs into the fluid part of the model.

## 4. A MC model for determining neutral density profiles

In order to increase the self-consistency of the model, the plasma dynamics needs to be coupled with the neutral kinetics. The H₂ vibrational kinetics is particularly important, since the negative ion production is heavily dependent on the molecular vibrational distribution function (VDF). For this purpose, an MC approach was used to describe the dynamics and kinetics of neutral gas particles (atoms and molecules) considering the plasma (electron and ion species) to be a fixed background with non-uniform density and temperature distributions (the inputs of the plasma fluid model). This gas-plasma decoupling is a necessary hypothesis due to the very different plasma and gas kinetic timescales; the first is on the scale of microseconds, while the latter is in the millisecond range.

The MC model [29] is one-dimensional in the radial direction and keeps the same geometrical reduction arguments as those used for the plasma fluid model, namely, the uniformity along the axial and azimuthal directions, as mentioned in the introductory paragraph. The simulation starts with a uniform radial distribution of molecules (corresponding to a gas pressure $P_{gas} = 0.3$ Pa and temperature $T_{H_2} = 300$ K); at every time step, a certain number of molecules are uniformly injected along the radial domain in order to keep the number of molecules fixed. This numerically mimics the continuous axial flow of molecules necessary to keep the gas pressure constant. Each particle contains information about its electronic (for atoms) and vibrational (for molecules) states. Only the fundamental electronic state of H₂ is taken into account. No information on the rotational state is considered, since the high ionization coefficient and the electron temperature involved guarantee a rotational equilibrium distribution. In the dissociative attachment reaction, the effect of the rotational excitation is smaller than that of the vibrational excitation [30]. The most relevant collisions with electrons and H⁺ ions (the dominant ion species, see section 5) are considered through the use of state-selective cross-sections [31–33], while neutral-neutral collisions are negligible due to the low-pressure regime. The equations of motion of the particles are solved, and the heterogeneous surface processes are taken into
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Radial averaged $n_{H_2(v)}$ VDF (normalized to the $v=0$ level density) result of the MC model, displaying a two-temperature distribution.</figcaption></figure>

account in the impacts with the wall by using appropriate energy-dependent and state-dependent coefficients. The cross-sections of neutral/ion-surface processes involve a complex kinetics, in which the cross-sections depend on neutral/ion and surface temperatures [34]. For atoms, the cross-sections depend on the reflection and recombination coefficients [35, 36], while for molecules, they depend on the relaxation and dissociation coefficients [37, 38]. After 4 ms, the gas subsystem reaches a steady state, exhibiting a global dissociation coefficient (equation (14)) of $D_d = 0.103$. The radial averaged VDF is reported in figure 3, which shows a typical parabolic Treanor distribution for the first four vibrational levels and a hyperbolic plateau for those $\ge 4$, corresponding to a two-temperature distribution $T_{01} = 2000$ K and $T_{07} = 7000$ K.

Once the $n_H$ and $n_{H_2}$ radial density profiles have been determined, the degrees of ionization and dissociation can be computed; these are shown in figure 4. The degree of ionization, defined as:

$$D_i = \frac{n_e}{n_H + n_{H_2}}, \qquad (13)$$

is peaked at the center where it attains about 2%. The degree of dissociation, defined as:

$$D_d = \frac{n_H}{n_H + 2n_{H_2}}, \qquad (14)$$

is slightly hollow at the center and around 10% across the vessel. The $D_d$ (dissociation degree) well at the center of the plasma column is mainly due to a combination of central plasma heating and H transport toward the walls. H atoms produced by dissociation take about 2.2 eV each; therefore, they are much hotter than the ~0.1 eV $H_2$ translational temperature.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Ionization and dissociation degrees derived from the MC simulation.</figcaption></figure>

# 5. Transport of ion species

## 5.1. Ion density equilibrium and potential profiles

By solving the system of coupled equations (2)–(4), the density of each ion species can be computed, as depicted in figure 5, together with the net source rate profiles (production—destruction) for $H^+, H_2^+, H_3^+$ and $H^-$. $H^+, H_2^+$, which show centrally peaked density profiles, while $H_3^+$ peaks at the edge. The $H^-$ density profile also shows a shell-like structure, whose width and absolute value are comparable to measurements performed with laser-based diagnostics [7]. $H^-$ net production takes place in a region with a width of $\approx 1$ cm on the edge of the plasma column at $\approx 5$ cm, and the $H^-$ density peaks at $\approx 5$ cm. This result suggests that transport plays an important role in the $H^-$ density equilibrium profiles and that a non-negligible quantity is transported to the center of the plasma column. At the position of the $H^-$ density peak, the $H_2^+$ and $H_3^+$ densities are also comparable to that of $H^+$. These results show that the plasma column is characterized by a hot and dense core region of ionizing plasma dominated by $H^+$ surrounded by a halo of $H_2^+$ and $H_3^+$ recombining plasma. In the latter region, $H^-$ is efficiently produced and can survive.

In figure 6 we show the computed profiles of the plasma potential, which is self-consistently obtained by the resolution of Poisson's equation with the constraint of potential grounding on the wall and the boundary targets. The presence of electric forces arising from charge displacement guarantees charge quasi-neutrality, so that an explicit charge conservation equation is not needed. Figure 6 shows the computed plasma potential profile. It shows a central region with a small peak and a monotonic decay to the wall. The absolute peak value is sensitive to axial losses ($\Gamma_z(r)$), which are not self-consistently calculated. The axial losses are finely tuned to match the computed potential profile with the experimental
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Net generation rate (on the left) and density profiles (on the right) of ion species computed by the fluid model, including experimental $n_{H^{-}}$ data.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Plasma potential calculated by the fluid model (dashed line) and experimental data (asterisks). Axial particle losses to the walls are finely tuned so that the computed $V_p$ matches the experimental points [26].</figcaption></figure>

$V_p$ profile [26]. However, the absolute profile of the plasma potential has only a weak influence on the computed plasma species profiles and transport. Since the plasma potential profile is strongly dependent on boundary conditions, a rigorous treatment would require us to include the physics of the sheath, which goes beyond the scope of this model.

To gain insight into radial plasma species transport across the vessel's section, it is instructive to plot the fluxes of different ions. Figure 7 shows the absolute values of the fluxes of each plasma species across the radius. The arrows show the directions of the fluxes. The $H^+$ and $H_2^+$ fluxes peak at about 2 cm from the center and are lost within a few centimeters due to volume destruction processes with the background gas. $H_3^+$ ions are transported from the shell region up to the vessel wall. $H^-$ ions are transported inward across the radius.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Fluxes of ions across the radius. Positive ions are transported outwards, while $H^-$ ions are transported inwards and destroyed by electron detachment and mutual neutralization with positive ions.</figcaption></figure>

## 5.2. Production and destruction mechanisms of $H^-$

Here, we study each generation and loss mechanism of $H^-$ to better understand the volume production and destruction processes. In conventional negative ion sources for HNB injectors, the production of negative ions is achieved in two steps that take place in two different volumes: a region called the *driver*, where molecular hydrogen (or deuterium) is dissociated, and an expanding region with a caesiated surface, where H/D atoms pick up electrons. In the case of a Cs-free negative ion source, such as RAID, where negative ions are only produced in the plasma volume, it is interesting to study the contributions of production and destruction processes along the radial direction.

In figure 8 we show the individual reaction rates for $H^-$ production (dashed lines) and destruction processes (solid lines) for electron collision reactions (left) and those for reactions involving only heavy species (right). These profiles are obtained by considering as input parameters the measured electron density and temperature profiles shown in figure 1. The dissociative attachment reaction ($e + H_2(\nu) \rightarrow H^- + H$) from $H_2$ ro-vibrational excited states is the main contributor to negative ion production. The main destruction processes in the plasma center are electron detachment ($e + H^- \rightarrow 2e + H$) and mutual neutralization with $H^+$ ($H^+ + H^- \rightarrow 2H$); however, at the edge of the plasma column, where the $H^-$ density is peaked ($r \approx 5$ cm), the most important destruction processes are mutual neutralization with $H_3^+$ and $H_2^+$ ($H^- + H_3^+ \rightarrow 0.5H_2 + 3H$, $H^- + H_2^+ \rightarrow H_2 + H$) and associative detachment ($H + H^- \rightarrow H_2 + e$). Thus, negative ions are mostly produced along the axis but are efficiently simultaneously
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Individual production (dotted lines) and destruction (continuous lines) processes of $H^−$ along the radius induced by electrons (left) and by ions and neutrals (right). The dissociative attachment of ro-vibrationally excited $H_2$ is the dominant production process. In the plasma center, $H^−$ is mainly detached by electron impact, while at the edge, mutual neutralization processes with positive ions and associative detachment are the main loss reactions.</figcaption></figure>

destroyed by electron detachment. They mostly survive on the edge of the column because $T_e$ is lower; at this location, net creation of $H^−$ occurs, as shown in figure 5, which gives the sum of all the production and destruction rates for $H^−$.

## 6. Discussion and conclusions

A 1.5D ion fluid model was developed to better understand ion transport and species reactions in a helicon hydrogen plasma. This method is based on the separation of a fluid treatment for ions and an MC model for the neutrals. On the basis of electron density and temperature measurements, the neutral equilibrium density profiles as well as the density of $H_2$ in each vibrational excited state can be found. These values are used as inputs for the 1.5D fluid transport model. Due to the interest in volume-negative ion production in RAID, it is helpful to investigate the interplay of the production and destruction processes of $H^−$. We compared the reaction rate profiles with the steady-state equilibrium densities to understand how ions are transported, and observed that negative ions have a net production rate in a shell at the edge of the plasma column, from where they diffuse toward the axis. We observe that $H^−$ ions are mainly destroyed in the plasma center by electron detachment and are mutually neutralized by $H_2^+$ and $H_3^+$ at the edge, close to the peak density position. Future possible developments in $H^−$ volume production in helicon sources should focus on finding the correct balance between these competitive mechanisms to maximize the volume density of $H^−$ close to the extraction region. Preliminary proof-of-principle negative ion extraction tests have recently been carried out using a radial geometry [39].

A first envisaged step of this model would be the inclusion of axial transport, which could be helpful to understand plasma composition in expanding plasmas, such as those produced in conventional negative ion sources for HNB injectors. An important improvement would consist in the inclusion of

a power balance equation and a power deposition profile of helicon waves to self-consistently calculate the electron density and temperature profiles. The inclusion of the full helicon wave equations and the self-consistent calculation of plasma formation would be a considerable leap, requiring a 3D geometry, since helicon waves have a full 3D structure [13].

## Data availability statement

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

## Acknowledgments

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work was supported in part by the Swiss National Science Foundation.

## ORCID iDs

R Agnello https://orcid.org/0000-0002-1462-6480

G Fubiani https://orcid.org/0000-0001-8920-9908

I Furno https://orcid.org/0000-0001-8348-1716

A Howling https://orcid.org/0000-0003-4417-4339

F Taccogna https://orcid.org/0000-0001-9361-9010

## References

[1] Hemsworth R S *et al* 2017 *New J. Phys.* **19** 025005

[2] Serianni G *et al* 2020 *Nucl. Fusion* **60** 066025

[3] Wünderlich D, Riedl R, Fröschle M, Fantz U and Heinemann B 2021 *Plasma* **2021** 172–82

[4] Fubiani G, Garrigues L, Hagelaar G, Kohen N and Boeuf J P 2017 *New J. Phys.* **19** 015002

[5] Hagelaar G J M and Oudini N 2011 *Plasma Phys. Control. Fusion* **53** 124032

[6] Furno I, Agnello R, Fantz U, Howling A, Jacquier R, Marini C, Plyushchev G, Guittienne P and Simonin A 2017 *EPJ Web Conf.* **157** 03014

[7] Agnello R *et al* 2020 *Nucl. Fusion* **60** 026007

[8] Pandey A *et al* 2017 *Rev. Sci. Instrum.* **88** 103509

[9] Pandey A *et al* 2019 *Plasma Phys. Control. Fusion* **61** 065003

[10] Santoso J *et al* 2015 *Phys. Plasmas* **22** 093513

[11] Santoso J *et al* 2018 *Plasma Sources Sci. Technol.* **27** 10LT03

[12] Guittienne P, Chevalier E and Hollenstein C 2005 *J. Appl. Phys.* **98** 083304

[13] Guittienne P *et al* 2021 *Plasma Sources Sci. Technol.* **30** 075023

[14] Fubiani G *et al* 2021 *Phys. Plasmas* **28** 063503

[15] Lishev S, Schiesko L, Wünderlich D, Wimmer C and Fantz U 2018 *Plasma Sources Sci. Technol.* **27** 125008

[16] Ricci P, Halpern F D, Jollie S, Loizu J, Mosetto A, Fasoli A, Furno I and Theiler C 2012 *Plasma Phys. Control. Fusion* **54** 124047

[17] Jacquier R *et al* 2019 *Fusion Eng. Des.* **146** 1140–4

[18] Chen F F 2006 *Plasma Physics and Controlled Fusion* 2nd edn (Berlin: Springer)
[19] Magee R M, Galante M, Carr E, Lusk G, McCarren D W and Scime E E 2013 *Phys. Plasmas* **20** 125311

[20] Lieberman M A and Lichtenberg A J 1994 *Principles of Plasma Discharges and Materials Processing* 2nd edn (New York: Wiley)

[21] COMSOL MULTIPHYSICS Version 5.6 (available at: www.comsol.com/)

[22] Santoso J 2018 *PhD Thesis* The Australian National University

[23] Shinohara S, Kaneda N and Kawai Y 1998 *Thin Solid Films* **316** 139–47

[24] Chabert P and Braithwaite N 2011 *Physics of Radio-Frequency Plasmas* (Cambridge: Cambridge University Press)

[25] Janev R K, Reiter D and Samm U 2003 Collision processes in low-temperature hydrogen plasmas (availale at: https://juser.fz-juelich.de/record/38224/files/Juel_4105_Reiter.pdf)

[26] Agnello R 2020 Negative hydrogen ions in a helicon plasma source *PhD Thesis* EPFL

[27] Bacal M 2006 *Nucl. Fusion* **46** S250–9

[28] Pigarov A and Krasheninnikov S 1996 *Phys. Lett.* A **222** 251–7

[29] Taccogna F, Schneider R, Longo S and Capitelli M 2007 *Phys. Plasmas* **14** 073503

[30] Horáček J et al 2004 *Phys. Rev. Lett.* A **70** 052712

[31] Celiberto R, Janev R K, Laricchiuta A, Capitelli M, Wadeha J M and Atems D E 2016 *New J. Phys.* **18** 125005

[32] Krstic P S and Schultz D R 2003 *J. Phys. B* **36** 575–88

[33] Hunter G and Kuriyan M 1977 *Proc. R. Soc. A* **353** 575–88

[34] Fridman A 2008 *Plasma Chemistry* (New York: Cambridge University Press)

[35] Kim Y H et al 1999 *Chem. Phys. Lett.* **314** 1

[36] Rutigliano M 2014 *Plasma Sources Sci. Technol.* **23** 045016

[37] Billing G D and Cacciatore M 1995 *AIP Conf. Proc.* **380** 118

[38] Hiskes J R and Karo A M 1989 *Appl. Phys. Lett.* **54** 508

[39] Taccogna F et al 2021 *Eur. Phys. J. D* **75** 227