• Home
• Sitemap
Article View

## Research Paper

Applied Science and Convergence Technology 2023; 32(3): 69-72

Published online May 30, 2023

https://doi.org/10.5757/ASCT.2023.32.3.69

## Numerical Method for Poisson-Boltzmann Equation Using the Lambert W Function

Nam-Sik Yoon*

School of Electrical Engineering, Chungbuk National University, Chungbuk 28644, Republic of Korea

Correspondence to:nsyoon@cbnu.ac.kr

Received: January 28, 2023; Revised: April 14, 2023; Accepted: April 15, 2023

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc-nd/4.0/) which permits non-commercial use, distribution and reproduction in any medium without alteration, provided that the original work is properly cited.

The Poisson-Boltzmann equation is a well-known nonlinear differential equation that is fundamental in plasma theory. However, since it has nonlinear characteristics, the conditions and situations for obtaining analytical solutions are limited. Therefore, research and development on ways to obtain solution through numerical analysis methods have been continued. In this paper, a numerical method using the Lambert W function to solve the nonlinear Poisson-Boltzmann equation is explained. To investigate the applicability of the current method, three cases were assumed: a uniform ion density, the collisionless sheath of cold ion plasma and a nonuniform ion density. Solutions can also obtained even if the density of the ions is given as a function of the potential. In this case, the iterative steps can be integrated and simplified further.

Keywords: Poisson-Boltzmann equation, Plasma, Numerical analysis, Lambert function

### 1. Introduction

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

$ne(z)=n0expeV(z)−eVp kTe ,$

where n0 is a constant density, e is the charge of an electron, V(z) is an electric potential distribution, k is the Boltzmann constant, and Te is the electron temperature. Vp 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]

$d2V(z)dz2=en0ε0expeV(z)−eVpkTe −n(z),$

where n(z) is a spatial distribution of ion density and ϵ0 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.

### 2. Lambert function method

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

$d2ϕ(x)dx2=eϕ(x)=n(x).$

And here,

$ϕ=V−VpVt,Vt=kTee,x=zλ,λ=ε0Ten0e2.$

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

$ϕj+1−2ϕj+ϕj−1=(Δx)2[eϕj−nj],$

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

$ϕj+1i+1−2ϕji+1+ϕj−1i+1=(Δx)2[eϕji−nji].$

where $ϕji$ 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.

$ϕj+1i+1−2ϕji+1+ϕj−1i+1=(Δx)2[eϕji+1−nji].$

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

$ϕj+1i−2ϕji+1+ϕj−1i=(Δx)2[eϕji+1−nji].$

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

$ϕji+1=D−W(Δx)22eD,$

where

$D≡ϕj+1i+ϕj−1i2+(Δx)22nji.$

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

$ϕj=lnnj.$

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.

$W(z)=L1−L2+L2L1+L2(−2+L2)2L12+L2(6−9L2+2L22)6L13+L2(−12+36L2−22L22+3L23)12L14𝒪(L2 L1 )5,$

where L1 = ln z and L2 = ln ln z [8]. In other words, it is not necessary to calculate the eD in Eq. (9) directly because W(z) is calculated as a function of ln z.

### 3. Calculation and analysis

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.

$L=0.2m, Te=5eV, ϕ(0)=0, ϕ(L)=−100.$

If n0 is 1016 m−3, 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 Vp. Due to the substitution in Eq. (4), Vp does not appear in this mathematical modeling. Vp 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

$ni(z)=n01−eV(z)E0 −1/2,$

where E0(= Te/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].

$d2V(z)dz2=en0ε0expeV(z)−eVpkTe −1−eV(z)E0 −1/2.$

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

$ϕj+1i−2ϕji+1+ϕj−1i=(Δx)2eϕji+1 −11−2ϕji.$

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

$D≡ϕj+1i+ϕj−1i2+(Δx)2211−2ϕji.$

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 Vp 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.

$n(z)=sinπzL+1.$

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 = Vp). Vp 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 = Vp 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 ϕ = −Vp (that is, V = 0). Figure 6 shows the sheath structure formed under such conditions. The value of Vp can be determined by the electron density at the boundary. Therefore, the value of Vp can be assumed based on the following Bohm sheath potential Vs.

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

$Vs=Te2elnM2πm,$

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

### 4. Conclusions

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.

### Acknowledgments

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

### Conflicts of Interest

The authors declare no conflicts of interest.

1. F. F. Chen, Introduction to Plasma Physics and Controlled Fusion (Plenum Press, 1984). p. 75.
2. K. Köhn, D. Krüger, E. Kemanceci, L. Xu, E. Eremin, and R. P. Brinkmann, Plasma Sources Sci. Technol. 30, 105014 (2021).
3. C. G. Gray and P. J. Stiles, Eur. J. Phys. 39, 053002 (2018).
4. L. G. D'yachkov, Tech. Phys. Lett. 31, 204 (2005).
5. D. Kim and D. J. Economou, J. Appl. Phys. 94, 2852 (2003).
6. Newton’s method https://en.wikipedia.org/wiki/Newton’sMethod (accessed May 20, 2023).
7. Lambert W function https://en.wikipedia.org/wiki/LambertWfunction (accessed May 20, 2023).
8. R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth, Adv. Comput. Math. 5, 329 (1996).
9. N. Sternberg and V. Godyak, IEEE Trans. Plasma Sci. 35, 1341 (2007).
10. D. Bohm, The Characteristics of Electric Discharges in Magnetic Fields (McGraw Hill, 1949).