The Poisson Solver in SIMION (Version 8.1 Experimental Feature)

DRAFT-20080228 (D.Manura)

Overview

SIMION's electric field solver (Refine function) has recently been extended to support solving the Poisson Equation. As before, the Refine function can calculate potentials inside a volume defined by electrodes, but the new Poisson solving capability allows it to also handle the case when a known, arbitrary distribution of space-charge density also fills that volume (such as due to high current charged particle beams or charged insulators)

Applications

Taking into account space-charge density can be important for certain problems. Space-charge effects tends to be most severe for a slowly moving, high current beam confined in a small volume, for this causes a large space-charge density and large periods of time over which the fields from that space-charge density can have an effect on charged particle motion. A typical example is the beam near the surface of a thermionic cathode operated at maximum (space-charge-limited) current, for here the beam is not just high current but also at low energy (thermal energies as not yet having been accelerated by the anode).

The Poisson solver can also be useful for other situations in which a charge density distribution is known a priori, such as in the case of a charged insulator material with known charge distribution. This is simpler to handle because it does not require particle flying to do the Refine.

Version Requirements

The Poisson solving capability is one of the major new features that will be formally included in 8.1. It is already partially implemented (currently for 2D planar and 2D cylindrical arrays but not yet for 3D arrays -- FIX-TODO). Prior to the 8.1 release, SIMION 8.0 users can preview this feature by downloading SIMION version 8.0.5-TEST5 or above and putting SIMION in "early access" mode. This is noted in the following issues reports:

http://simion.com/issue/482 - Refine extended to solve the Poisson Equation
http://simion.com/issue/421 - Early Access mode for accessing 8.1 (or 9.0?) features ).

Previous versions of SIMION (8.0 and below) assumed that space-charge density was negligible. That is, it solved only the Laplace equation. However, 7.0 and 8.0 do have some rough "charge repulsion" methods that are useful in approximating space-charge effects under certain conditions (see the repulsion example in SIMION 8.0.4 or above for full details on this).

Math

The Poisson equation is formally

$\nabla^2 \Phi = -\rho / \epsilon_0$

where $\phi(x,y,z)$ is electrostatic potential as a function of position $(x,y,z)$ , $\rho(x,y,z)$ is space-charge density as a function of position, and $\epsilon_0$ is the permittivity of free space constant. To solve for the potentials ( $\phi$ ) inside some volume (i.e. Refine function), SIMION must in general know not only $\phi$ on the boundary of the volume but also $\rho$ throughout the inside of that volume.

In the case of no space charge ( $\rho = 0$ ), the equation reduces to the Laplace equation, which earlier versions of SIMION solved.

SIMION solves the Poisson equation with a finite difference scheme. The equations are similar to those in the SIMION 7.0/8.0 but with the addition of a charge term. For example, for 2D planar and 2D cylindrical systems, the finite difference equations are respectively

$P0 = \frac{1}{4}(P1 + P2 + P3 + P4 + (\bar{\rho}/\epsilon_0) h^2)$

$P0 = \frac{1}{4}(P1 + P2 + (1-1/2r)P3 + (1+1/2r)P4 + (\bar{\rho}/\epsilon_0) h^2)$

where $\bar{\rho}$ is the average space-charge density inside the square of area $h^2$ centered at P0. The formulation is based mathematically on a finite element assumption that charge density is uniform inside each of the four squares having a vertex at P0. Except in special cases, this assumption is only met approximately, but it improves as the grid spacing h approaches zero.

SIMION Array Modeling

To use the Poisson solver, you will need to create an electrostatic PA, as usual, plus an additional PA in which to hold charge densities. This "charge density PA" has the same format as electrostatic PAs, but we are storing charge densities, rather than potentials, in it.

Currently, if the size of the electrostatic PA is (wx+1, wy+1) grid points, then the size of the charge density PA must be $(wx/2^N,
wy/2^N)$ grid points for non-negative integer N and dimensions. N > 0 allows you to optionally have a charge density PA at a lower mesh resolution than the electrostatic PA, which might be useful in cases when the charge density is less critical or less known. The absence of the "+1" in the charge density array size occurs because the cells in the charge density PA are interpreted differently than in an electrostatic PA. In an electrostatic PA, points (0,0), (1,1), (1,0), and (0,1) represent the vertices of the lower-left volume cell: if all four points are marked as electrode points, then the cell volume between them is interpreted as a solid electrode. In the charge density PA, point (0,0) represents the entire volume of the lower-left cell itself. This distinction is made because in the charge-density PA, the volume in a cell itself holds charge density, while in an electrostatic PA we are concerned with representing boundaries (surfaces) of cells on which to apply potentials (including zero width ideal grid boundaries).

Note: When SIMION's solver is also extended to support dielectrics (or other material properties), an additional "dielectrics PA" of size $(wx/2^M, wy/2^M)$ , integer M >= 0, would be used to hold relative dielectric constants, and point (0,0) will also represent the lower-left cell itself holding dielectric material.

The symmetry of your electrostatic array in general must match the symmetry of your charge density array. If you have a geometry with 2D cylindrical symmetry and 3D asymmetrical space-charge inside of it, then the potentials inside are in general 3D asymmetrical and cannot be modeled with simply a 2D array.

For 2D arrays, the charge density values (Coulombs per volume) stored in the space-charge array must be specified in units of C gu^-2 mm^-1 (Coulombs per square grid unit per mm). For 2D planar arrays, each grid cell is interpreted as being extruded (projected) infinitely in the Z direction, and this value represents the Coulombs per grid cell for each mm of that extrusion. For 2D cylindrical arrays, each grid cell is interpreted as being rotated (projected) around the X axis, and this value represents the Coulombs per grid cell for each mm of that projection around the axis (note: the length of that projection is $2 \pi r$ for a grid cell centered at radius r, so the total charge represented by a grid cell is $2 \pi r$ times the charge density value). You may alternately think of these units as C gu^-3 (mm gu^-1)^-1, where (mm gu^-1) is the mm-per-grid-unit scaling factored defined on the PAs tab in View. Embedding this mm-per-grid-unit scaling constant into the units of charge density for the values in the space-charge array avoids having to specify the mm-per-grid-unit scaling during the Refine process. Note also that for planar (2D or 3D) arrays with a mirror plane, the actual charge represented by a grid cell is twice that charge value stored in the grid cell since the array is projected across the mirror plane. It is critical to understand how these units are defined so that charge density is not accidentally scaled by some factor.

Usage

Currently, there is no user interface for using the Poisson solver. The Poisson solver must be invoked from Lua code. See the Poisson examples for examples on how to do this. To use the Poisson solver, SIMION must be switched into the "early access" mode (http://simion.com/issue/421). One must create the charge density PA object and fill its cells appropriately with charge densities (rather than potentials). Then the electrostatic PA object is refined while passing it the space-charge array as an option.

Computing charge density from trajectories

Typically, in practice, the charge density distribution is not known a priori but is rather calculated from particle trajectories. Often, the particle trajectories and fields are mutually dependent and so are solved alternately, in iteration. For a continuous beam, you might solve the fields and trajectories alternately until a mutual convergence is reached. For a pulsed beam, you might re-refine the field on every time-step, or some multiple of time-steps (e.g. whenever the particles advance by one grid unit).

Each particle traced shall represent some fraction of either the beam current (for continuous beams) or beam charge (for pulsed beams). For example, a 1 mA continuous beam traced with 100 particles, gives 0.01 mA per particle. Alternately, a 1 mA (i.e. 10^-3 C/s) beam pulsed for 1 millisecond represents 10^-6 Coulombs. If the pulse is modeled with 100 particles throughout the volume of the bunch, then each particle represents 0.01*10^-6 Coulombs.

Additional care is needed when the beam is modeled with symmetry (e.g. 2D planar, 2D cylindrical, or mirroring). For a 2D planar beam (e.g. the 2D planar Pierce gun), where the 2D image of the trajectories is interpreted to be extruded (projected) infinitely in Z, each particle represents some current per distance of the extrusion (i.e. mA/mm) or charge per distance (C/mm). For a 2D cylindrical beam (e.g. the 2D cylindrical Pierce gun), where the 2D image of the trajectories is interpreted to be rotated (projected) around the X axis, each particle represents some current per distance of the projection (i.e. mA/mm) or charge per distance (C/mm). Similarly, mirror planes cause a beam to represent twice its current or charge otherwise.

2D cylindrical beams can require special care since the proportion of charge that a given cell in the charge density array represents is proportional to the distance of the center of the cell from the X axis. The 2D cylindrical Pierce gun example is a good example of how to handle this. If a traced particle that starts at a large distance from the axis is a assigned a certain mA or C (mA/mm or C/mm) but later moves closer to the X axis (e.g. beam narrows), then the mA or C is conserved but the mA/mm or C/mm increases since the charge is now distributed around a smaller circle.

Some custom Lua code is required to deposit the space-charges from the time-step of each trajectory integration into the charge density array. This can occur in the other_actions segment of a Lua workbench program. A refine operation on the PA must be triggered at appropriate stages of the operation, such as between runs in a terminate or initialize segment (continuous beam) or between time-steps in an other_actions segment (pulsed beam). Also, reruns might be scheduled for the iterations, i.e. sim_rerun_flym=1 in a terminate segment (in a continuous beam). Over time, some library code and GUI interfaces will likely be further developed to make parts of this easier. The pierce2dplanar and pierce3dcylindrical examples both makes use of a picutil.lua module with some functions that assist in this. For example, there is code to deposit space-charge from the particles trajectories into the charge density array cells, as well as to perform smoothing and display of the charge density array.

Performance

The performance of the Poisson solver will naturally likely be slower than the regular Laplace solver, though not substantially so. The charge density array also uses addition memory. The Pierce gun examples with a non-continuous beam have run times on the order of a few minutes. Using a lower density space-charge array might speed things up, or at least use less memory.

Examples

See the SIMION Poisson Solving Examples for examples on using the Poisson solver. Note that the Poisson solver is an advanced feature. It currently requires Lua user programming to effectively make use of it, but it can be very powerful and flexible. It is recommended that you study the source code.

Also, the Gauss' Law example (example/gauss_law) may be useful when experimenting with the Poisson solver.

Current Limitations / TODO

The Poisson solving is currently only known to work on basic arrays (.PA files). Using it with fast adjust arrays (.PA# files) or with time-dependent voltages has not yet been studied.

Currently the Poisson solving is only implemented for 2D planar arrays. Support for 3D arrays will added next (FIX - TODO).

The GUI and programming interfaces to the Poisson solver are currently limited but functional.