Self-consistent fluid model for simulating power coupling in hydrogen ICPs at 1 MHz including the nonlinear RF Lorentz force

Radio frequency (RF) power coupling in inductively coupled plasmas is investigated numerically using a self-consistent fluid model. Hydrogen discharges are simulated at pressures from 0.3–10 Pa and at RF powers of around 1 kW. At the low excitation frequency of 1 MHz a high magnetic RF field of around 30 G is generated by the RF coil, meaning that discharges at low pressures are in the nonlinear skin effect regime. Therefore, a description of the RF power coupling by simple collisional Joule heating is not appropriate. Moreover, models that account for collisionless heating by means of a stochastic collision frequency or as diffusion of the RF current density (as is state of the art for discharges operated in the anomalous skin effect regime at higher frequencies of e.g. 13.56 MHz) are incapable of describing the RF power coupling in the nonlinear skin effect regime properly. This is due to their total neglect or simplified treatment of the RF Lorentz force. Instead, this work demonstrates that the RF power coupling mechanism for discharges operating at low RF in the nonlinear skin effect regime can be described by an electron momentum balance retaining the nonlinear RF Lorentz force as well as electron inertia and advection. The crucial role of the RF Lorentz force in generating the RF plasma current density and thus in shaping the plasma parameter profiles is validated successfully with experimentally obtained electrical and spatially resolved plasma parameters for pressures as low as 0.5 Pa. Below this pressure the results obtained from the model and the ones from the experiment diverge. Most likely this is caused by a sudden change in the electron distribution function at the lowest pressures.


Introduction
Inductively coupled plasmas (ICPs) operated in hydrogen are widely used in the plasma processing industry and in fusion * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. science because of their favorable properties such as easy setting up and low maintenance requirements, as well as good control over the plasma parameters [1,2].
In these ICPs typically only a certain fraction of the radio frequency (RF) generator output power-called RF power transfer efficiency η-is absorbed by the plasma. Consequently, undesired heat is produced by Joule heating in the RF coil and by eddy currents induced in metallic components of the setup. Therefore, considerable experimental [3][4][5][6][7][8] and numerical modeling [9][10][11] efforts have been undertaken to improve the understanding of the RF power coupling mechanism and subsequent optimization of η. In these studies it was found that η depends in a complicated way on external parameters such as the applied frequency, RF generator output power, discharge geometry and RF coil geometry, gas pressure, and on the plasma parameters, in particular their spatial profiles; most importantly on those of the electrons, since they are directly heated by the electric component of the RF field. In order to study how the electrons gain energy from the electric component of the RF field, their power balance has to be solved together (i.e. self-consistently) with Maxwell's equations for the electric and magnetic RF fields. In hydrogen plasmas the electrons lose their energy not only during elastic and inelastic collisions with molecules, where a multitude of processes such as rotational and vibrational excitation are important [12], but also during collisions with dissociated atoms, which exist in non-negligible amounts. In addition, energy losses by convection to the discharge walls are important. This is due to the high surface to volume ratio typically encountered in ICPs operated in the plasma processing industry or in fusion science.
In both applications, discharge sizes of several liters are not uncommon [1,13]. Kinetic methods are too numerically expensive to simulate discharges of this size, for which reason spatially resolved multi-species fluid models are usually used to study different physics aspects of these discharges. In these fluid models, the most simplified treatment of the RF power coupling is to assume a fixed power deposition profile, i.e. the RF fields are not spatially resolved [14,15]. Hence the impact of the RF fields on the plasma parameters and thus the power coupling is not modeled. It is thus more common to spatially resolve the RF fields by solving Maxwell's equations and assuming a local dependency between the electric RF field component and the plasma current density yielding Ohm's law, the development of an RF sheath and collisional Joule heating [16]. However for pressures below 1 Pa collisional Joule heating cannot explain the experimentally observed heating, which is generally accepted to take place via the kinetic collisionless heating mechanism [17,18]. There are fluid models accounting for this heating mechanism either by introducing a stochastic collision frequency [19] or by introducing diffusion of the RF current density [20,21]. The derivations of these approaches are closely related to the anomalous skin effect regime, where the nonlinear RF Lorentz force is usually neglected [22]. At higher radio frequencies (e.g. at the industry standard frequency of 13.56 MHz) this assumption is well justified, because here the typical RF magnetic field strengths are rather small. However, at lower driving frequencies and thus comparatively high magnetic RF field strengths experimental results obtained by Godyak et al [23] show that the nonlinear RF Lorentz force generates higher harmonics in the induced RF current density and in the plasma parameters. This in turn strongly affects the discharge dynamics, wherefore the RF Lorentz force has to be modeled properly to produce results that agree with the experimental measurements, as stated in [21,22,24,25].
Hence a self-consistent fluid model is introduced. Its aim is to describe the RF power coupling in hydrogen discharges that operate in a regime where the RF Lorentz force cannot be neglected, i.e. at pressures from 0.3-10 Pa and RF generator output powers of around 1 kW and, most importantly, at a low RF of 1 MHz. This paper is structured as follows: section 2 provides a brief review of the different RF skin effect regimes and their relation to the power coupling mechanism. Section 3 introduces the experimental ICP setup CHARLIE that operates at 1 MHz. It is used to validate the self-consistent model that is derived in section 4. Section 5 validates the model and discusses the role of the nonlinear RF Lorentz force in sustaining discharges. Final concluding remarks are made in section 6.

RF power coupling regimes
Froese et al have shown that beyond the local and anomalous skin effect regimes [26] there exists the 'nonlinear regime'. The three regimes can be distinguished depending on the pressure p and on the magnetic RF field strength B RF [22]. To show this, plasma and electrical parameters measured in the CHARLIE experimental setup (as described in section 3) are used for estimation of the boundaries between the different regimes as well as for the experimental operating points in figure 1. To establish the relation between the electron temperature and the neutrals pressure (that is plotted on the abscissa), a global model particle balance is used (cf [15]). Note that depending on the ICP geometry (see e.g. ferromagnetically enhanced ICPs), plasma can also be sustained in regions of negligible magnetic RF fields [25,27]. However, this special case is not in the focus of this work, wherefore it is not included in figure 1.
At high pressures above p ≈ 1 Pa the discharges tend to be in the local regime. In this regime, electron heating is well described by collisions. Ohm's law is derived from the electron momentum balance, where the Lorentz force and advection are neglected. This yields a local RF conductivity that connects the RF plasma current density to the RF electric field. In this local regime, Joule heating fully accounts for the power that is absorbed by the plasma [1].
The anomalous skin effect regime is encountered at rather low magnetic RF field strengths, for which reason the nonlinear RF Lorentz force is negligible. For instance, discharges at p = 0.3 Pa with B RF < 0.6 G are in this regime. Here electrons are primarily heated by a wave particle interaction known as collisionless heating [18]. This mechanism can be described by quasilinear theory [1], where collisionless heating is modeled as diffusion of the RF current density caused by viscosity. Two limits can be considered: (i) the limit of locality, for which Chang and Bose presented an expression for the viscosity [28] and (ii) the opposite limit of non-locality, for which Hagelaar derived an expression [20].
For higher magnetic RF field strengths (at p = 0.3 Pa for B RF > 0.6 G), there exists a third skin effect regime, which Froese et al termed the 'nonlinear regime' [22]. This regime differs from the other two in that here the magnetic RF field strength becomes non-negligible. Consequently, the nonlinear Lorentz force should play an important role in facilitating the generation of the RF plasma current density and also contribute Three different RF skin effect regimes depending on the pressure p and on the magnetic field amplitude B RF . Electrical and plasma parameters from the investigated experiment CHARLIE (described in section 3) are used to obtain the regime boundaries by equating the different approximations for the skin depths, as suggested in [22]. Plotted along are typical experimental operating points at a fixed applied frequency f RF = 1 MHz and 800 W generator output power.
significantly to sustaining the discharge and shaping the spatial profiles of the plasma parameters. Interestingly, at even higher magnetic RF field strengths (at p = 0.3 Pa for B RF > 60 G), the Larmor radius of the electrons becomes so small that particles start to behave as if they are once again in the local regime.
Plasma and electrical parameters measured in the CHAR-LIE experimental setup (as described in section 3) are used in figure 1 to show that a typical discharge at an applied frequency of 1 MHz and a generator output power of 800 W are within the local regime for pressures above 1 Pa, whereas they are in the nonlinear regime for pressures below 1 Pa. Figure 2 is a schematic view of the experimental setup CHAR-LIE used for validating the RF power coupling model as described in section 4. The discharge consists of a cylindrical quartz tube with an inner radius of 4.5 cm and an axial length of 40 cm. At z = ±20 cm, stainless steel transition pieces are installed that effectively serve as end plates for the discharges. A helical copper RF coil with 5 windings is connected via a matching network to an RF generator that produces output powers of up to 1 kW at a RF of 1 MHz. A constant inflow of 5 sccm of molecular hydrogen is maintained and by changing the speed of the vacuum pump the pressure can be varied between 0.3 and 10 Pa.

Experimental setup
The experimental setup is equipped with several diagnostic devices for determination of electrical and plasma parameters. The generator output power P generator and the RF coil current amplitude I 0 are independently determined using a V -I probe and a current transformer respectively. A subtractive method is applied to determine the RF power transfer efficiency Here P plasma is the power that is absorbed by the plasma and R network characterizes the losses in the RF coil and the RF network. The efficiency η is determined using the method described in [4,7]. During all measurements perfect matching is ensured by adjusting the matching capacitors C 1 and C 2 , i.e. P generator equals the forward power and the reflected power is zero. The Balmer lines and the Fulcher band are measured with optical emission spectroscopy (OES) along the radial line-ofsight indicated in figure 2. The collisional radiative models Yacora H and H 2 are then used to determine electron density n e , temperature T e , atomic and molecular particle densities n a and n m and the temperature of the H 2 molecules T m [29,30]. In the evaluation, opacity effects are considered by means of the escape factor method [31]. By shifting the OES line-of-sight in the z direction, axial plasma parameter profiles are obtained, where each point of the profile is radially line-of-sight averaged.

RF power coupling model
The simulation domain is based on the experimental ICP setup as described in section 3. As indicated in figure 3, cylindrical symmetry ∂/∂ϕ = 0 is used, i.e. every variable depends only on the two spatial coordinates r, z and not on ϕ. For a vector valued dependent variable u, this means in component notation u = u r (r, z, t)e r + u ϕ (r, z, t)e ϕ + u z (r, z, t)e z , where every component also depends on the time t. To further reduce the numerical effort by 1/2, axial symmetry in the plane that is spanned by the e r and e ϕ axes at z = 0 is exploited, i.e. only z 0 is considered. To focus on the self-consistent RF power coupling mechanism, a uniform background of neutral atoms and molecules throughout the whole discharge is assumed. For a specific pressure, the input number densities n a and n m as well as the temperature T m are obtained from OES evaluation, as described in section 3. It is further assumed that the atomic temperature T a = T m .
The fluid approach involves making assumptions about the velocity distribution function of the charged species (most importantly about the electron velocity distribution function) for closure purposes [32], for deriving boundary conditions and for calculating reaction rate coefficients for the various elastic and inelastic processes. Experimentally validated kinetic simulations of the CHARLIE experiment revealed an electron velocity distribution function of Maxwell-Boltzmann type for pressures larger than 0.5 Pa [33]. This assumption is used here as well.

Gas discharge model
In the region occupied by the plasma, particle balances are solved for each charged species yielding its respective number  [29,30]. Reproduced from [7]. © IOP Publishing Ltd. All rights reserved.
where the notation ∂ t = ∂/∂t is used for brevity. Negative hydrogen ions H − are not included in the model because the process of electron detachment by electron impact is very efficient [12] at the typical T e > 5 eV in the CHARLIE experiment and consequently, n H − n j . Nine processes are considered in the volume reaction rates R j for the creation and destruction of each species: ionization of H, (non-dissociative) ionization and dissociation of H 2 , dissociative excitation and recombination of H + 2 and H + 3 , and creation of H + 3 via the collision of H + 2 with H 2 . The respective rate coefficients for these inelastic collisions are taken from [12]. This selection of processes is well suited to describing hydrogen plasmas in the investigated pressure and power range and has already been used in [15,34,35].
For each of the three positive ion species i ∈ {H + , H + 2 , H + 3 } a momentum balance is solved for the fluid velocities u i = u i,r e r + u i,z e z , i.e.
where m i denotes the ion mass, e the elementary charge and k B Boltzmann's constant. The advection term m i n i (u i · ∇)u i is retained since at low pressures it becomes comparable to the forces per unit volume that act on the ion fluid elements. These are due to the electrostatic field E = E r e r + E z e z , the pressure gradient, where the ion temperature T i measured in units of Kelvin is assumed to be equal to the experimentally obtained T m , elastic collisions between ions and neutrals (μ i,n being the reduced mass and ν i,n the momentum transfer frequency) and the momentum loss associated with a newly created ion m i S i u i , where S i comprises only the positive contributions of R i . As a boundary condition, it is assumed that the ions reach the wall with an effective velocity u i,eff , i.e.
that is calculated by assuming that the velocity distribution function of the ions is Maxwellian with a drift u i,⊥ toward the walls. From this follows where erf denotes the Gauss error function and v i, p = 2k B T i /m i 1/2 the most probable speed.
The electron momentum balance looks similar to that of the ions but with the important difference that the electrons are-in contrast to the ions-magnetized by the RF magnetic field. This is why the Lorentz force (6) Inertia and advection can be neglected for the r and z components of the electron momentum balance [28]. Hence, these components can be simplified to a drift-diffusion approximation for the electron fluid velocity with the electron mobility μ e = e/m e ν e,n , the electron mass m e and the momentum transfer frequency between electrons and neutrals ν e,n .
To derive a boundary condition for the electron velocity, a Maxwellian velocity distribution function without drift is assumed. From this follows, in analogy to equation (5), with u e,⊥ = 0 where v e, p = 2eT e /m e 1/2 . An energy balance is solved for the electron temperature T e , measured in electron volts where the total electron energy flux Q e is due to convection and conduction: The convective energy flux is affected by the RF magnetic field via F L in equation (7). For the conductive heat flux, a separate transport equation is solved [36]: In its r and z components the time derivative can be dropped, as was already done in the momentum balance. Using this approximation, equation (11) is solved for q e : where the thermal conductivity coefficient is defined as n e eT e m e ν e, n .
Note that in analogy to the electron fluid velocity from equation (7) there is a contribution q RF to the heat flux due to the RF magnetic field. This contribution is calculated in section 4.3. On the right-hand side of the electron energy balance equation (9), E RF is the source term describing the RF heating, as defined in section 4.2 and equation (22). The losses E loss are due to elastic collisions of electrons with the neutral species as well as due to inelastic processes such as electronic excitation and ionization of H as well as ionization, dissociation, singlet and triplet state and vibrational excitation of H 2 , where excitation of the first and second vibrational levels are included. The rate coefficients for the elastic processes are taken from [37,38] while those for the inelastic processes are taken from [12]. Electrons lose energy in sustaining an electrostatic field if u e · E > 0, or they gain energy by Joule heating, if u e · E < 0. As a boundary condition for the total electron energy flux, a Maxwellian velocity distribution function without drift toward the wall yields To obtain the electrostatic potential φ, that results from the space charge, Poisson's equation is solved within the plasma domain. Here, 0 is the vacuum permittivity. The electrostatic field is calculated as E = −∂ r φe r − ∂ z φe z . To ensure proper development of an electrostatic plasma sheath, φ = 0 is used as a boundary condition at the dielectric discharge walls. Even though this boundary condition is only strictly true for a metallic wall, it has been shown that its use for dielectric walls results in only minor changes to the spatial profiles of the plasma parameters inside the discharge [39].

Electromagnetic model
Maxwell's equations are solved for the electric E RF and magnetic B RF components of the RF fields in the plasma and in the surrounding air domain. The plasma current density is denoted by J RF . At the low applied frequency of 1 MHz it is necessary to resolve the discharge dynamics during the RF period, as was demonstrated by Kolobov and Godyak [25]. For this reason Maxwell's equations are solved in the time domain. In the Ampère's law equation (16) the quasi-magnetostatic approximation is used, i.e. the displacement field is neglected [40]. This approximation is well fulfilled in the CHARLIE experiment due to the low driving frequency. The system is excited by a surface current that is applied at the outer surface of the RF coil. The interior of the coil wire is thus excluded from the simulation domain. Assuming a surface current for the excitation is a good approximation because the skin depth in the copper wire is around 65 μm, which is small compared to the diameter of the coil wire, d coil = 6 mm, used here. The amplitude value of the surface current density J surf,ϕ is directed in ϕ direction, i.e.n where a sinusoidal excitation at ω RF = 2π f RF is assumed. As a result of the excitation in ϕ direction and because of the assumed cylindrical symmetry B RF = B RF,r e r + B RF,z e z , J RF = J RF,ϕ e ϕ and E RF = E RF,ϕ e ϕ . This implies that capacitive coupling, i.e. electric RF field components in the r and z directions, is precluded from the simulation. This is a good approximation for n e > 10 16 m −3 [41], which is true for the CHARLIE experiment. J surf,ϕ is related to the total RF current amplitude I 0 by The RF coil does not produce a relevant far field, which explains why the generator output power P generator is either absorbed by the plasma or lost in the RF network. Therefore, the equation must be true. Here is calculated self-consistently from the RF-averaged absorbed power per unit volume The network resistance R network in equation (20) is a measured model input as described in section 3 that accounts for the losses in the RF network. For the CHARLIE experimental setup at f RF = 1 MHz with a five turn coil as depicted in figure 2, its value is R network = 0.09 Ω, as determined in [7]. An integral controller controls the value of I 0 at any time t, i.e.
Here P generator,set is the set generator output power that has to be specified by the user before each simulation run. The integral gain K i is tuned such that the system is smoothly ramped up to the steady state value of I 0 , where the actual generator power P generator,is equals the set value. Here K i has no influence on the steady state that is obtained, rather it only influences behavior during the ignition phase.
The boundary conditions at the simulation domain boundaries are those of a perfect electric conductor, where no electric RF field component parallel to the boundary is allowed, i.e. E RF,ϕ = 0. The dimensions of the simulation domain r = z = 100 cm are chosen, such that the effect of this boundary condition on the RF fields in the plasma and in the air surrounding the RF coil is negligible.

Nonlinear terms
At higher driving frequencies it is often possible to only account for the ponderomotive force (i.e. the RF averaged Lorentz force) [27]. However, at 1 MHz, the full time dependency in all equations has to be retained [25]. The RF plasma current density J RF,ϕ , used in equation (22) to obtain E RF , is calculated from J RF,ϕ = −en e u e,ϕ . The RF electron fluid velocity u e,ϕ is determined from the ϕ component of the electron momentum balance. The components of the RF Lorentz force are determined as follows: The r and z components of equation (24) are used in equation (7). The ϕ component of the electron momentum balance equation (6) is obtained as m e ∂ t u e,ϕ + u e,r ∂ r u e,ϕ + u e,z ∂ z u e,ϕ + u e,r u e,ϕ r = −eE RF,ϕ − e(u e,z B RF,r − u e,r B RF,z ) − m e ν e, n u e,ϕ .
Note that herein the nonlinear terms resulting from the Lorentz force as well as the ones from the advection are retained. The r and z components of q RF in equation (12) Considering the ϕ component of equation (11) yields m e ∂ t q e,ϕ = −e(q e,z B RF,r − q e,r B RF,z ) − m e ν e, n q e,ϕ .
State-of-the-art ICP fluid models that focus on the nonlinear effects at low applied frequencies, such as the one by Si et al [42], account for the RF Lorentz force as well as for inertia and advection in the electron momentum balance. However, the RF contribution to the heat flux, as given by equations (12) and (27), is to the knowledge of the author retained for the first time in this work.

Implementation
Weak formulations of the equations in section 4 are derived and implemented in COMSOL Multiphysics R that is used for mesh generation and for solving the system of partial differential equations (PDE) [43]. The typical mesh cell length in the plasma domain is h ≈ 10 −3 m and there are smoothly decreasing boundary layers down to 10 −4 m directly at the walls to ensure proper plasma sheath development. For the time stepping, an implicit backward-difference formula is used. In comparison with explicit time stepping methods, this has the important advantage that no Courant-Friedrichs-Lewy condition must be obeyed, which would vastly increase the time needed to reach a steady state that is typically obtained at around 15 RF cycles. However, since the time stepping has to properly resolve the RF time scale, time steps are typically restricted to 10 −8 s. The resulting linear systems at each time step are solved with a direct solver.
The ion particle and momentum balances are dominated by advection and therefore have to be stabilized. This is done by adding numerical diffusion to the ion particle and momentum balances. It is ensured that the amount of numerical diffusion is kept as low as possible such that the ion densities and fluxes are not affected by it. The electron momentum balance with its nonlinear terms is stabilized in a similar fashion. An exponential formulation for all particle densities is used to avoid negative densities that can appear for numerical reasons especially near the walls where steep gradients are common during the plasma ignition phase. The integral controller from equation (23) is rewritten as an ordinary differential equation that is solved along with the PDE system. Linear shape functions are used for all dependent variables resulting in typically around 60 000 degrees of freedom. The model is solved on a 6-core workstation (Intel64 Family 6 Model 62 processors running at 3.5 GHz) with 128 GB RAM. One run takes approximately six hours to reach a steady state solution.

Model validation
The model inputs and outputs are summarized in table 1, where the values of the calculated output quantities are obtained as a steady state solution for the model input parameters that are obtained from experimental evaluations. Figure 4 shows typical spatial distributions of the steady state quantities (averaged over one RF period) of the RF magnetic field, electron density, Lorentz force and electron temperature. These values are obtained from the model at a pressure of 1 Pa, i.e. n m = 6.8 × 10 19 m −3 , n a = 3.0 × 10 19 m −3 and T m = T a = 740 K, P generator = 800 W and R network = 0.09 Ω. In this case, the plasma absorbs P plasma = 705 W and thus η = 0.88. An RF coil current amplitude of I 0 = 46 A develops generating a high magnetic RF field amplitude of around 20 G directly at the RF coil. The magnetic RF field (shown in the first plot) decreases radially inwards in the plasma to values of around 10 G in the center of the discharge. Further away from the RF coil for z > 10 cm, the RF field strength is almost zero. Therefore, electrons are only strongly magnetized in the vicinity of the RF coil. The electron density (shown in the second plot) has its maximum value of around 4.5 × 10 17 m −3 in the center of the discharge and the density decreases sharply in a radially outward direction. The main reason for this sharp decrease is the Lorentz force, which is shown in the third plot. In the region of the sharp density decrease, F L compensates the radially outward pointing pressure gradient force, for which reason the plasma is pushed away from the RF coil. The RF magnetic field also affects the heat flux via q RF in equation (12), which mainly impacts on the electron temperature profile (shown in the fourth plot). T e has its maximum value of around 9 eV at the discharge wall in the vicinity of the RF coil and it sharply decreases radially inwards to a value of around 6 eV. Note that quantities such as the plasma potential, electron temperature and density as well as the axial and radial current densities have a DC component as well as an oscillating one including several higher harmonics. These are generated by the RF Lorentz force, as shown by Si et al [42]. For the quantities above the second harmonic dominates over all higher harmonics in agreement to what was found in [42]. Figure 5 shows at a fixed P generator = 800 W for pressures in the range 0.3-10 Pa a comparison between the experimentally obtained and modeled steady state parameters RF power transfer efficiency η, electron temperature T e and electron density n e . In all cases shown in this section, T e and n e are averaged in time over one RF period as well as in space along the radial line-of-sight at z = 0, as depicted in figure 2. To show the importance of the nonlinear RF effects induced by the Lorentz force F L , as given by equation (24), by the advection term in the electron momentum balance equation (25) and by the q e × B RF term, as given by equation (26) in the heat flux balance, numerical calculations with and without including these terms are performed. I.e. the above terms are either included (for brevity labeled as 'model' in figure 5) or set to zero (labeled as 'model w/o RF effects' in figure 5).

Pressure variation
For pressures above 3 Pa, whether or not the RF effects are taken into account yields similar quantitative results: as the pressure increases, η and T e decrease, whereas n e remains substantially unchanged. Note that the trend of an almost invariant n e fits better with the experiment when the RF effects are included. The decrease in η is more pronounced in the model than in the experiment. The resulting absolute difference in η of roughly 10% at 10 Pa is not fully understood yet because the absolute values of n e fall within the error bars and those of T e form a close fit. For the model to produce values of η as high as in the experiment, an unrealistically high n e of 7 × 10 17 m −3 would be necessary at 10 Pa.
For pressures below 1 Pa, the qualitative behavior of the quantities obtained from the model and the ones from the experiment are again similar: as the pressure decreases, η and n e decrease, whereas T e increases. However, quantitatively, the model results change substantially: when the RF effects are taken into account, solutions obtained at 1 and at 0.5 Pa produce values of η, T e and n e in good agreement with the experimentally measured values. Note that at low pressures, the experimental error for the evaluation of T e is quite high because the rate coefficients for the processes that are used in the evaluation of the OES are rather flat for T e > 10 eV. When the RF effects are not taken into account, a steady state solution of the model is only possible for pressures as low as 0.55 Pa.
Here T e is already close to 20 eV, so that it is scarcely possible for the electron energy losses that are mainly at the surfaces to be balanced by the RF input power. This coincides with an electron density of around 10 16 m −3 , i.e. almost one order of magnitude smaller than the experimentally obtained electron density.
For pressures below 0.4 Pa, the model is not able to reproduce the experimental trends of the plasma parameters and η. This discrepancy might be caused by the RF Lorentz force,  which is known to affect the electron distribution function at these low pressures by depleting slow electrons, as was measured by Godyak et al in a similar experimental setup [23].

Axial electron density profiles
To further validate the model, the axial n e profiles obtained with and without including the RF effects are compared to the experimental profile. Accordingly, the steady state values of n e are averaged in time over one RF period and in space along the radial line-of-sight at each axial position z. z = 0 is located at the symmetry plane of the RF coil, as depicted in figure 2. This comparison is performed at P generator = 800 W and p = 1 Pa. The left plot in figure 6 shows the relative axial decrease n e /n e (z = 0). Without including the RF effects in the model the decrease in the electron density is not sufficiently pronounced compared to the decrease in the experimentally measured electron density for almost all z. On the other hand, including the RF effects yields a relative decrease that is well comparable to the experimental one. The right plot of figure  6 shows the absolute value of n e (z). On the one hand, the increase in n e toward z = 0 is almost not apparent when the RF effects are not included. Therefore the absolute values of n e are not in agreement with the experimentally obtained values for z < 11. As already shown in figure 5, this missing increase in n e results in a situation where for p 0.55 Pa the electron density and thus the RF power transfer efficiency become so low that no steady state solution can be achieved. On the other hand, including the RF effects reproduces the increase in n e toward z = 0.

RF coupling mechanism
The model fails to reproduce the experimentally observed electrical and plasma parameters at low pressures when the RF effects are not included, as shown in figures 5 and 6. What is  commonly done in this case is to include collisionless heating at pressures below 1 Pa. Two different approaches, both closely related to the anomalous regime, are typically followed in the literature. One, proposed by Vahedi et al [19], uses a stochastic frequency ν st in a local expression for the RF plasma conductivity. One crucial assumption in this approach is that the RF fields decrease exponentially in an RF skin layer δ that has to be much smaller than the discharge radius. As is shown in the plot of the RF magnetic field at 1 Pa, in figure 4, the RF skin effect is not fully established, i.e. the magnetic component of the RF field does not fully decay radially inside the discharge and the skin depth δ ≈ r plasma = 4.5 cm. An estimation of the collisionless skin depth (as defined in [19]) using experimentally obtained plasma parameters confirms this result together with a stochastic frequency that is roughly equal to the applied frequency, i.e. ν stoc ≈ ω RF at 1 Pa. This was also found by Rauner et al [7]. Since the decrease of the electric field is thus also not exponential in the negative radial direction (but rather linear due to the assumed cylindrical symmetry), it follows that in this case ν st cannot be used in a spatially resolved fluid model to calculate the electric component of the RF field.
An approach for including collisionless heating in spatially resolved fluid models was proposed by Hagelaar [20], where diffusion of the RF current density caused by thermal motion of the electrons is described by an effective viscosity coefficient μ eff that is incorporated in the RF component of the electron momentum balance. Using this approach in the model (regardless of whether the RF effects are included or not) has the outcome that no steady state solution can be found because it appears that the amount of diffusion introduced by the effective viscosity coefficient μ eff 5 × 10 4 m 2 s −1 is very high, which in turn increases the electric RF field needed to sustain the discharges at low pressures. The overestimated diffusion of the RF current density might be due to the underlying assumption of this approach that the RF Lorentz force is small, which is invalid in the nonlinear regime.  RF field amplitude, i.e. μ eff ∝ B x RF , where x < 0 in the nonlinear regime rather than with an inverse power law of the RF frequency μ eff ∝ ω −1/3 RF , as is the case in the anomalous regime. However, the author is not aware of any known analytic method for determining the effective viscosity coefficient in the nonlinear skin effect regime. The fact that the model reproduces the experimental parameters at low pressures when the RF effects are included, is a first indication that the influence of the effective RF current density diffusion might be rather small compared to the RF effects.

RF Lorentz force and RF-magnetized heat flux
As is evident from the pressure variation shown in figure 5, for pressures below 3 Pa the nonlinear RF effects have an important impact on the electrical and on the plasma parameters. It has been checked by excluding the RF-magnetization of the heat flux (i.e. by setting q e × B RF from equation (26) to zero) that the major influence on the electron density and on the RF power coupling stems from the Lorentz force and the advection, whereas the RF-magnetized heat flux mainly affects the electron temperature profile by creating a transport barrier for the heat flux in the radial direction and thus steepening the radial electron temperature profile, as shown in the fourth plot in figure 4. The focus of the discussion in the remainder of this paper is thus on the Lorentz force and not on the RF-magnetized heat flux.
The r and z components of the Lorentz force change the electron momentum balance in the vicinity of the RF coil, as a result of which the resulting electron flux is significantly altered in the whole discharge. This can be seen by comparing the two plots in the upper row of figure 7 that are obtained at P generator = 800 W and p = 1 Pa without (upper left plot) and with (upper right plot) including the RF effects. When the RF effects are not included, the electrons flow counterclockwise in a semicircle. Including the RF effects reverses the direction of the electron flux and additionally a vortex forms in the vicinity of the RF coil. Note also that the magnitude of the electron flux increases by almost two orders of magnitude, when the RF effects are included. The theoretical possibility of non-ambipolar vortex type electron fluxes in ICPs was also discussed by Shivarova et al [44]. The radially inwards pointing electron flux in the vicinity of the RF coil due to the r component of the Lorentz force compresses the plasma and causes an increase in n e . The spatial profiles of n e are shown in the lower row of figure 7. Here the steepness of the electron density profile in radial direction increases when the RF effects are included (lower right plot). Quantitatively, this yields a value of the electron density in the discharge center that is around seven times higher than when the RF effects are not included (lower left plot). The increased steepness of the density profile in axial direction that was already observed in figure 6 is a consequence of the radial density increase that is caused by the Lorentz force.

Conclusion
A self-consistent fluid model for description of the RF power coupling in hydrogen ICPs at a low excitation frequency of 1 MHz has been established. The model is capable of reproducing experimentally obtained trends of the parameters RF power transfer efficiency as well as electron density and temperature for ICPs operated at power levels of around 1 kW in a broad pressure range from 0.5-10 Pa.
Due to the low excitation frequency, discharges at pressures below 1 Pa are shown to be in the nonlinear skin effect regime, where the nonlinear Lorentz force is important for the RF power coupling. Hence the model focuses on the selfconsistent description of the Lorentz force in the electron momentum balance. It has been shown that due to a high magnetic RF field and thus a high plasma current density in this regime the RF Lorentz force plays an important part in shaping the electron flux and density profiles by compressing the plasma in the radial direction toward the discharge axis. For pressures below 0.4 Pa, the results obtained with the model substantially deviate from the experimentally obtained ones. A likely cause is that the electron distribution function becomes non-Maxwellian at these lowest pressures. Hence, to investigate discharges at pressures below 0.4 Pa, a kinetic description of the electrons might be necessary.

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