In plasma theory, electron density distribution is often expressed by the Boltzmann equation as

where n₀ is a constant density, e is the charge of an electron, V(z) is an electric potential distribution, k is the Boltzmann constant, and T_e is the electron temperature. V_p is the electric potential at a particular position (usually plasma potential). This Boltzmann relation between charged particles and electric potential is not only based on statistical theory but also on the balance between electric force and pressure gradient force [1].

If the electron density is expressed as Eq. (1) and an ion density is given then the Poisson equation becomes the following form of the Poisson-Boltzmann equation [2]

where n(z) is a spatial distribution of ion density and ϵ₀ is the permittivity of vacuum. This kind of the Poisson-Boltzmann equation in which representing charged particle density as a Boltzmann relation is important not only in plasma but also in many other scientific and engineering fields [3]. As in the case of hot plasma, there are many cases that also the ion density (hole in semiconductors) is expressed by the Boltzmann relation. However, ions are cold in a processing plasma so that Boltzmann relation usually can not be a good approximation. Ion density will be given as numerical values or in the form of another function of the potential.

The Poisson-Boltzmann equation is a nonlinear differential equation, and the situations and conditions that can obtain an analytical solution are limited [4]. Therefore, there are many cases where the solution is obtained by transforming it into a numerical approximation formula such as a finite difference form. In particular, an analytical solution cannot be expected when the density of ions is given as only numerical values. The Poisson-Boltzmann equation, transformed into a numerical approximation, becomes a matrix equation, and the solution is usually obtained through an iteration method such as the Newton-Raphson method [5,6]. However, as is well known, the Newton-Raphson method has several failure factors such as overshoot, stational point, fore initial estimate, and so on. Of course, programs exist and used as a result of solving such problems, but research and development are continuing to find new methods.

This research proposes a method to obtain the solution of the Poisson- Boltzmann equation using the Lambert W function [7]. The most important advantage of the current method is that it is simple and easy to apply.

In the subsequent section, details of the new method will be explained and the results applied to several situations will be analyzed and discussed.

The Poisson-Boltzmann equation, Eq. (2), can be normalized by thermal potential V_t and the Debye length λ as

And here,

In order to deal with Eq. (3) numerically, an approximated equation will be obtained through various numerical methods such as the finite difference method (FDM), the finite element method, the finite volume method, and so on. The present method is not applicable only to the FDM. However, it will be explained based on the FDM for simplicity.

The finite difference form of Eq. (3) becomes

where Δx is the size of the spatial grid space and ϕ_j, n_j are the j-th element of ϕ, n, respectively. In the explicit scheme, Eq. (5) will be given as

where ${\varphi }_j^i$ is the i-th iteration step value of ϕ_j. But it is well known that the explicit scheme is hard to be numerically stable. Therefore, the following fully implicit scheme can be attempted.

As the next step, the Newton-Raphson method can be adopted for linearization of iteration process. And then, a matrix inversion or a method for linear systems will be used to obtain the solution of Eq. (7).

In this study, we develop a semi-implicit method that assumes only the center value as the implicit term as

The values in the next step are then determined through the Lambert function W(x) as

where

The current semi-implicit Lambert function method is very simple and easy to deal with because it does not require a matrix inversion or a numerical solution of linear system. Furthermore, the linearization approximation, which is used in the Newton-Raphson method, is not necessary.

The initial values of ϕ is determined on the assumption that the electron density equals the ion density as

Also, as can be seen in Eq. (9), an overflow may occur due to an argument in which the Lambert function contains an exponential function term. In that case, it can be solved through the following asymptotic form of the Lambert function.

where L₁ = ln z and L₂ = ln ln z [8]. In other words, it is not necessary to calculate the e^D in Eq. (9) directly because W(z) is calculated as a function of ln z.

Three cases were assumed to investigate the applicability of the current method: an uniform ion density, the collisionless sheath of cold ion plasma and a nonuniform ion density.

3.1. Uniform ion density

Figure 1 shows the results of assuming a spatially uniform ion density. The assumed conditions in this calculation are as follows:

Figure 1. The calculated ϕ(z) and normalized densities in the case of uniform ion density.

If n₀ is 10¹⁶ m⁻³, the Debye length will be about 1.7 mm.

Analyzing the numerical stability conditions of the current method in a general range is not easy because of the term of the Lambert function of the exponential function which can be seen in Eq. (9). However, the solution was stably obtained under wide range of condition variables.

The boundary condition inside the plasma body, that is, the boundary condition at z = 0, may be replaced by the condition that the electric field is 0. In both cases, the results were the same. If one want to deal with a floated plasma rather than this grounded plasma, ϕ will be replaced with the value added with the plasma potential V_p. Due to the substitution in Eq. (4), V_p does not appear in this mathematical modeling. V_p is determined by an additional physical condition and it becomes the plasma potential in the current case.

Even when the electrode is positively biased, mathematical solutions can be obtained stably. However, in such a case, the electron flux and electron density on the electrode surface is limited by thermal diffusion.

3.2. Collisionless sheath of cold ion plasma

Within a collisionless sheath of cold ion plasma, uniform flux and ion energy conservation can be assumed [9]. Threfore, the density of ions can be expressed as a function of the electric potential as

where E₀(= T_e/2) is the ion energy at entrance of the sheath. Thus, the Poisson equation of Eq. (2) becomes the following Bohm’s sheath equation form [10].

In this case, the ion density is given as another function of the potential rather than fixed values. Since the ion density can not be treated implicitly, two iterations are required now. If ϕ is determined by the first iteration, the next iterative processes can be attempted based on the ion density through that value. However, it is also possible to unify the two iterations into one. Therefore, the finite difference form, Eq. (8) becomes

Now, D in Eq. (9) is replaced by the following form.

Figure 2 shows ϕ and the normalized densities approaching a convergence and Fig. 3 is the converged densities. It shows the characteristics of a typical sheath with a reduced ion density and a very low electron density. The assumed conditions are the same as in the case of uniform ion density. The same is also true about the plasma potential V_p and the electrode.

Figure 2. ϕ(z) and normalrized densities approaching a convergence in the Bohm sheath case.

Figure 3. The converged normalized densities in the Bohm sheath case.

As such, the present method is not only possible if the ion density is given as some fixed values. It is possible also when the ion density is a function of electric potential. If the nonlinear dependence on the potential of ion density is not severe compared to the exponential function term of electron density, it can be treated as an explicit term. If the dependence of ion density on the potential is linear, it can be included in the process of obtaining the Lambert function result through variable substitution.

3.3. Non-uniform ion density

As discussed in the introduction, the Boltzmann relation of electron density means the force balance between the electric force and the pressure gradient force. So, in the case of uniform ion density, when both boundaries are grounded, ϕ becomes zero at all points. In order to generate a potential difference (or an electric field) when both boundaries are grounded, an ion density gradient must exist. Therefore, to examine such an ion density gradient effect, the following form of non-uniform ion density is assumed.

Figure 4 shows the results when different potentials were applied to the two electrodes. It can be seen that asymmetric electron densities are formed in the sheath regions. In the plasma body, ϕ becomes 0 (that is, V = V_p). V_p is not determined only by the Poisson Boltzmann equation and boundary conditions for ϕ, but additional conditions are required. We can see the result in Fig. 5 when ϕ is zero at both boundaries, in which case the electron density the ion density are nearly the same at all points and no prominent sheath structure is seen. That is, V = V_p at both boundaries, and the formed potential is just the presheath potential generated by the balance between the electric force and the pressure gradient force.

Figure 4. The results when ϕ(0) = −50, ϕ(L) = −100 in the case of nonuniform ion density.

Figure 5. The results when ϕ(0) = ϕ(L) = 0 in the case of nonuniform ion density.

The floated sheath structure by an actually grounded electrode surface will be implemented under the condition that ϕ = −V_p (that is, V = 0). Figure 6 shows the sheath structure formed under such conditions. The value of V_p can be determined by the electron density at the boundary. Therefore, the value of V_p can be assumed based on the following Bohm sheath potential V_s.

Figure 6. The results when ϕ(0) = ϕ(L) = −V_p in the case of nonuniform ion density.

where m and M are the masses of electron and ion, respectively.

A numerical method using the Lambert W function to solve the nonlinear Poisson-Boltzmann equation was developed. The main focus of the present method is to put the exponential function term, which is the main cause of numerical instability, as an implicit value, and the remaining terms as explicit values to avoid complexity. As the result of that, a very simple scheme without matrix inversion is obtained. It was shown that solutions can be obtained even if the density of ions is given as a function of potential. And in that case, it was found that the iterative steps could be integrated and further simplified.

In terms of the calculation speed and stability, the present method will not always be superior to other methods. In many cases, Newton’s method can ensure sufficient stability and there are methods of treating a linear system of size m with the calculation time of the order m. However, also, for certain problems, the current method may be efficient.

The present method will not only be applicable to one-dimensional problems, nor will it be based solely on the finite difference method.

Also, the Poisson-Boltzmann equation, which includes Boltzmann relation for charged particles, is not limited to cold ion plasmas, but is also used in other research areas. According to a survey by Gray and Stiles [3], there are more than 20 such fields in science and engineering. The current method is also expected to be applicable to certain problems in such areas.

This work was supported by a funding for the academic research program of Chungbuk National University in 2022.

The authors declare no conflicts of interest.