Numerically solving for the floating potential

In this example we consider how to numerically obtain the floating potential by means of a root solver. We shall consider a sphere situated in a hydrogen plasma, which is known to have a floating potential of -2.5kT/e, where k is Boltzmann’s constant, T is the temperature, and e is the elementary charge [Whipple]. The code for the example is given below:

from langmuir import *
from scipy.constants import value as constants
from scipy.optimize import root_scalar

n=1e11
T=1000
e = constants('elementary charge')
k = constants('Boltzmann constant')

plasma = [Electron(n=n, T=T),
          Hydrogen(n=n, T=T)]

geometry = Sphere(r=0.2*debye(plasma))

def res(V):
    return OML_current(geometry, plasma, V)

sol = root_scalar(res, x0=-3, x1=0)

print(sol.root)
print(-2.5*k*T/e)

The requirement for a steady floating potential is that the net current into the object, by all species, is zero. Otherwise the potential would increase or decrease. Thus, in our case we seek to find the root of the characteristic, which we take as OML_current in the example above. The function root_scalar from SciPy iterates on the value V of res within the range (-3,0), until it returns a value that is sufficiently close to zero, i.e., the root of OML_current. The script returns:

-0.2157892685234552
-0.21543333155362945

and hence is in excellent agreement with previous results.