## Introduction

In 1961, I was awarded a research grant by the DSIR (Department of Scientific and Industrial Research) to carry out work on electrical discharges for a PhD degree, under the supervision of Professor P M Davidson. It will be useful to have a short background of this work to see the role that computing was going to play. Readers who are not interested in these details can jump to the point where the hardware is introduced.

## Growth of ionisation in gases

### Steady state equations

Gases are normally good insulators - a property that is of huge importance in the power distribution network. If the electric field in a gas is sufficiently large, however, then a rapid transition to a conducting state can occur (an electric spark) and currents of thousands of amps can flow, usually limited only by the ability of the power supply and the connecting leads to maintain the current.

In the simplest case, we can assume that the current in the gas is carried by free electrons and positive ions. One electron, travelling through an electric field, can gain enough energy to ionise neutral atoms, so producing an additional electron and a positive ion. This mechanism is called the primary ionisation process (although in the treatment of detectors like Geiger-Müller tubes, it is sometimes called secondary ionisation, to contrast it with the ionisation produced by the incident high-energy particle). The primary process is described mathematically by the equation

*n*(*x*) = *n*(0)e^{ax}, (1)

where *n*(*x*) is the number of electrons reaching a distance *x*, *n*(0) is the initial number of electrons at position 0, and a is the primary ionisation coefficient, sometimes called the Townsend primary coefficient in honour of J S Townsend who developed this theory in about 1903^{1}. [Superscript numbers refer to the publication list at the end of this document.] A group of electrons arising from one initiating electron is called an avalanche. (We note that, if electrons are taken to travel in the positive *x* direction, then the electric-field vector points in the negative *x* direction, because of the negative charge of the electron. It is convenient, however, to consider the electrons to move in the positive *x* direction, and this should not cause any problem.) The exponential growth of electron numbers described in equ. (1) cannot continue indefinitely and in a well-controlled experiment the growth is terminated when the electrons reach the positive electrode (anode) at *x* = *d*. Similarly the positive ions generated by the primary ionisation process travel in the negative *x* direction, and end up at the negative electrode (cathode) at *x* = 0.

If only the primary process was active, then the current would stop flowing as soon as all the charged particles were collected at the electrodes. Townsend introduced several secondary processes to explain why the current can, in practice, continue to flow. The simplest one describes how new electrons can be generated at the cathode by the impact of positive ions there, and the secondary coefficient g is defined as the probability that a new electron is produced following the arrival of each positive ion. Putting these two processes together, Townsend derived the equation

*I*(*d*) = *I*_{0}e^{ad}/[1 - g(e^{ad} - 1)], (2)

where *I*(*d*) is the current arriving at the anode at a distance *d* from the cathode, while *I*_{0} is a small initiating current supplied at the cathode by some external means, e.g. by irradiating it with ultra-violet light. (If this initiating current was not supplied, then we would have to wait for the random appearance of an electron, from cosmic rays for example.) Equ. (2) has been the basis of many experiments. Keeping *I*_{0} constant, *I*(*d*) is measured as a function of *d* (i.e. for different anode-cathode separations), while the voltage is increased in proportion to *d*, so as to keep the electric field, and therefore a, constant. A typical graph of *I*(*d*) against *d* is shown in Fig. 1. The initial part of the graph is exponential, and therefore appears as a straight line when a logarithmic scale is used on the *y*-axis. For large *d*, however, the graph becomes increasingly steep, and tends to infinity when the denominator of equ. (2) becomes zero, *i.e.* at the value of *d* = *d _{b}* given by

g(e^{ad} - 1) = 1. (3)

This is known as the breakdown criterion, and describes a transition to a regime that is not accessible to steady state experimental measurement.

**Fig. 1.** Evaluation of equ. (2) for the parameters *I*_{0} = 1, a = 5, g = 0.0065

## Temporal growth

Up to this point in the theory, the mathematics does not appear sufficiently complicated to require the assistance of a computer. There is some scope for computation in applying non-linear curve-fitting to data described by equ. (2), so that the experiments yield values of a and g, together with estimates of their uncertainty (errors). This is an interesting topic in its own right, but is not the main computational project that we are interested in. We ask what happens when the electrode separation is greater than *d _{b}*, or when the voltage is increased until it is just greater than the value needed to satisfy equ. (3). In this case, a steady state cannot exist. As soon as the initiatory current

*I*

_{0}is switched on (or as soon as the voltage is applied between the electrodes), the current grows extremely rapidly, eventually exceeding the capacity of the power supply to keep the voltage constant. To study this process (temporal growth of ionisation), we need to introduce time into the equations, together with the average velocities of the electrons,

*w*, and positive ions,

_{e}*w*. The equations basically describe the steady motion of particles in the positive and negative directions, together with the primary ionisation process, which increases the numbers of both types of particles. Representing the number densities of electrons and ions by

_{p}*n*(

_{e}*x*,

*t*) and

*n*(

_{p}*x*,

*t*), the equations take the form

(¶/¶*t*)*n _{e}* = - (¶/¶

*x*)(

*w*) + a

_{e}n_{e}*w*, (4)

_{e}n_{e}(¶/¶*t*)*n _{p}* = (¶/¶

*x*)(

*w*) + a

_{p}n_{p}*w*. (5)

_{e}n_{e}In addition, we have a boundary condition at the anode:

*n _{p}*(

*d*,

*t*) = 0, (6)

indicating that no positive ions are emitted from the anode, while at the cathode we have, for example

*w _{e}n_{e}*(0,

*t*) = g

*w*(0,

_{p}n_{p}*t*), (7)

(or a similar equation representing secondary emission due to the arrival of photons at the cathode). There is also a need for some initiatory mechanism, which could be as little as a single electron, otherwise equ. (4) - (7) have the trivial solution *n _{e}* º

*n*º 0.

_{p}Equ. (4) and (5) show a partial symmetry between the electrons and the positive ions, but the ionisation term is the same in both equations, which breaks the symmetry and makes the solution more difficult. Various attempts^{2,3,4} were made in the 1930s and 1940s to solve these equations. We can obtain a rough estimate as follows^{5}: Equ. (3) can be regarded as a replacement condition. If one "primary" electron leaves the cathode, then m, defined by the left-hand side of equ. (3) gives the average number of secondary electrons produced by the resulting avalanche. Clearly, if m > 1, subsequent generations of electrons will increase in the ratio 1 : m : m^{2}, … The time between generations, t, is of order *d*/*w _{e}*, when the secondary emission is due to photons, and

*d*/

*w*when it is due to positive ions. (It is not difficult to obtain more accurate expressions, but these are near enough for the present discussion.) At a time

_{p}*t*, the number of generations will be

*t*/t, and the current will be proportional to m

^{t/t}= e

^{lt}where l = (lnm)/t. This is an "exponential" increase (with the correct meaning of the word!). Inserting typical values of the parameters into this formula, we can easily see how a significant increase in current can occur in times of order µs to ms.

The definitive work in solving equ. (4) - (7) was done by P M Davidson in a series of papers beginning in 1947^{6,7,8,9,10}. He used the method of Laplace transforms, and obtained the solution in the form of a contour integral. This could be evaluated by summing the residues at all the singularities of a complex function, leading to an infinite series of exponential terms. Evaluating the solution numerically is a significant problem and needs the help of a computer, but this has generally not been done. The reason is that the exponential terms contain one real, increasing term, and all the others are damped relative to this one. In the experiments being done at the time, it was not possible to observe small, rapidly changing, currents. By the time the current had increased to a measurable value, the solution was described to a good approximation by the real exponential term alone, and the evaluation of this is a much simpler problem.