# Nonlinear Poisson Solver

One of the governing equations for electrostatic plasma simulations is the Poisson’s equation,

$$\nabla^2\phi=-\dfrac{\rho}{\epsilon_0}=-\dfrac{e}{\epsilon_0}\left(Z_in_i-n_e\right)$$

In the types of discharges we typically consider here at PIC-C, the ion density is obtained from kinetic particles (the particle-in-cell method) while electrons are approximated as a fluid. The electron thermal velocity assuming 5eV electrons is over one million meters per second. On the other hand, typical ion velocities in plumes are around 15 km/s. This is a speed ratio of three orders of magnitude. The PIC method requires that particles travel less than a cell length per time step. As such, the simulation time step will be determined by the fastest species. Hence, if want to include electrons, we need to run the simulation at 1000 times smaller time step than what would be satisfactory for the ions. In other words, the simulation will take 1000 times longer to complete the same physical time span.

### Boltzmann relationship

Luckily, we can often get away with not simulating electrons as kinetic particles. This is true especially when magnetic fields are not present or we are dealing with motion in direction along the magnetic field lines. Consider the electron momentum equation,

$$mn\left[\dfrac{\partial \mathbf{v}}{\partial t} + \left(\mathbf{v}\cdot \nabla\right)\mathbf{v}\right]=qn\mathbf{E}-\nabla p$$

In the frame of reference of ions, electrons are seen to respond instantaneously to a disturbance, and thus the time dependent term can be neglected. The convective term can also be generally neglected since its magnitude is small (or zero if we consider spatially uniform velocity). This then leaves us with the force balance, \(q\mathbf{E} = \dfrac{1}{n}\nabla p\). We can next make three substitutions, \(q=-e\), \(\mathbf{E}=-\nabla \phi\), and \(p=nkT_e\) (assuming isothermal electrons), which gives us \(e\nabla\phi = \dfrac{kT_e}{n}\nabla n\). This equation can be easily integrated along some arbitrary path to obtain

\(e\left(\phi-\phi_0\right)=kT_e\ln\dfrac{n}{n_0}\). This very important relationship is known as the Boltzmann relationship. It is encountered frequently in plasma physics. The importance of the Boltzmann relationship lies in the fact that it relates local potential and the local electron density. Given a set of reference parameters, electron potential can be computed everywhere in the domain (where the governing assumptions hold) from simply knowing the electron density.

### Nonlinear Poisson’s Equations

The Boltzmann relationship can also be inverted to obtain an expression for electron density,

\(n_e = n_0\exp\left(\dfrac{\phi-\phi_0}{kT_e}\right)\)

This is the missing piece needed for our electrostatic potential solver! Instead of simulating electrons as particles, we can use this relationship directly in the Poisson’s equation,

$$\nabla^2\phi = -\dfrac{e}{\epsilon_0}\left[Z_i n_i – n_0\exp\left(\dfrac{\phi-\phi_0}{kT_e}\right)\right]$$

Unfortunately, this relationship results in a non-linear equation, since now the expression we are solving for (potential) appears on the right hand side as well. In matrix form, we have \(\mathbf{A}\phi=f(\phi)\). Luckily, solving non-linear equations is no rocket science – we just need little help from Mr. Newton.

### Newton’s Method for Systems

Newton’s method is a simple algorithm for finding function roots. Consider some arbitrary equation \(f(x)=0\). We can use Taylor’s expansion to compute the value of this function in the vicinity of a known point, \(f(x)=f(\bar{x}) + (x-\bar{x})f'(\bar{x}) + O(2)\). For small differences the second order terms drop out and we can write

$$f(x)=0\approx \bar{x}-\dfrac{f(\bar{x})}{f'(x)}$$

Numerically, the above equation is used to get consecutively closer approximations to the function root, `x_new = x_old - f(x_old)/f_prime(x_old)`. We can use a similar approach to treat a **system of equations**. The matrix analogue to the single equation is \(\mathbf{x}^{k+1} = \mathbf{x}^{k} – \mathbf{J}^{-1} \mathbf{F(x}^{k}\mathbf{)}\), where

$$

\mathbf{J}(x) = \dfrac{\partial f_i}{\partial x_j}(\mathbf{x}) = \left[\begin{array}{ccc} \dfrac{\partial f_1}{\partial x_1}(\mathbf{x}) & \dfrac{\partial f_1}{\partial x_2}(\mathbf{x}) & \ldots \\[0.3cm] \dfrac{\partial f_2}{\partial x_1}(\mathbf{x}) & \dfrac{\partial f_2}{\partial x_2}(\mathbf{x}) & \ldots \\[0.3cm] \vdots & \vdots & \ddots \end{array}\right]

$$

is the Jacobian. This algorithm is generally implemented in two steps two avoid the need to compute a matrix inverse. First, the \(\mathbf{J}(\mathbf{x})\mathbf{y}=\mathbf{F}(x)\) system is solved for \(\mathbf{y}\) using your favorite linear solver. The solution is then updated as \(\mathbf{x}^{k+1}=\mathbf{x}^{k}-\mathbf{y}\). The process continues until \(||y||<\epsilon_{tol}\), some prescribed tolerance.

### The non-linear Poisson solver algorithm

So how does this method look in practice when applied to the Poisson’s equation? Scroll down for the entire Matlab source code (note, the coding was actually done in Octave since we don’t have access to a Matlab license). To keep things simple, we apply the algorithm to a 1D system bounded on both ends with fixed potential nodes. The potential is -5V on the left boundary and +5V on the right boundary. Uniform plasma density of \(n_0=10^{16}\; \text{m}^{-3}\) and temperature \(kT_e=5\;\text{eV}\) is used.

Using the standard finite difference algorithm, we can write our governing equations as

$$\begin{array}{rcl}\phi_1&=&-5\\\phi_1-2\phi_2+\phi_3&=&Cn_{i,2} – C n_0\exp[(\phi_2-\phi_0)/kT_e]\\

\phi_2-2\phi_3+\phi_4&=&Cn_{i,3} – C n_0\exp[(\phi_3-\phi_0)/kT_e]\\

\vdots & = & \vdots\\

\phi_{n} &=& 5

\end{array}$$

Here the constant \(C=-e/\epsilon_0 \Delta^2h\). In other words, we have the system \(\mathbf{A}\mathbf{x} = \mathbf{b_0} + \mathbf{b_x}(\mathbf{x})\) (we are using \(\phi\equiv x\) for generality). The two vectors on the right hand side correspond to the constant and the non-linear term. As the first step in our algorithm, we compute the non-linear vector (the linear doesn’t change between iterations).

%1) compute bx for n=1:nn if (fixed_node(n)) bx(n)=0; else bx(n) = -C*den0*exp((x(n)-phi0)/kTe); end end

The example is made more general by using a `fixed_node` array to indicate which nodes are fixed and which are fluid. We compute electron density only on the fluid nodes. Next we add the two right hand side vectors and subtract to obtain \(\mathbf{F(x)} = \mathbf{A}\mathbf{x} – \mathbf{b} = 0\).

%2) b=b0+bx b = b0 + bx; %3) F=Ax-b F = A*x-b;

We next need to compute the Jacobian. Our function \(\mathbf{F(x)} = \mathbf{A}\mathbf{x} – \mathbf{b_0} – \mathbf{b_x}\) consists of three terms: the linear system given by the matrix multiplication, a vector which is not a function of `x`, and another vector which is. The derivative of the non-linear RHS vector is \(\partial (b_x)_i/\partial x_j = -\delta_{i,j}Cn_0\exp[(\phi_i-\phi_0)/kT_e]/kT_e\). The Jacobian is thus given by \(\mathbf{J}=\mathbf{A}-\partial (b_x)_i/\partial x_j=\mathbf{A}-\text{diag}(\mathbf{P})\). We first compute this `P` vector on fluid nodes,

%4) compute P=dbx(i)/dx(j), used in computing J for n=1:nn if (fixed_node(n)) P(n)=0; else P(n) = -C*den0*exp((x(n)-phi0)/kTe)/kTe; end end

and set the Jacobian,

%5) J=df(i)/dx(j), J=A-diag(dbx(i)/dx(j)) J = A - diag(P);

The rest is easy. We first solve for `y` and then use it to update the solution. The final profile is plotted in Figure 1.

%6) Solve Jy=F y = J\F; %7) update x x = x - y;

The iterations continue until convergence

%8) compute norm; l2 = norm(y); if (l2<1e-6) disp(sprintf("Converged in %d iterations with norm %g\n",it,l2)); break; end

### Complete Source Listing

Here is the complete source code in Matlab. The code was actually tested using Octave, so let us know if it doesn’t work in Matlab proper. Also, feel free to contact us if you need a Java implementation.

%Demo non-linear Poisson solver % we are solving: % d^2phi/(dh^2) = -rho/eps0*(ni-n0*exp((phi-phi0)/kTe)) % with phi=-5, 5 on left and right boundary % See https://www.particleincell.com/2012/nonlinear-poisson-solver/ % for additional info %clear screen clc %constants EPS0 = 8.85418782e-12; Q = 1.60217646e-19; %setup coefficients den0=1e16; phi0=0; kTe=5; %precomputed values lambda_D = sqrt(EPS0*kTe/(den0*Q)); %for kTe in eV dh = lambda_D; %cell spacing C = -Q/EPS0*dh*dh; %setup matrixes nn=15; %number of nodes A=zeros(nn,nn); fixed_node = zeros(nn,1); b0=zeros(nn,1); %left boundary A(1,1)=1; b0(1)=-5; fixed_node(1)=1; %right boundary A(nn,nn)=1; b0(nn)=5; fixed_node(nn)=1; %internal nodes for n=2:nn-1 %FD stencil A(n,n-1)=1; A(n,n)=-2; A(n,n+1)=1; fixed_node(n)=false; b0(n)=C*den0; end %initial values bx = zeros(nn,1); P = zeros(nn,1); J = zeros(nn,nn); x = zeros(nn,1); y = zeros(nn,1); %--- Newton Solver ---- for it=1:20 %1) compute bx for n=1:nn if (fixed_node(n)) bx(n)=0; else bx(n) = -C*den0*exp((x(n)-phi0)/kTe); end end %2) b=b0+bx b = b0 + bx; %3) F=Ax-b F = A*x-b; %4) compute P=dbx(i)/dx(j), used in computing J for n=1:nn if (fixed_node(n)) P(n)=0; else P(n) = -C*den0*exp((x(n)-phi0)/kTe)/kTe; end end %5) J=df(i)/dx(j), J=A-diag(dbx(i)/dx(j)) J = A - diag(P); %6) Solve Jy=F y = J\F; %7) update x x = x - y; %8) compute norm; l2 = norm(y); if (l2<1e-6) disp(sprintf("Converged in %d iterations with norm %g\n",it,l2)); break; end end disp(x'); plot(x);

You may also be interested in the article on solving potential n composite dielectrics.

Code Optimization: Speed up your code by rearranging data access

Particle In Cell Method in Cylindrical Coordinates

Multiscale Modeling of Hall Thrusters

Thanks so much for posting this! It was extremely helpful, and is working perfectly. Thanks – Isaac

Great post, very helpful, thanks!

I will just offer one quick (constructive) criticism. The line “Here the constant C=-e/\epsilon_0 \Delta^2h” is confusing. Perhaps it’s just my misinterpretation of the nomenclature, but \Delta^2h looks like “delta squared times h”, not “delta h squared”. The meaning becomes clear in the code, but the text is confusing.

Thanks Glen! This has been the standard syntax through my various numerical methods courses here in the USA. Perhaps things are written differently in Australia (judging by your email address)? Basically \(\Delta^2 h \equiv (\Delta h)^2\).

Thanks for that. I wondered if I had worded what I said a bit strongly. I’ve just never seen it written that way before and the document is obviously so close to perfect I thought I’d ask. Should have used a question mark, there’s my mistake. ðŸ™‚

This is the first document of yours that I have read, I had some difficulty working out what was what, but I’ll get better, and I look forward to reading more. Thanks for the reply.

No worries, that was a completely valid question. One of the great things about blogging is that it is a great two way learning curve. Before your comment it didn’t even occur to me that was not a universal symbol, now I know… So thank you!

Well I checked with an Aussie who works in the field and he said he’d never seen it done that way, but right now I’m reading an American text book online and I just saw that notation for the first time. Maybe there is a difference? Or maybe we read different textbooks? ðŸ™‚

Also I wanted to say that this page has been a constant reference for the past month or so and it is one of the best pieces of explication I’ve read. I’ve not done anything in modelling before but this is right at the level I need and on the topic I need. Actually the only reason I commented the first time was that the symbol(s) in question was (I thought) the only part of the theory that was in any way ambiguous, the rest of it was all totally clear, and it seemed strange and trivial. So thanks again.

Thanks for a very helpful post! Is there a 2D matlab implementation too?

Thanks a lot for this useful tutorial.. Very well-written and clear. I was wondering if you could help me with solving the same nonlinear Poisson equation using the finite volume (i.e. box integration) method.

When you write: “The Boltzmann relationship can also be inverted to obtain an expression for electron density”

n = n_o * exp( (phi-phi_o)/(kT) )

why is the charge “e” disappeared?

Shouldn’t it be:

n = n_o * exp( e(phi-phi_o)/(kT) ) ?

Thank you.

His kTe is in electron volts so “e” is cancelled in the numerator.

Thanks for the excellent example! The bx term relating to the electron charge (exp[(phi-phi0)/kTe]) can get pretty large when dealing with high voltages and low temperatures – even running into overflow limits (I should note that I’m not working in octave or matlab, so I don’t know if they have special ways of dealing with very large numbers.) Do you know if there is a good way of dealing with this issue? Perhaps it is possible to treat the system in log space?

You can use Taylor Series to rewrite \(e^x = \sum_{n=0}^{\infty}\frac{x^n}{n!}=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\cdots\)