Kinetic analysis of negative power deposition in inductive low pressure plasmas

Low temperature plasmas have found widespread use in materials processing. Particularly inductively coupled plasmas (ICPs), which are in the spotlight of this paper, play a crucial role for a number of applications, e.g., for anisotropic etching of sub-micron patterns [1–5]. Surface modifications like anisotropic etching require that the impinging ions do not significantly experience collisions with the neutrals of the background gas, in particular in the region above the substrate or wafer. This can only be achieved in the low gas pressure regime, where the gas pressure is not higher than approximately 10 mTorr. In this so called non-local regime the local relationship between current density and electric field in the plasma described by a simple Ohm's law given by

$$\begin{eqnarray}&&\vec{j}(\vec{r},t)=\sigma (\vec{r},t)\vec{E}(\vec{r},t)\end{eqnarray} \tag{ 1 }$$

is not valid anymore. The reason is that the thermal motion of electrons is assumed to be negligible and that the plasma is heated by local collisional dissipation of electromagnetic energy. The collisions of electrons with the background gas perturb their regular motion which leads to frequent randomization, effectively erasing the electrons' memory. The phase is reset on every collision encounter and there is no phase shift between the current density and the electric field in between two encounters. The heating therefore results in local Ohmic power dissipation.

The situation at low gas pressures is completely different. Since the electron mean free path is comparable or even larger than the size of the plasma itself, collisions of electrons with the background gas become rare. Consequently, local collisional heating ceases to be an effective mechanism. The power dissipation is mainly due to collisionless electron heating [6–9]. A phase shift establishes due to the finite thermal velocity of the electrons. That is, thermal electrons move through the plasma and gain electromagnetic energy from the electric field particularly in the confined region close to the boundary, i.e., the so-called skin layer. If the thermal velocity of the electrons is high enough as to pass through the skin layer on a time scale short compared to the time scale of the oscillating electromagnetic field, then the moving electrons may gain net energy from the electromagnetic field. The required randomization of the regular electron motion is again due to collisions, however, the heating is non-local.

In a non-local picture, one can find that the electrical field $\vec{E}$ at the position ${\vec{r}}^{\prime }$ inside the plasma at an instant of time in the past ${t}^{\prime }$ influences the current density $\vec{j}$ at the position $\vec{r}$ at present time t. The local version of Ohm's law has to be rewritten

$$\begin{eqnarray}&&\vec{j}(\vec{r},t)={\int }_{V}{\int }_{-\infty }^{t}\sigma (\vec{r},{\vec{r}}^{\prime },t-{t}^{\prime })\vec{E}({\vec{r}}^{\prime },{t}^{\prime }){{\rm{d}}}^{3}{r}^{\prime }{{\rm{d}}{t}}^{\prime }.\end{eqnarray} \tag{ 2 }$$

Here, σ is now a distributed conductivity containing all the necessary information from all past instants of time and at all positions in space. It is thus both non-local in space and non-local in time allowing for causality. In fact, this is the most universal way to couple two vector fields by a scalar field.

It has been long acknowledged that the spatial and temporal dispersion of the plasma conductivity is due to a 'warm plasma' (i.e., the thermal motion of electrons), analogous to the anomalous skin effect in metals [10, 11]. In the context of plasmas the anomalous skin effect is first discussed by Weibel [12]. Certain peculiarities are observed in both experiments and theory: the distribution of the electromagnetic field shows a non-monotonic decrease, rather than an exponential decrease observer for the classical skin effect, so that the appearance of local minima and maxima of the electromagnetic field  is  observed. The most peculiar phenomenon reported is, however, negative power deposition. This means that the electrons gain kinetic energy from the electric field at some position and time and later return a fraction of their kinetic energy to the electrical field at a different position [13–18].

The anomalous skin effect in gaseous plasmas is particularly studied using one-dimensional models for both, bounded and infinite domains. Most results reported for bounded plasmas assume a uniform plasma density. Kaganovich et al developed self-consistent one-dimensional models in order to study the power deposition and electron heating in ICPs [19–22]. Particularly, the effect of the electron energy distribution function on the power deposition and the plasma density is studied by means of a self-consistent one-dimensional model [23, 24]. Quite recently, Hagelaar  reported  on a fluid description of non-local electron kinetics in ICPs in order to  obtain  scaling laws for the non-monotonic spatial structure of the electromagnetic field, the anomalous skin effect, and a condition for the appearance of local negative power deposition [25, 26].

To the best of our knowledge, non of  these  approaches provide a self-consistent analytic solution for a one-dimensional bounded domain, unsymmetrical with respect to the boundary conditions (i.e., with excitation only from one side). In this work we focus exactly on this problem. We analyze the anomalous skin effect and negative power deposition in an unsymmetrical plasma at low  pressures  using an analytical kinetic model of a finite plasma slab. The paper is organized as follows: we first revisit the basic equations derived from consistently coupled Maxwell's equations and Boltzmann's equation, and motivate the spatially one-dimensional description of the problem. We formulate a boundary value problem of Sturm–Liouville type for the electron distribution function and the penetrating electric field. Boltzmann's equation is solved using an analytical ansatz similar to the approach proposed by Tyshetskiy et al [27]. However, we do not assume an exponentially decreasing electrical field in an infinite plasma, but we develop self-consistent analytical solutions to both problems—the distribution function and the electrical field—by means of Hilbert space methods. We finally discuss the results for certain generic plasma parameters (for a low-pressure case) and compare them with results from particle in cell simulations. Particularly, we discuss phase mixing, negative power deposition, and the effect of a finite thermal electron velocity on the distribution of the electrical field and the power deposition.

In order to describe the anomalous skin effect, which is a spatially and temporally non-local effect, Maxwell's equations have to be coupled consistently to Boltzmann's equation for each constituent of the plasma. Of course, this set of equations is not tractable. It has to be simplified by means of a scale analysis.

We assume the angular (radio-)frequency of the penetrating electromagnetic field ω to lie significantly above the ion plasma frequency ${\omega }_{{\rm{pi}}}$ and significantly below the electron plasma frequency ${\omega }_{{\rm{pe}}}$. For the time scales we assume the ordering ${\omega }_{{\rm{pi}}}\ll \omega \ll {\omega }_{{\rm{pe}}}$. The ions are not affected by the rf modulation and react only on the phase-averaged electric field. The radio-frequency current is therefore carried by electrons alone. To describe the electrodynamic effects it is therefore justified to formulate Boltzmann's equation solely for the electrons

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}+\vec{v}\cdot {\rm{\nabla }}f-\displaystyle \frac{e}{m}(\vec{E}+\vec{v}\times \vec{B})\cdot {{\rm{\nabla }}}_{v}f={\left.\displaystyle \frac{\partial f}{\partial t}\right|}_{{\rm{col}}},\end{eqnarray} \tag{ 3 }$$

where f is the distribution function of electrons and m and $-e$ are the electron mass and charge, respectively. The plasma is assumed to be quasi-neutral and (for simplicity) homogeneous, ${n}_{{\rm{e}}}={n}_{{\rm{i}}}=n$. We therefore neglect the plasma boundary sheaths, which are anyhow small compared to the discharge length. The plasma density is given by

$$\begin{eqnarray}&&n=\int {{f}{\rm{d}}}^{3}v.\end{eqnarray} \tag{ 4 }$$

Maxwell's equations governing the electric field $\vec{E}$ and the magnetic field $\vec{B}$ can then be written in the quasi-stationary approximation

$$\begin{eqnarray}&&{\rm{\nabla }}\times \vec{B}={\mu }_{0}\vec{j},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\times \vec{E}=-\displaystyle \frac{\partial \vec{B}}{\partial t},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \vec{B}=0,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \vec{E}=0.\end{eqnarray} \tag{ 8 }$$

Here, the displacement current on the right-hand side of (5) is neglected [28]. This approximation is justified since it scales with ${\omega }^{2}/{\omega }_{{\rm{pe}}}^{2}\ll 1$. In addition, the space charge on the right-hand side of (8) vanishes due to the assumption of a quasi-neutral plasma, since it scales in the same way. The current density is defined  by  the first velocity moment of the distribution function

$$\begin{eqnarray}&&\vec{j}=-e\int \vec{v}{{f}{\rm{d}}}^{3}v.\end{eqnarray} \tag{ 9 }$$

This system of coupled nonlinear partial differential equations is still very difficult to solve. We therefore study the perturbation of the equilibrium state and apply a linearization of the dynamical (and time-harmonic) quantities (indicated by tilde). We set $f=\bar{f}+\tilde{f}$ with the Maxwellian distribution $\bar{f}$, the magnetic field $\vec{B}=\bar{\vec{B}}+\tilde{\vec{B}}=\tilde{\vec{B}}$, and the electrical field $\vec{E}=\bar{\vec{E}}+\tilde{\vec{E}}=\tilde{\vec{E}}$. $\tilde{f}$ is assumed small compared to the equilibrium distribution $\bar{f}$. Since we are interested in the high-frequency dynamics of the system, we set the equilibrium electromagnetic field given by $\bar{\vec{B}}$ and $\bar{\vec{E}}$ to zero. We then obtain a system of linear equations for the dynamical quantities

$$\begin{eqnarray}&&\vec{v}\cdot {\rm{\nabla }}\tilde{f}+(\nu +{\rm{i}}\omega )\tilde{f}=\displaystyle \frac{e}{m}\tilde{\vec{E}}\cdot {{\rm{\nabla }}}_{v}\,\bar{f},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\times {\rm{\nabla }}\times \tilde{\vec{E}}={\rm{i}}\omega e{\mu }_{0}\int \vec{v}\,\tilde{f}{{\rm{d}}}^{3}v.\end{eqnarray} \tag{ 11 }$$

This system of differential equations has to be solved subject to appropriate boundary conditions for both, the distribution function $\tilde{f}$ and the electromagnetic field. In the given system the collisions of electrons with the neutrals of the background gas are taken into account using an effective collision rate for momentum transfer described by ${\nu }_{{\rm{m}}}$. The divergence equations (7) and (8) act as constraints. For the dynamical quantities they translate to

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \tilde{\vec{B}}=0,\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \tilde{\vec{E}}=0.\end{eqnarray} \tag{ 13 }$$

Additionally, we have

$$\begin{eqnarray}&&\int \tilde{f}{{\rm{d}}}^{3}v=0.\end{eqnarray} \tag{ 14 }$$

In order to analytically study the anomalous skin effect and  the  related dissipation mechanisms we investigate a one-dimensional plasma slab as depicted in figure 1. An electromagnetic field given by $\tilde{\vec{B}}=\tilde{B}(x){\vec{e}}_{z}$ and $\tilde{\vec{E}}=\tilde{E}(x){\vec{e}}_{y}$ prescribed at x = 0 penetrates into a semi-infinite bounded homogeneous low-pressure plasma of length L. It is important to notice that the plasma is infinite in the lateral plane, but bounded in x direction. For this one-dimensional case the differential equations become

$$\begin{eqnarray}&&\left(\displaystyle \frac{\partial }{\partial x}+\displaystyle \frac{\alpha }{{v}_{x}}\right)\,\tilde{f}=\displaystyle \frac{\beta }{{v}_{x}}\tilde{E}(x),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}\tilde{E}}{\partial {x}^{2}}={\rm{i}}\omega {\mu }_{0}\,\tilde{j}(x;\tilde{E})\end{eqnarray} \tag{ 16 }$$

with $\alpha ={\rm{i}}\omega +{\nu }_{m}$ and $\beta =-e\bar{f}{v}_{y}/{{mv}}_{{\rm{th}}}^{2}$. The thermal velocity of the electrons is ${v}_{{\rm{th}}}={({T}_{{\rm{e}}}/m)}^{1/2}$ with the electron temperature T_e (in eV; the electron mean velocity is $\bar{v}={(8{T}_{{\rm{e}}}/\pi m)}^{1/2}$). The high-frequency current density $\tilde{\vec{j}}=\tilde{j}(x;\tilde{E}){\vec{e}}_{y}$, which is implicitly a function of the electrical field (and the position) is given by the first moment of the dynamical distribution function $\tilde{f}$

$$\begin{eqnarray}&&\tilde{j}=-e\int {v}_{y}\tilde{f}(x,\vec{v};\tilde{E}){{\rm{d}}}^{3}v.\end{eqnarray} \tag{ 17 }$$

As boundary conditions for (16), we assume a time-harmonic magnetic field B₀ oscillating with angular frequency ω given at x = 0 and let the electric field vanish at the conducting boundary at x = L. The boundary conditions for the electric field $\tilde{E}$ are therefore given by

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial \tilde{E}}{\partial x}\right|}_{x=0}=-{\rm{i}}\omega {B}_{0},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\tilde{E}| }_{x=L}=0.\end{eqnarray} \tag{ 19 }$$

Gathering equations (16) to (19), we have a well-posed boundary value problem for the electric field, i.e., a homogeneous differential equation and inhomogeneous boundary conditions.

Zoom In Zoom Out Reset image size

For practical reasons it is convenient to have homogeneous boundary conditions for the electric field. We thus apply for the electric field the ansatz

$$\begin{eqnarray}&&\tilde{E}={E}_{0}+{E}_{h},\end{eqnarray} \tag{ 20 }$$

with ${E}_{0}={\rm{i}}\omega (L-x){B}_{0}$. The boundary value problem for E_h becomes

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{E}_{h}}{\partial {x}^{2}}={\rm{i}}\omega {\mu }_{0}\,\tilde{j}(x;{E}_{0}+{E}_{h}),\end{eqnarray} \tag{ 21 }$$

subject to homogeneous boundary conditions

$$\begin{eqnarray}&&{\displaystyle \frac{\partial {E}_{h}}{\partial x}| }_{x=0}=0,\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{{E}_{h}| }_{x=L}=0.\end{eqnarray} \tag{ 23 }$$

Since (15) is linear, the right-hand side of (21) can be modified. We can decompose the current density in two components, $\,\tilde{j}(x;{E}_{0}+{E}_{h})={j}_{0}(x;{E}_{0})+{j}_{h}(x;{E}_{h})$, where j₀ is due to E₀ and j_h is due to E_h. The two current densities are then given by

$$\begin{eqnarray}&&{j}_{0}=-e\int {v}_{y}{f}_{0}(x,\vec{v};{E}_{0}){{\rm{d}}}^{3}v,\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{j}_{h}=-e\int {v}_{y}{f}_{h}(x,\vec{v};{E}_{h}){{\rm{d}}}^{3}v,\end{eqnarray} \tag{ 25 }$$

with the corresponding decomposition $\tilde{f}={f}_{0}+{f}_{h}$.

The distribution functions f₀ and f_h, which implicitly depend on the electric field E₀ and E_h, are the solutions to the linearized Boltzmann equation (15) given the appropriate right-hand side inhomogeneity, therefore  

$$\begin{eqnarray}&&\left(\displaystyle \frac{\partial }{\partial x}+\displaystyle \frac{\alpha }{{v}_{x}}\right)\,{f}_{0}(x,\vec{v};{E}_{0})=\displaystyle \frac{\beta }{{v}_{x}}{E}_{0}(x),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\left(\displaystyle \frac{\partial }{\partial x}+\displaystyle \frac{\alpha }{{v}_{x}}\right)\,{f}_{h}(x,\vec{v};{E}_{h})=\displaystyle \frac{\beta }{{v}_{x}}{E}_{h}(x).\end{eqnarray} \tag{ 27 }$$

The (spatial) boundary conditions for the distribution functions are so-called 'reflecting' boundary conditions with

$$\begin{eqnarray}&&{{f}_{0}({v}_{x})| }_{x=0,L}={{f}_{0}(-{v}_{x})| }_{x=0,L},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{{f}_{h}({v}_{x})| }_{x=0,L}={{f}_{h}(-{v}_{x})| }_{x=0,L}.\end{eqnarray} \tag{ 29 }$$

We now have exactly one first order differential equation for each distribution function f₀ and ${f}_{h}$, each of which subject to exactly two boundary conditions, one at either sides. This is in fact a well-posed mathematical problem which can be solved analytically, as shown  in  the subsequent section.

With the knowledge on the distribution functions, and thus with the knowledge on the current densities, we are able to reformulate the desired inhomogeneous boundary value problem for the electrical field (21)

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{E}_{h}}{\partial {x}^{2}}-{\rm{i}}\omega {\mu }_{0}\,{j}_{h}(x;{E}_{h})={\rm{i}}\omega {\mu }_{0}\,{j}_{0}(x;{E}_{0}),\end{eqnarray} \tag{ 30 }$$

subject to homogeneous boundary conditions (22) and (23). This boundary value problem is of Sturm–Liouville type and can also be solved analytically by means of standard Hilbert space methods.

The first step toward the analytical solution of (30) for the electric field is to calculate the current densities j₀ and j_h from the solutions of the kinetic equations (26) and (27). These can be obtained by means of the method of integrating factors. Hence, the general solution to differential equations of the given type is

$$\begin{eqnarray}f(x,\vec{v};E) & = & f({x}_{0}){{\rm{e}}}^{-\alpha (x-{x}_{0})/{v}_{x}}\\ & & +\displaystyle \frac{\beta {{\rm{e}}}^{-\alpha x/{v}_{x}}}{{v}_{x}}{\displaystyle \int }_{{x}_{0}}^{x}{{\rm{e}}}^{\alpha {x}^{\prime }/{v}_{x}}E({x}^{\prime }){{\rm{d}}{x}}^{\prime }.\end{eqnarray} \tag{ 31 }$$

Here, we have one unknown constant $f({x}_{0})$, but two boundary conditions ${f| }_{x=0,L}$. This inherent problem can be overcome using the appropriate ansatz $f={f}^{+}+{f}^{-}$ for the distribution function. ${f}^{+}$ represents the solution for ${v}_{x}\geqslant 0$, while vanishing for ${v}_{x}\lt 0$, and satisfying the boundary condition at x = 0. Accordingly, ${f}^{-}$ is the solution for ${v}_{x}\leqslant 0$, while vanishing for ${v}_{x}\gt 0$, and satisfying the boundary condition at x = L. We consequently obtain for ${f}^{+}$ and ${f}^{-}$, respectively

$$\begin{eqnarray}&&{f}^{+}=\left\{\begin{array}{l}f(0){{\rm{e}}}^{-\alpha x/{v}_{x}}+\displaystyle \frac{\beta {{\rm{e}}}^{-\alpha x/{v}_{x}}}{{v}_{x}}{\displaystyle \int }_{0}^{x}{{\rm{e}}}^{\alpha {x}^{\prime }/{v}_{x}}E({x}^{\prime }){{\rm{d}}{x}}^{\prime }\\ \quad \quad \mathrm{for}\quad {v}_{x}\geqslant 0,\\ 0\\ \quad \quad \mathrm{for}\quad {v}_{x}\lt 0,\end{array}\right.\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{f}^{-}=\left\{\begin{array}{l}f(L){{\rm{e}}}^{-\alpha (x-L)/{v}_{x}}+\displaystyle \frac{\beta {{\rm{e}}}^{-\alpha x/{v}_{x}}}{{v}_{x}}{\displaystyle \int }_{L}^{x}{{\rm{e}}}^{\alpha {x}^{\prime }/{v}_{x}}E({x}^{\prime }){{\rm{d}}{x}}^{\prime }\\ \qquad \mathrm{for}\quad {v}_{x}\leqslant 0,\\ 0\\ \qquad \mathrm{for}\quad {v}_{x}\gt 0.\end{array}\right.\end{eqnarray} \tag{ 33 }$$

The unknown coefficients f(0) and f(L) can now be calculated explicitly from the boundary conditions for ${f}^{+}$ and ${f}^{-}$, which are

$$\begin{eqnarray}&&{{f}^{+}({v}_{x})| }_{x=0,L}={{f}^{-}(-{v}_{x})| }_{x=0,L}.\end{eqnarray} \tag{ 34 }$$

We obtain the two expressions

$$\begin{eqnarray}f(0) & = & \displaystyle \frac{\beta }{(1-{{\rm{e}}}^{2\alpha L/{v}_{x}}){v}_{x}}\left({{\rm{e}}}^{2\alpha L/{v}_{x}}{\displaystyle \int }_{L}^{0}{{\rm{e}}}^{-\alpha {x}^{\prime }/{v}_{x}}E({x}^{\prime }){{\rm{d}}{x}}^{\prime }\right.\\ & & \left.-{\displaystyle \int }_{0}^{L}{{\rm{e}}}^{\alpha {x}^{\prime }/{v}_{x}}E({x}^{\prime }){{\rm{d}}{x}}^{\prime }\right),\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}f(L) & = & \displaystyle \frac{\beta {{\rm{e}}}^{\alpha L/{v}_{x}}}{(1-{{\rm{e}}}^{2\alpha L/{v}_{x}}){v}_{x}}\left({\displaystyle \int }_{L}^{0}{{\rm{e}}}^{-\alpha {x}^{\prime }/{v}_{x}}E({x}^{\prime }){{\rm{d}}{x}}^{\prime }\right.\\ & & \left.-{\displaystyle \int }_{0}^{L}{{\rm{e}}}^{\alpha {x}^{\prime }/{v}_{x}}E({x}^{\prime }){{\rm{d}}{x}}^{\prime }\right).\end{eqnarray} \tag{ 36 }$$

Herewith, a general analytic solution as a function of an arbitrary electric field E(x) and subject to reflecting boundary conditions (34) is found. The respective solutions to (26) and (27) namely f₀ and f_h have to be phrased explicitly in terms of E₀ and the unknown ${E}_{h}$ using expressions (32)–(36).

For the solution of the boundary value problem (30), the current density ${j}_{h}$ can now be calculated explicitly from (25) (except for the integration over ${v}_{x}$). It is convenient to again redefine the current density as

$$\begin{eqnarray}&&{j}_{h}=-{e}{\int }_{-\infty }^{\infty }{ \mathcal L }(x,{v}_{x};{E}_{h}){{\rm{d}}{v}}_{x}\end{eqnarray} \tag{ 37 }$$

with

$$\begin{eqnarray}&&{ \mathcal L }={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{v}_{y}{f}_{h}(x,\vec{v};{E}_{h}){{\rm{d}}{v}}_{y}{{\rm{d}}{v}}_{z}.\end{eqnarray} \tag{ 38 }$$

It is found that the integral over ${v}_{x}$ converges and can thus be calculated by means of standard numerical methods [29].

To obtain the right hand side of (30), the current density j₀ has to be calculated as well. This can be done by replacing the electric field E in the general solution (32)–(36) by the known expression ${E}_{0}={\rm{i}}\omega (L-x){B}_{0}$. From (24) we consequently obtain

$$\begin{eqnarray}&&{j}_{0}=-e{\int }_{-\infty }^{\infty }{ \mathcal R }(x,{v}_{x};{E}_{0}){{\rm{d}}{v}}_{x}\end{eqnarray} \tag{ 39 }$$

with

$$\begin{eqnarray}&&{ \mathcal R }={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{v}_{y}{f}_{0}(x,\vec{v};{E}_{0}){{\rm{d}}{v}}_{y}{{\rm{d}}{v}}_{z}.\end{eqnarray} \tag{ 40 }$$

This expression can be written down explicitly. However, since it is quite cumbersome, we abstain from doing so [29]. Now all constituents of expression (30) are known and the boundary value problem for the electrical field ${E}_{h}$ can be solved.

The boundary value problem for the electric field (30) and (22)–(23) is of Sturm–Liouville type. With the current densities given by (37) and (39), we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{E}_{h}}{\partial {x}^{2}}+{\rm{i}}\omega {\mu }_{0}e{\displaystyle \int }_{-\infty }^{\infty }{ \mathcal L }(x,{v}_{x};{E}_{h}){{\rm{d}}{v}}_{x}\\ &&\quad =-{\rm{i}}\omega {\mu }_{0}e{\displaystyle \int }_{-\infty }^{\infty }{ \mathcal R }(x,{v}_{x};{E}_{0}){{\rm{d}}{v}}_{x}.\end{eqnarray} \tag{ 41 }$$

with homogeneous boundary conditions ${\partial }_{x}{E}_{h}{| }_{x=0}=0$ and ${E}_{h}{| }_{x=L}=0$. Such Sturm–Liouville problems are mathematically well characterized and can be solved by means of standard Hilbert space methods: it is known that (i) the eigenfunctions of a Sturm–Liouville operator form a complete set of orthogonal functions, and (ii) the eigenvalues form a set of non-degenerated positive numbers. For the problem at hand, the assigned eigenvalue problem (differential equation subject to decoupled boundary conditions) reads explicitly

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{u}_{k}(x)}{\partial {x}^{2}}={\lambda }_{k}{u}_{k}(x),\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial {u}_{k}}{\partial x}\right|}_{x=0}={{u}_{k}| }_{x=L}=0.\end{eqnarray} \tag{ 43 }$$

The associated eigenfunctions ${u}_{k}(x)$ and eigenvalues ${\lambda }_{k}$ are defined for integer k and given as follows

$$\begin{eqnarray}&&{u}_{k}(x)=\cos \sqrt{{\lambda }_{k}}x,\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&{\lambda }_{k}=\displaystyle \frac{{(2k-1)}^{2}{\pi }^{2}}{4{L}^{2}}.\end{eqnarray} \tag{ 45 }$$

Since the eigenfunctions form a complete set of orthogonal functions we can expand ${E}_{h}$ into an infinite series

$$\begin{eqnarray}&&{E}_{h}(x)=\displaystyle \sum _{k=1}^{\infty }{a}_{k}\cos \sqrt{{\lambda }_{k}}x.\end{eqnarray} \tag{ 46 }$$

In order to calculate the unknown coefficients ${a}_{k}$ and thus the analytical solution of (41) it is convenient to expand the kernel of the integrals in terms of the same set of orthogonal functions. After plugging in the expansion of the electrical field (46) into the kernel ${ \mathcal L }$ given by (38) we obtain

$$\begin{eqnarray}{ \mathcal L }(x,{v}_{x}) & = & \displaystyle \sum _{k=1}^{\infty }{a}_{k}\left({b}_{k}({v}_{x}){{\rm{e}}}^{-\tfrac{x\alpha }{{v}_{x}}}\right.\\ & & +{c}_{k}({v}_{x})\sin \sqrt{{\lambda }_{k}}x+{{d}}_{k}({v}_{x})\cos \sqrt{{\lambda }_{k}}x).\end{eqnarray} \tag{ 47 }$$

The coefficients ${b}_{k}$, ${c}_{k}$, ${d}_{k}$ are themselves functions of ${v}_{x}$. Finally, this expression for ${ \mathcal L }$ has to be expanded into the eigenfunctions ${u}_{k}(x)$, such that

$$\begin{eqnarray}{ \mathcal L }(x,{v}_{x}) & = & \displaystyle \sum _{k=1}^{\infty }{a}_{k}\displaystyle \sum _{l=1}^{k}[{\underline{b}}_{l,k}({v}_{x})\\ & & +{\underline{c}}_{l,k}({v}_{x})+{\underline{d}}_{l}({v}_{x})\,{\delta }_{l,k}]\cos \sqrt{{\lambda }_{k}}x\end{eqnarray} \tag{ 48 }$$

with an alternative set of velocity ${v}_{x}$ dependent coefficients ${\underline{b}}_{l,k}$, ${\underline{c}}_{l,k}$, and ${\underline{d}}_{l}$. These coefficients can also be written down explicitly, as  given  below. For the kernel ${ \mathcal R }$ in (40) we straightforwardly obtain

$$\begin{eqnarray}&&{ \mathcal R }(x,{v}_{x})=\displaystyle \sum _{k=1}^{\infty }{\underline{r}}_{k}({v}_{x})\cos \sqrt{{\lambda }_{k}}x.\end{eqnarray} \tag{ 49 }$$

Exploiting the orthogonality of the eigenfunctions, we ultimately find a linear system of equations for the unknown coefficients ${a}_{k}$

$$\begin{eqnarray}&&\displaystyle \frac{{\lambda }_{k}{a}_{k}}{{\rm{i}}\omega {\mu }_{0}e}+\displaystyle \sum _{l=1}^{k}\left[{a}_{k}{\displaystyle \int }_{-\infty }^{\infty }{\underline{b}}_{l,k}({v}_{x})\,{{\rm{d}}{v}}_{x}+{\underline{D}}_{l}\,{\delta }_{l,k}\right]\\ &&\quad =-{\displaystyle \int }_{-\infty }^{\infty }{\underline{r}}_{k}({v}_{x}){{\rm{d}}{v}}_{x}.\end{eqnarray} \tag{ 50 }$$

The known coefficients are cumbersome, but are given explicitly by

$$\begin{eqnarray}&&\,{\mathop{b}\limits_{\_}}_{l,k}({v}_{x})\\ &&=-\displaystyle \frac{e\,n\,{v}_{x}^{3}\,\sqrt{{\lambda }_{l}{\lambda }_{k}}\,{{\rm{e}}}^{-{v}_{x}^{2}/(2{v}_{{\rm{th}}}^{2})}\coth (\alpha L/{v}_{x})\sin \sqrt{{\lambda }_{l}}L\sin \sqrt{{\lambda }_{k}}L}{\sqrt{\pi /2}\,L\,m\,{v}_{{\rm{th}}}({\alpha }^{2}+{v}_{x}^{2}{\lambda }_{l})({\alpha }^{2}+{v}_{x}^{2}{\lambda }_{k})},\end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray}&&{\mathop{D}\limits_{\_}}_{l}={\int }_{-\infty }^{\infty }{\mathop{{d}}\limits_{\_}}_{l}({v}_{x})\,{{\rm{d}}{v}}_{x}=-\displaystyle \frac{e\,n\,{{\rm{e}}}^{{\alpha }^{2}/(2{\lambda }_{l}{v}_{{\rm{th}}}^{2})}\,\mathrm{erfc}[\alpha /(\sqrt{2{\lambda }_{l}}{v}_{{\rm{th}}})]}{\sqrt{2{\lambda }_{l}/\pi }\,m\,{v}_{{\rm{th}}}},\end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray}&&{\mathop{r}\limits_{\_}}_{k}({v}_{x})=-\displaystyle \frac{{\rm{i}}\omega {B}_{0}\,{e}\,n\,{{\rm{e}}}^{-{v}_{x}^{2}/(2{v}_{{\rm{th}}}^{2})}[(1+{{\rm{e}}}^{\alpha L/{v}_{x}}){\alpha }^{3}+({{\rm{e}}}^{\alpha L/{v}_{x}}-1){({v}_{x}\sqrt{{\lambda }_{k}})}^{3}\sin \sqrt{{\lambda }_{k}}L]}{\sqrt{\pi /2}\,L\,m\,{v}_{{\rm{th}}}\,{\lambda }_{k}\,{\alpha }^{2}(1+{{\rm{e}}}^{\alpha L/{v}_{x}})({\alpha }^{2}+{v}_{x}^{2}{\lambda }_{k})}.\end{eqnarray} \tag{ 53 }$$

Here, the fact that the symmetric integration of any arbitrary function $g({v}_{x})$ over ${v}_{x}\in (-\infty ,\infty )$ is equal to the integral over its even contribution only has been used. In consequence, the integrals ${\underline{C}}_{l,k}$ vanish after integration of ${\underline{c}}_{l,k}({v}_{x})$ over ${v}_{x}$, because ${\underline{c}}_{l,k}({v}_{x})$ is an odd function in ${v}_{x}$. In addition, the coefficients ${\underline{d}}_{l}({v}_{x})$ have already been integrated analytically in ${v}_{x}$, thus represented by ${\underline{D}}_{l}$.

After numerical integration of the coefficients ${\underline{b}}_{l,k}({v}_{x})$ and ${\underline{r}}_{k}({v}_{x})$ over ${v}_{x}$ and solving for ${a}_{k}$, we obtain ${E}_{h}$, which is the solution to the electric field equation (30) subject to (22) and (23). The total electrical field $\tilde{E}$ can then be calculated from (20) and consequently the magnetic field $\tilde{B}$ and the current density $\tilde{j}$. We have

$$\begin{eqnarray}&&\tilde{E}={\rm{i}}\omega (L-x){B}_{0}+\displaystyle \sum _{k=1}^{\infty }{a}_{k}\cos \sqrt{{\lambda }_{k}}x,\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&\tilde{B}=-\displaystyle \frac{1}{{\rm{i}}\omega }\displaystyle \frac{\partial \tilde{E}}{\partial x},\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&\tilde{j}=\displaystyle \frac{1}{{\rm{i}}\omega {\mu }_{0}}\displaystyle \frac{{\partial }^{2}\tilde{E}}{\partial {x}^{2}}.\end{eqnarray} \tag{ 56 }$$

5.1. Analytical results

In this subsection, the analytic solution to the one-dimensional plasma slab domain depicted in figure 1 is self-consistently evaluated for generic parameters of an ICP discharge. The length of the plasma is set to L = 6 cm. The plasma density and the electron temperature is chosen to be $n=5\times {10}^{11}$ cm⁻³ and 4 eV, respectively. We chose a simple argon plasma model following [5], where the electron-neutral momentum transfer collision frequency is specified as ${\nu }_{m}={K}_{m}{n}_{g}$, with the rate coefficient approximated by ${K}_{m}={10}^{-13}$ m³ s⁻¹ and the neutral gas density as a function of gas pressure ${n}_{g}$(m⁻³) = $3.3\times {10}^{22}p$(Torr). The magnetic field at the interface between the vacuum and the plasma at x = 0 is set to ${B}_{0}={10}^{-6}$ V s m⁻². The angular driving frequency ω and the pressure p (Torr) are parameters which are varied, and thus specified in the respective paragraphs.

Figure 2 shows the magnitude of both, the electrical field (above) and the current density (below) for $\omega =2\pi \ \times 13.56\,\mathrm{MHz}$ and two different gas pressure. In the high pressure regime p = 45 mTorr (i.e., the local regime) indicated by dashed lines an exponential decrease of both fields is observed, which is in fact typical for the classical skin effect. In the non-local regime at p = 1 mTorr, the electrical field as well as the current density reveal a non-monotonic decrease, which is typical for the anomalous skin effect. It is characterized by the appearance of local minima and maxima in the fields' magnitudes.

Zoom In Zoom Out Reset image size

The appearance of the non-monotonic spatial distribution of the fields, translates into the complex phase of the time-harmonic observables. This can be seen in figure 3, where the absolute values (above) as well as the phases (below) of both the electric field (left) and the current density (right) are shown for  varied  gas pressures. For all pressures the more moderate decay of the current density distribution compared to the electric field decay is obvious. This demonstrates the effect of the current diffusion due to thermal electron motion (i.e., a warm plasma effect). The individual phases are themselves directly linked to their penetration depth. The profiles of the phases are correspondingly different for the electric field and current density. The higher the pressure (i.e., more local) the more similar they are. The fact that they are different suggests distinct mechanisms of propagation of the electric field and the current density.

Zoom In Zoom Out Reset image size

It should be recalled that the spatial distributions of the electric field and the current density are comparable for the  classical  skin effect—scaled in magnitude, but identical in phase. This is indicated by figure 3 (red, long-dashed line). Both phases are simply shifted and their phase difference ${\rm{\Delta }}\phi$ remains less than $\pm \pi /2$. It is zero in the limit of purely Ohmic conduction and close to zero in the local regime (see, p = 45 mTorr).

The total power density ${p}_{{abs}}$ deposited can be generally calculated from

$$\begin{eqnarray}&&{p}_{{abs}}=\displaystyle \frac{1}{2}{\rm{Re}}\ \tilde{j}{\tilde{E}}^{* }.\end{eqnarray} \tag{ 57 }$$

In the local regime only positive (Ohmic) power deposition is to be recognized. In contrast, in the non-local regime phase differences ${\rm{\Delta }}\phi$ between electric field and current density with values in the range $\pi /2\lt {\rm{\Delta }}\phi \lt 3\pi /2$ are observed. This is depicted in figure 4 (above). Here the above mentioned range is indicated by the gray region. A phase shift within these limits represents negative power deposition (figure 4, below). Whenever the phase difference crosses the specified interval limits, the power deposition changes its sign. In the logarithmic plot given below, the zero-crossing of the power density diverges. The negative power deposition is clear from a mathematical point of view. The distribution of power density deposited throughout its propagation into the plasma is   given  by ${p}_{{abs}}=1/2\,| \tilde{j}| \,| \tilde{E}| \,\cos {\rm{\Delta }}\phi$. In terms of a polar coordinate representation within the complex plane the real part is clearly negative for the above mentioned phase difference interval. Although clear from a mathematical standpoint, negative power deposition is a quite remarkable physical phenomenon nevertheless.

Zoom In Zoom Out Reset image size

Ohmic power deposition due to collisions of electrons with the neutral gas background is always positive. Given a conductivity σ, the purely Ohmic power density contribution can be approximated. In the limit of a cold collisional plasma (assuming locality), the plasma conductivity is $\sigma ={{\rm{e}}}^{2}{n}_{{\rm{e}}}/{m}_{{\rm{e}}}({\nu }_{m}+{\rm{i}}\omega )$. Correspondingly, from the local version of Ohm's law the power density can in principle be evaluated from either the current density or the electric field. However, given that Ampere's law ties the current density to the curl of the magnetic field irrespective of any material model (i.e., Ohm's law), it seems more reasonable to approximate the Ohmic power density from the non-local current density (in favor of the non-local electric field). It can thus be calculated from

$$\begin{eqnarray}&&{p}_{{\rm{ohm}}}=\displaystyle \frac{1}{2}{\rm{Re}}({\sigma }^{-1})\ \tilde{j}\ {\tilde{j}}^{* }.\end{eqnarray} \tag{ 58 }$$

The power density distribution deposited is given in figure 5 for various gas pressures at a constant frequency of the electromagnetic field of $\omega =2\pi \times 13.56$ MHz. The total deposited power density is represented by solid lines, while the Ohmic contribution is indicated by dashed lines. Clearly, the higher the pressure, the better the local power density calculated from expression (58) approximates the exact non-local power density deposited as obtained from equation (57). For the case of p = 30 mTorr both are nearly the same. Only in the low pressure cases $p\,\leqslant$ 15 mTorr the remarkable feature of negative power deposition is apparent. In these situations, power is removed from certain regions of negative power deposition and transferred to regions of positive power deposition. This implies that a mechanism acts to turn the positive sign of net power deposition negative. This is exactly due to collisionless heating. Electrons gain energy throughout certain periods of time, but experience a randomizing collision only much later after having traversed to a different position within the discharge. The collective effect is clearly non-local in position space. For the most obvious case of p = 1 mTorr, negative power deposition is observed in the interval $1.2\lesssim x\lesssim 3\,\mathrm{cm}$. Of course this phenomenon is only visible for the total power density ${p}_{{abs}}$, the Ohmic power density is ${p}_{{\rm{ohmic}}}\geqslant 0$ in the whole discharge.

Zoom In Zoom Out Reset image size

For all cases it can be seen that the energy is predominantly absorbed near the plasma boundary and that deeper within the plasma alternating regions of negative and positive power deposition appear. The penetration depth as well as the phase structure are non-trivially linked to the degree of non-locality. A phenomenon previously characterized (among others) by [6] who discusses the peculiar parameter space in which collisionless heating efficiently takes place. In particular: (i) that the inhomogeneity within the discharge be large enough allowing the electrons to experience regions of different field amplitude. And more importantly, (ii) it is essential that within one rf period the electron motion   is  sufficiently fast, such that   electrons  traverse a distance larger then the depth of field penetration. This is intrinsically linked to the mean electron energy (e.g., the electron temperature ${T}_{{\rm{e}}}$) and the collision frequency ${\nu }_{m}$.

Because the time for a 'warm' electron to proceed without experiencing a significant change of the electromagnetic field is largely dependent of the rf frequency, a similar phenomenon as for different pressures can be observed for various frequencies ω of the penetrating electromagnetic field while keeping the pressure fixed at p = 15 mTorr. This is depicted in figure 6. Also in this case, the occasions and extent of regions with negative power deposition strongly vary,   depending  on the actual frequency and corresponding phase relation of the penetrating electric field and current density. As apparent from the top left subfigure in figure 6, for the case of $\omega =2\pi \times 6.78\,\mathrm{MHz}$ no negative power deposition is observed. The net power   deposition  is very close to the Ohmic power density   deposited  in local approximation. Only for larger frequencies of $\omega \gtrsim 2\pi \times 13.56\,\mathrm{MHz}$ negative power deposition is present. With increasing frequency the regions of negative power deposition are quenched closer to the driven boundary. These results are in good agreement with previous model results and experimental data [7, 17, 18, 27].

Zoom In Zoom Out Reset image size

5.2. Comparison with particle in cell

For comparison, it is instructive to review the identical discharge setup also in terms of an   alternative  kinetic simulation approach. For this purpose the particle in cell/Monte Carlo collisions (PIC/MCC) approach is utilized, taking into account the full-wave transversal electromagnetic fields (i.e., ${E}_{y}$ and ${B}_{z}$). An extended version of the yapic code, which has been previously benchmarked against other PIC implementations, is utilized [30]. Within the simulation the configuration space is reduced to one dimension, but all three components in velocity space are retained (i.e., 1D3V). The essentials of the proposed simulation model are presented elsewhere [30, 31]. The employed modifications can be summarized as follows: (i) while the transient electron dynamics are maintained, for the ions a homogeneous and constant density is imposed within the plasma slab. (ii) With the same reasoning, imposing the assumption of a quasi-neutral plasma bulk, the longitudinal electric field ${E}_{x}$ is virtually zero using the expression for the ambipolar field $\vec{E}=-\tfrac{{T}_{{\rm{e}}}}{{n}_{i}{\rm{e}}}{\rm{\nabla }}{n}_{i}\equiv 0$ for ${n}_{i}={\rm{const}}$ [32, 33]. (Principally, the longitudinal electric field could have been also obtained from Poisson's equation based on the electrostatic approximation.) Both assumptions assure that the plasma conditions are as similar as possible to the ones assumed within the analytic model. A fortuitous side effect is a greatly diminished statistical noise level that is inherent to PIC. (iii) The full-wave transversal electromagnetic fields E_y and B_z are incorporated by means of a standard finite difference time domain approach [34, 35].

Next, PIC simulation results for the mentioned discharge scenario with a driving frequency $\omega =2\pi \times 13.56$ MHz and at a pressure of p = 1 mTorr shall be discussed. Figure 7 presents the correspondingly obtained spatial profile of the relative phase (solid red) between the transversal electric field and the current density (i.e., their parallel contribution) in comparison with the result from the analytical model (dashed black). Again shaded in gray is the interval $\pi /2\lt {\rm{\Delta }}\phi \lt 3\pi /2$ of negative power deposition. As stands out, a remarkable similarity in the results is observed, despite the very different level of complexity of the compared kinetic models. The spatial profile is qualitatively nearly equivalent. The observed physical consequences due to the anomalous skin effect are clearly inherent to both models. Interestingly, there are slight deviations in numbers, in particular as regards the position of phase transition from positive to negative power deposition at $x\approx 1.3\,\mathrm{cm}$. It is, however, important to acknowledge that the phase and therewith the exact position of the transition is a rather sensitive measure. This will be addressed shortly.

Zoom In Zoom Out Reset image size

The phase transition is intrinsically connected to the power deposition as depicted in figure 8 for both PIC (solid red) as well as the analytical model (absolute power density: dashed black; Ohmic power density: dashed–dotted blue). It is apparent that the position of phase transition toward negative power deposition has a one to one correspondence with the phase difference between the electric field and the current density. The transition to a negative real part in polar representation is reflected by a zero crossing of the absolute deposited power (i.e., a singularity in the semi-log plot shown). As expected from the phase difference, the singularity obtained from PIC is slightly shifted with respect to the analytical model results. Besides this shift, the principle physical dynamics as well as the quantitative dependence are resembled closely. In the frame of a comparison, two aspects seem quite peculiar:

Zoom In Zoom Out Reset image size

Firstly, the level of statistical noise observed in the PIC results. Being inherent to all Monte Carlo models, accounting for the essential stability and accuracy considerations, these are of rather cosmetic concern. In terms of an evaluation and interpretation, caution is suggested when the signal is too noisy as for instance in figure 8, for distances larger than $x\approx 3\,\mathrm{cm}$. The electric field strength is maximum at the left interface and strongly damped while penetrating the plasma. This is due to the low excitation frequency in comparison with the electron plasma frequency. Consequently, the signal greatly diminishes and thus the signal to noise ratio substantially increases from the plasma interface into the plasma (i.e., left to right).

Secondly, deviations due to contrasting model assumptions. While the goal was to differentiate the two models as least as possible, a few differences are intrinsic to the respective model. In approximate order of importance: (i) within the analytic model collisions are irrespective of energy taken into account by a constant collision frequency. In comparison to the collisional processes and energy dependent cross sections utilized in the PIC model, this is a rather crude approximation [36, 37]. (ii) In the analytic model the force term within the Boltzmann equation is reduced to the electric field force. In contrast, in the PIC model the complete Lorentz force term $\vec{F}=q(\vec{E}+\vec{v}\times \vec{B})$ is retained. (iii) The analytical model assumes the displacement current ${\varepsilon }_{0}\tfrac{\partial \vec{E}}{\partial t}$ to be negligible. It is, however, included in the full-wave PIC model. As the total current density within the plasma (the sum of the conduction and the displacement current density) is vastly dominated by the conduction current, this assumption seems quite unimportant. (iv) The analytic model is developed around a linear perturbation. Nonlinear effects are consequently excluded. Within PIC nonlinear effects are intrinsically included. However, for the investigations performed the field strength has been intentionally chosen small enough, so that nonlinear effects do not play any role [24].

To conclude with the comparison of the analytic and the PIC model, both appear to be in remarkable agreement given the detailed differences between the two model approaches. It is to notice, however, that a great difference between the two models lies within the time and computational effort required for the solution. While the analytical model is evaluated within seconds, typical runs for the (non optimized) PIC code took about two months to only provide a fair statistical basis. The analytical model therefore facilitates exceptional accuracy of the calculations results paired with a tremendously reduced computational burden.

Negative power deposition in low pressure radio frequency plasmas is investigated by means of an analytical model which couples Boltzmann's equation and the quasi-stationary Maxwell's equations. The formulated boundary value problem is of Sturm–Liouville type and is solved exploiting standard Hilbert space methods. An explicit solution for both, the electric field and the distribution function of the electrons has been found for the first time for a bounded unsymmetrical discharge configuration. Particularly, the anomalous skin effect with its peculiarities, e.g., the non-monotonic decay of the field distribution and the effect of phase mixing is discussed. The analytical solution is compared with results from particle in cell simulations. Here a one-dimensional electromagnetic full wave model has been developed and implemented. A comparison of the analytical and the numerical results show an excellent agreement.

The authors gratefully acknowledge the support by the Deutsche Forschungsgemeinschaft in the frame of Collaborative Research Centre TRR 87. Valuable discussions with and support from Prof. T Eisenbarth, Dr M Lapke, and Prof. R P Brinkmann are gratefully acknowledged.