getelec
GETELEC -- General Tool for Electron Emission Calculations.
Thermal-field electron emission current density and Nottingham heat for metallic and semiconducting emitters.
Quick start
The common case is one line, and any argument accepts an array:
>>> import getelec # doctest: +SKIP
>>> getelec.current_density(field=5.0, work_function=4.5) # doctest: +SKIP
>>> getelec.current_density(field=np.arange(3, 8, 0.5)) # doctest: +SKIP
A sweep is solved as a single batched call, not a Python loop, so it is substantially faster than calling the scalar version repeatedly.
For distributions, or to reuse one configuration across many calls:
>>> emitter = getelec.metal_emitter(work_function=4.5, # doctest: +SKIP
... fermi_level=7.5, temperature=300)
>>> emitter.calculate_current_density() # doctest: +SKIP
>>> energies, ted = emitter.calculate_total_energy_distribution() # doctest: +SKIP
Full control
The convenience layer just assembles four interchangeable components. Build them yourself whenever you need something the shortcuts do not expose:
potential_barrier -- SchottkyPotential, SmallRadiiPotential, TriangularPotential,
Customised (a potential of your own)
band_structure -- Metal, SmartMetal, Semiconductor, SmartSemiconductor,
CustomMetal, CustomSemiconductor,
DensityOfStatesMetal (a tabulated DOS)
electron_supply -- FermiDirac, LogFermiDirac
transmission_solver -- Noumerov (exact), NoumerovReference (written out
plainly, for checking), NeuralSolver (learned)
transmission_solutions -- WKB (fast, semiclassical), AiryTriangular
>>> from getelec.electron_emitter import MetalEmitter # doctest: +SKIP
>>> MetalEmitter(potential, solver, supply, band) # doctest: +SKIP
Units
Inputs are eV, nm, fs, K, and V/nm for the field. Current density comes back in A/cm^2, Nottingham heat P_N in W/cm^2, and distributions in A/(eV cm^2).
Citing
If you use GETELEC, please cite the software, S. Barranco Cárceles, A. Kyritsakis and A. Ayari, GETELEC: General Tool for Electron Emission Calculations, Zenodo, https://doi.org/10.5281/zenodo.23093209, and the papers listed in README.md and CITATION.cff.
Build a metal emitter with sensible defaults.
Parameters
work_function : float
Surface work function, eV.
fermi_level : float
Fermi level measured from the bottom of the conduction band, eV. Set
here once; it is propagated to the barrier and the supply function
together, which is the pairing most easily got wrong by hand.
temperature : float
Emitter temperature, K.
field : float
Local electric field at the surface, V/nm.
barrier : {'schottky', 'sharp_tip', 'triangular'} or potential_barrier.Barrier
Planar image-charge barrier, one corrected for tip curvature, or the
triangular barrier with no image charge, for which method='airy'
is exact. A barrier object is used as it is given, most usefully
~getelec.potential_barrier.Customised, which wraps a potential
of your own; fermi_level, work_function, field and
temperature are then set on it, as on every other component.
method : {'noumerov', 'ml', 'wkb', 'airy'}
'noumerov' integrates the Schrodinger equation, with no
approximation beyond discretisation. 'ml' uses ~getelec.transmission_solver.NeuralSolver, a
trained network with models shipped for both barriers. It pays off for
barriers with several parameters; for the planar barrier it is a
worked example, and getelec.training trains one for your own
conditions. 'wkb' is semiclassical: cheapest, but accurate only to
a factor of order unity deep in the tunnelling regime.
fast : bool, default False
Use Noumerov.fast(), which is within a small factor of the WKB cost
while keeping current density and energy distributions converged. See its
docstring for the exact accuracy contract.
reference : bool, default False
Use ~getelec.transmission_solver.NoumerovReference, the same
integration written out plainly, one energy at a time, to check a
result: about half a second for a metal current density. It takes
x_start, x_end, h and seed. Cannot be combined with
fast.
radius, gamma : float
Tip radius (nm) and field enhancement, used only by 'sharp_tip'.
The barrier is valid for radii of 20-1000 nm.
energy_resolution : float
Energy grid spacing, eV. Ignored when band is given, since a band
structure carries its own.
band : band_structure.BandStructure, optional
The band structure to use instead of the default
~getelec.band_structure.SmartMetal. Pass a
~getelec.band_structure.DensityOfStatesMetal to weight the
emission by a tabulated density of states; fermi_level and
work_function are then set on it as well, so its normalisation
follows a sweep.
**solver_kwargs
Passed to the solver, e.g. h=5e-4 to refine the spatial grid.
Returns
MetalEmitter
Examples
>>> em = metal_emitter(work_function=4.5, temperature=1000) # doctest: +SKIP
>>> em.calculate_current_density() # doctest: +SKIP
>>> em.update_params(field=6.0) # doctest: +SKIP
Build a semiconductor emitter with sensible defaults.
Defaults describe silicon. Parameters shared with metal_emitter()
behave the same way, with one exception: a barrier object passed as
barrier is used exactly as it is, and none of the Fermi level, work
function, field or temperature is set on it. Give a
~getelec.potential_barrier.Customised barrier those values
yourself.
Parameters
top_valence : float
Energy of the valence band maximum, eV, measured from the same zero as
fermi_level.
band_gap : float
Band gap, eV.
electron_eff_mass, hole_eff_mass : float
Effective masses relative to the free electron mass.
energy_resolution : float, default 0.002
Energy grid spacing, eV. Five times finer than for a metal, because
each band's grid ends at its band edge, where the TED has a finite
slope: the trapezoid error in the current is then ~(h / k_B T)^2 / 12,
about 1.1% at 0.01 eV and 300 K, and 0.05% at 0.002 eV. A metal's grid
ends in tails where the TED is flat, which is why 0.01 suffices there.
Returns
SemiconductorEmitter
Emitted current density in A/cm^2.
Any of field, work_function, fermi_level and temperature may
be arrays; they are broadcast against each other and the whole set is solved
in one batched pass.
Parameters
field : float or array_like
Local surface field, V/nm.
work_function : float or array_like
Work function, eV.
fermi_level : float or array_like
Fermi level, eV.
temperature : float or array_like
Temperature, K.
**kwargs
Forwarded to metal_emitter() (barrier, method, h, ...).
Returns
float or np.ndarray Scalar if every input was scalar, otherwise the broadcast shape.
Examples
>>> current_density(field=5.0) # doctest: +SKIP
>>> current_density(field=np.linspace(3, 8, 50)) # doctest: +SKIP
>>> current_density(field=5.0, temperature=[300, 800, 1500]) # doctest: +SKIP
Nottingham heat P_N in W/cm^2: negative when the surface heats (electrons leave from below E_F and are replaced by hotter ones), positive when it cools.
Broadcasts over its arguments exactly like current_density().
Returns
float or np.ndarray
Thermal-field emission from a metal surface.
Parameters
potential : potential_barrier.Barrier Barrier profile (work function, field, Fermi level). solver : transmission_solver.TransmissionSolver How transmission probabilities are evaluated. supply : electron_supply.Supply Electron supply function. band : band_structure.BandStructure Energy grid generator.
Examples
>>> from getelec import potential_barrier, band_structure # doctest: +SKIP
>>> from getelec import transmission_solver, electron_supply # doctest: +SKIP
>>> emitter = MetalEmitter( # doctest: +SKIP
... potential_barrier.SchottkyPotential(7.5, 4.5, 5.0),
... transmission_solver.Noumerov(),
... electron_supply.LogFermiDirac(7.5, 300.0),
... band_structure.SmartMetal(),
... )
>>> emitter.calculate_current_density() # doctest: +SKIP
Transmission probability D(E) on the emitter's energy grid.
Returns
energies, transmission : np.ndarray
Supply function l(E) on the emitter's energy grid, in eV.
Returns
energies, supply : np.ndarray
Total energy distribution, in A/(eV cm^2).
TED(E) = R(E) f(E) * integral of D(E_z) dE_z up to E -- the bare
occupancy times the transmission integrated over normal energy. See
~getelec.electron_supply.Supply.get_occupancy() for why this uses
f while the normal distribution uses l.
R(E) is the band structure's state weight, one everywhere for a free
electron gas and so for every band structure but
~getelec.band_structure.DensityOfStatesMetal.
Returns
energies, ted : np.ndarray
Normal energy distribution, in A/(eV cm^2).
NED(E_z) = S(E_z) * D(E_z). For a metal the effective mass equals
the free mass, so the transverse window is unrestricted and this is the
whole of it.
For a free electron gas S is the log supply l(E_z). When the
band structure weights its states -- a tabulated density of states --
the supply is instead integrated over the weighted occupancy above
E_z, which is what keeps integral NED dE_z == integral TED dE.
See get_weighted_normal_supply().
Returns
energies, ned : np.ndarray
Nottingham heat P_N, in W/cm^2.
The net power carried away by the emitted electrons relative to the replacement energy, taken at the Fermi level:
P_N = integral of (E - E_F) * TED(E) dE .
Negative means the emitter heats, positive that it cools. An electron leaving from below E_F is replaced by a hotter one at E_F, so cold field emission heats the tip (P_N < 0); at high temperature the emission moves above E_F and the tip cools (P_N > 0).
Returns
float
Notes
Reported in W/cm^2, the same area unit as the current density: the TED is in A/(eV cm^2), so its energy moment is in W/cm^2 directly.
Thermal-field emission from a semiconductor surface.
Conduction and valence band channels are treated separately, each with its own effective mass, and summed.
Parameters
potential : potential_barrier.Barrier solver : transmission_solver.TransmissionSolver supply : electron_supply.Supply band : band_structure.Semiconductor
The wavefunction for each band, on that band's energy grid.
Parameters
energies : tuple of array_like, optional
(conduction, valence). Defaults to the emitter's own grids.
x_points : array_like, optional
Positions in nm, shared by both bands.
Returns
energies_cb, x_cb, psi_cb, energies_vb, x_vb, psi_vb : np.ndarray Each band gets its own spatial grid: the domain is sized to enclose the barrier down to the lowest energy requested, and the two bands start at different energies.
Probability current for each band, on that band's grid.
A stationary scattering state carries a current that does not depend on position, and the Noumerov scheme does not enforce that -- so how much this varies across the grid is an independent check on the integration, band by band.
Returns
energies_cb, x_cb, current_cb, energies_vb, x_vb, current_vb : np.ndarray
Transmission for each band, on its own energy grid.
D(E) itself, the tunnelling probability, in [0, 1]. The effective masses
do not appear here: they set the window limits of Eq. (8) that
calculate_total_energy_distribution() integrates D between.
Returns
energies_cb, trans_cb, energies_vb, trans_vb : np.ndarray
Transmission integrated over the emission window, for each band.
g(E) = integral of D(E_z) dE_z from lower(E) to E, the
window of Eq. (8), in eV -- a transmission times the width of the
window it was integrated over, so it is not a probability and is not
bounded by 1.
Why it exists. calculate_transmission_coefficient() returns
D(E), which does not depend on the effective masses at all: the
barrier contains no mass. The masses enter only through the limits this
integral runs between, so g(E) is the smallest object that shows
how the transmission enters the current as the masses change it. A
heavier mass opens the window wider and raises g, which is the same
direction the current moves in.
TED(E) = f(E) g(E) -- this is the same g the distributions and
the current are built from, read out rather than recomputed, so the
plotted curve cannot drift from the physics.
Not to be confused with the by-parts integrand g'(E) of Eq. (15),
which is negative throughout the valence band for every hole mass and
is not returned anywhere.
Returns
energies_cb, window_cb, energies_vb, window_vb : np.ndarray
Supply function l(E) for each band, on its own energy grid.
Returns
energies_cb, supply_cb, energies_vb, supply_vb : np.ndarray
Total current density summed over both bands, in A/cm^2.
The integral of the total energy distribution, band by band.
Returns
float
Total energy distribution for each band, in A/(eV cm^2).
TED(E) = f(E) * integral of D(E_z) dE_z over the band's window,
evaluated directly rather than recovered from its derivative. See
_band_window_distributions().
Returns
energies_cb, ted_cb, energies_vb, ted_vb : np.ndarray
Normal energy distribution for each band, in A/(eV cm^2).
NED(E_z) = D(E_z) * [ l(E_z) - l(E_upper) ] -- Stratton Eq. (95) for
the valence band, Eq. (49) for the conduction band. Non-negative by
construction.
Returned on its own E_z grid, which extends below the
total-energy grid. For the conduction band with m_e* > m this puts
part of the NED in the band gap, and that is correct: E_z is the
normal energy in vacuum, and the conserved parallel momentum carries
up to (m_e*/m)(E - E_C) of transverse energy there, more than the
electron's whole kinetic energy in the band. Truncating at E_C
would lose that part of the current. See the GUIDE section
"Semiconductors: the two energy distributions" and
_band_window_distributions().
Returns
energies_cb, ned_cb, energies_vb, ned_vb : np.ndarray
Nottingham heat P_N summed over both bands, in W/cm^2.
P_N = integral of (E - E_R) * TED(E) dE with the replacement energy
E_R at the Fermi level, per Eqs. (13) and (19). Negative means the
emitter heats, positive that it cools, exactly as for a metal: an
electron leaving from below E_F is replaced by a hotter one.
At low field the valence band leads: electrons leave from below E_F, so the tip heats (P_N < 0). As the field rises the conduction band, above E_F, takes over and the sign flips to cooling (P_N > 0).
Returns
float
Inherited Members
Transmission from a trained neural network.
The network never predicts D. It predicts the residual
R = ln D_noumerov - ln D_kemble
against a semiclassical reference for the same barrier. The reference
already carries the exponential -- hundreds of e-folds -- so the network
only learns a correction of order unity, and the error budget is explicit:
relative error in D is exactly absolute error in R.
When it is worth it. For the planar Schottky barrier the network is
several times faster than Noumerov.fast(), but that calculation was
already cheap: the model there is the worked example, not the reason. Its cost does not grow with the number
of barrier parameters, while exact solves or a table of them do, so it pays
off for barriers with several -- tip radius and field enhancement, as in the
shipped curved-tip model. The shipped Schottky model is there as a running
example of the workflow; getelec.training trains one for your own
barrier.
Parameters
model : None, str, Path or NeuralModel, optional
None picks the shipped model for the barrier's class (see
SHIPPED_MODELS). A path or a NeuralModel uses that
model, for barriers of the class it was trained on.
fallback : {'exact', 'wkb', 'error'}, default 'exact'
What to do for a barrier with no model, or outside the trained domain.
'exact' solves it with Noumerov, keeping the answer right at the
cost of speed. 'wkb' uses the numerical semiclassical result.
'error' raises, which is what a pipeline that must not silently
change accuracy wants.
exact_solver : optional
Used by the 'exact' fallback. Defaults to Noumerov.fast().
energy_nodes : int, default 48
Evaluate on this many energies and spline ln D onto the rest. Both
the reference and the forward pass cost time proportional to the number
of energies, and ln D is smooth in energy, so this makes the cost
flat in grid size and changes the answer by well under 1%. Set to 0 to
evaluate at every energy.
Examples
>>> solver = NeuralSolver() # doctest: +SKIP
>>> barrier = SmallRadiiPotential(7.5, 4.5, 5.0, radius=50.0, gamma=100.0)
>>> solver.calculate_transmission(barrier, energies) # doctest: +SKIP
Whether a model exists for this barrier and its parameters were trained on.
Energies whose barrier height falls inside the trained range.
Transmission probability at each energy.
Inherited Members
Semiclassical (WKB) transmission for the Schottky-Nordheim barrier.
The Gamow exponent is
G = sqrt(2m)/hbar * integral from x1 to x2 of sqrt(V(x) - E) dx,
with V(x) = E_F + phi - F x - k_e / (4x) and x1, x2 the classical
turning points, i.e. the roots of F x^2 + (E - E_F - phi) x + k_e / 4.
Transmission then follows from the Kemble form, T = 1 / (1 + exp(2G)).
The integrand is sqrt(F (x2 - x)(x - x1) / x), which vanishes like a
square root at both turning points. Gauss-Chebyshev quadrature of the second
kind carries exactly that weight, so it converges geometrically where an
adaptive general-purpose rule fights the endpoint behaviour: 64 nodes give
about 12 significant figures. It is also fully vectorised over energy, with
no Python-level loop.
Above the barrier top the turning points become complex. There the exponent
is continued using the inverted-parabola (Kemble) approximation at the
barrier maximum, so T -> 1 smoothly instead of the formula breaking down.
Parameters
fermi_level, work_function, electric_field : float Barrier parameters (eV, eV, V/nm). n_nodes : int, default 64 Quadrature nodes. 64 is converged to round-off for typical parameters.
Gamow exponent G at each energy (dimensionless).
Exposed separately because G, not T, is the quantity that appears in Fowler-Nordheim analysis and in slope/intercept fitting.
Parameters
energies : np.ndarray potential : potential_barrier.Barrier, optional If given, its parameters take precedence over the ones stored on this solver.
Transmission probability at each energy, via the Kemble form.
Inherited Members
Exact transmission through a pure triangular barrier, via Airy functions.
For V(x) = E_F + phi - F x outside the metal and V = 0 inside, with
no image term -- the barrier of
~getelec.potential_barrier.TriangularPotential -- the Schrodinger
equation outside is Airy's equation, so the transmission has a closed form.
A plane wave in the metal, matched at the surface to the outgoing wave
Bi + i Ai outside, gives
T = (4k / (pi beta)) / [ (k/beta)^2 (Ai^2 + Bi^2) + Ai'^2 + Bi'^2 + 2k / (pi beta) ]
with k the wavevector in the metal, beta = (F / (hbar^2/2m))^(1/3),
and the Airy functions evaluated at the surface, zeta = beta (E_F + phi - E) / F.
The last term is the Wronskian Ai Bi' - Ai' Bi = 1/pi. It is negligible in
deep tunnelling, where Bi is exponentially large, and it is what takes
the result to the step-barrier value 4 k q / (k + q)^2 far above the
barrier, q being the wavevector just outside the surface.
This is not the barrier a metal presents -- dropping the image charge removes the Schottky lowering, so it overestimates the barrier and underestimates emission. Its value is as a check: it is the one case where an exact answer exists, so it pins the numerical solver against something that is not another numerical solver.
Parameters
fermi_level, work_function, electric_field : float eV, eV, V/nm. Read from the barrier when one is passed.
Transmission probability at each energy.
Inherited Members
Transmission coefficient D(E), the bare tunnelling probability.
Parameters
field, work_function, fermi_level, temperature : float
Barrier and emitter parameters. temperature only sets the default
energy grid; D itself does not depend on it.
energies : array_like, optional
Energies at which to evaluate D, in eV. Defaults to the emitter's own
grid, which is chosen to cover everything that carries current.
**kwargs
Forwarded to metal_emitter() (method, barrier, ...).
Returns
energies, transmission : np.ndarray
Examples
>>> E, D = transmission_coefficient(field=5.0) # doctest: +SKIP
>>> E, D = transmission_coefficient(field=5.0, method="ml") # doctest: +SKIP
Electron supply function N(E).
Electrons arriving at the barrier per unit energy. Independent of the field and of the barrier shape; it depends only on temperature and Fermi level.
Returns
energies, supply : np.ndarray
Put the GETELEC attribution in the bottom-right corner of a plot.
Parameters
target : matplotlib Axes or Figure, optional
Where to draw. Defaults to the current axes. Passing a Figure places the
mark once for the whole figure rather than once per panel.
text : str, optional
Overrides WATERMARK_TEXT.
alpha, size : float
Opacity and font size.
Returns
matplotlib.text.Text
Examples
>>> import matplotlib.pyplot as plt # doctest: +SKIP
>>> fig, ax = plt.subplots() # doctest: +SKIP
>>> ax.plot(fields, current) # doctest: +SKIP
>>> getelec.watermark(fig) # doctest: +SKIP