Theory
******

Kinetic Monte Carlo (kMC) estimates the properties of the OLED stack
through stochastic sampling. A model of the stack is formulated as a
3D grid, with each node representing a charge carrier site. Electrode
contacts are placed at the edges of the device to allow for current
flow. The probability for a particular opto-electronic process to
occur, such as exciton generation or charge transfer, is taken to be
proportional to the rate of that process. Sampling of the various
mechanisms mimics the time-evolution of the OLED, allowing kMC to
explore the relevant regimes encountered during device operation.

In order to accurately model real OLED devices, relevant processes
must be included in the simulation. Various mechanisms are included in
Bumblebee.

* Charge transfer

* Exciton generation and emission

* Exciton annihilation

* Charge injection/removal at the electrodes

* Intersystem crossing

* Exciplexation

* Degradation


Exciton Generation
==================

The HOMO and LUMO energy levels are considered as the energy levels of
the charge carriers. Due to the inhomogeneous nature of solid-state
organic electronics, the HOMO and LUMO values tend to vary between
sites. This leads to a layer-wide distribution of energy levels.
Bumblebee includes various model distributions to represent this
variance. A standard deviation can be specified for either level.

Excitons are restricted by Bumblebee to either be in a singlet or
triplet state, as higher spin states are short-lived in typical OLED
materials. The ratio of singlet to triplet excitons can be adjusted
for each layer.

The singlet and triplet species each have a binding energy. This value
is combined with the HOMO and LUMO levels to determine the singlet and
triplet energy levels.

   E_{ex} = E_{LUMO} - E_{HOMO} - E_b

With E_{ex} the exciton energy level and E_{b} the binding energy.
Similar to the HOMO and LUMO energies, a distribution of exciton
energy levels may also be specified.

By default, the distribution in the energy levels is assumed to be
uncorrelated. A correlation between the HOMO and LUMO levels can be
specified to preserve the band gap. Anti-correlation is also
available, wherein the energy shift in the HOMO is opposite that of
the LUMO.

The exciton distribution may in turn be correlated to the band gap:

   \epsilon_{ex} = \sigma_{ex}
   \frac{\epsilon_{gap}}{\sqrt{\sigma_{HOMO}^2+\sigma_{LUMO}^2}}

The binding energy of the excitons is, by default, unaffected by the
energy level shifts. Automatic propagation of the exciton shifts to
the binding energies may be enabled.

Exciton generation may occur either through charge transfer or thermal
excitation. The latter mechanism, in the absence of photoluminescence,
is described through a Boltzmann process. The associated prefactor is
considered a static parameter, which is constant throughout the stack.


Energy Distribution Functions
-----------------------------

Several model distributions are available when sampling the energetic
disorder in energy levels.

   E = E_0 + P\left(X\right)

* Dirac delta

   P(X) = 0

* Gaussian

   P\left(X\right) =
   \frac{\sigma}{\sqrt{2\pi}}\exp\left(-\frac{1}{2}X^{2}\right)

* Exponential

   P\left(X\right) = \frac{q}{\left|q\right|}k_{b}T\exp\left(-X\right)

* Mirrored exponential

   P\left(X\right) = \frac{\left(X_1 - \frac{1}{2}\right)}{\left|X_1 -
   \frac{1}{2}\right|}k_{b}T\exp\left(-X_2\right)

* Chi-squared-extended Gaussian

   P\left(X\right) = P_{\textrm{Gaussian}}\left(X_1\right) + \frac{q}
   {\left|q\right|}W\frac{\sqrt{X_2}\exp\left(-\frac{1}{2}X_2\right)}
   {\Gamma\left(\frac{3}{2}\right)2^{3/2}}

* Dipole-correlated Gaussian

Bumblebee supports the inclusion of explicit dipole fields associated
with the molecular dipole moments in the stack. When a dipole-
correlated energy-level distribution is requested, the average dipole
potential is evaluated at each site. This field distribution is then
normalized to obtain the site-specific energy shifts:

   \sigma_A = \frac{\sum_j \frac{\bar{p_j}\cdot\left(\bar{p_A}-\bar{p
   _j}\right)}{\left[\left(\bar{p_A}-\bar{p_j}\right)\cdot\left(\bar{
   p_A}-\bar{p_j}\right)\right]^{3/2}}}{\sum_i \sum_j \frac{\bar{p_j}
   \cdot\left(\bar{p_i}-\bar{p_j}\right)}{\left[\left(\bar{p_i}-\bar{
   p_j}\right)\cdot\left(\bar{p_i}-\bar{p_j}\right)\right]^{3/2}}}


Charge Hopping
==============

The transport of charge carriers through the layer is modeled as a
hopping process between sites. The probability of electron hopping is
described using a Boltzmann factor.

   r = \nu \exp\left( \frac{-\Delta E}{k_b T} \right)

\Delta E represents the change in energy following the hop, k_b T is
the thermal energy and \nu is the hopping frequency. The frequency is
assumed to decay exponentially with the inter-site distance. For hops
between layers, a geometric average of the layer parameters is
assumed:

   \nu = \sqrt{\left( \nu_I\nu_F \right)} \exp\left( -\left(R-a\right
   )\left(\frac{1}{\lambda_I}+\frac{1}{\lambda_F}\right) \right)

\nu_l represents the reference hopping frequency between adjacent
sites at distance a in each layer. \lambda_l is the decay length
characterizing the decrease in hopping frequency with distance.

Note:

  Since many opto-electronic processes occur at very high frequencies,
  a reference timescale is introduced to re-scale the rate constants.
  See the section on timescales for more info.

Due to the exponential decay, hopping events will be localized near
each charge carrier site. To accelerate calculations, the number of
hopping events is restricted to a local environment. Per default, a
cubic volume with an edge-length of 4a is employed.

For charge exchange with the electrodes, the electrode Fermi level is
used to determine the change in energy. Charge transfer with the
electrodes is restricted to the boundary sites. The rate of charge
exchange with the electrodes may be modified in order to bias the
sampling distribution. For typical simulations, frequent electrode
exchange events may result in long computation times. Reducing
exchange rates then enhances the sampling of events in the stack. Care
should be taken that the exchange adjustments do not affect the
physical output of the simulation.

In order to investigate the effect of trap states, dopants or
contaminants, device-specific prefactors may be appended to the
hopping frequencies. These prefactors allow modification of the
hopping rates across the stack without having to modify individual
material properties. For more detailed studies, trap states or dopants
can also be explicitly included in the simulation.

Anisotropic charge transport may additionally be specified by
introducing an anisotropy vector:

   \nu_{ANI} = \left(\frac{\left|R_x\right|\mu_x^{ANI} +
   \left|R_y\right|\mu_y^{ANI} + \left|R_z\right|\mu_z^{ANI}}{|R_x| +
   |R_y| + |R_z|}\right)\nu

In excitonic simulations, charge hopping may result in formation of a
new exciton following polaron collision. An additional collision
prefactor can be introduced to alter the frequency of these events.


Exciton Hopping
===============

Hopping of an exciton may occur through both Dexter and Förster
mechanisms. The Dexter process involves charge transfer through an
overlap of the molecular wavefunctions, and is therefore necessarily
localized. Dexter transfer is described through the same hopping
mechanisms as the charge carrier species, with an associated hopping
frequency and decay length.

The Förster mechanism involves exciton transfer through dipolar
coupling:

   \nu = \left(\frac{F}{R}\right)^6

F denotes the Förster radius, characterizing the coupling strength.

Note:

  All hopping parameters are unique for singlet and triplet excitons.

The thermal dependence of the Förster rates can also be disabled,
turning the Förster reactions into fully electronic processes.


Exciton Decay
=============

Excitons may decay radiatively (resulting in photonic emission) or
non-radiatively (generating heat). Rates are specified for both
processes. These rates may differ between singlet and triplet species.


Intersystem Crossing
====================

Intersystem crossing and reverse intersystem crossing processes can be
included in the Bumblebee simulation to account for singlet-triplet
interconversions. By default, a constant rate is specified for these
crossing events. Thermally activated intersystem crossing can be
enabled when studying, for instance, temperature-dependent
photoluminescence.


Exciton Quenching and Annihilation
==================================

Excitons may be converted through various mechanisms.

* Exciton-polaron quenching

* Singlet-singlet annihilation

* Singlet-triplet annihilation

* Triplet-triplet annihilation

Each of these reactions can be specified using both Dexter and Förster
models. The Dexter parameters are assumed to be material-specific,
therefore using the same values as for the exciton hopping. Förster
reactions meanwhile have to be specified for individual processes.

Two reaction pathways can be defined for exciton-polaron quenching:

* Polaron transfer to the exciton site

* Polaron-induced exciton recombination

Bumblebee allows specification of distinct rates for these steps.

Triplet-triplet annihilation may lead to formation of both singlets
and triplets. By default, a singlet-triplet ratio is specified for the
annihilation product. This behavior may be altered to specify singlet-
selective or triplet-selective products.

As with the polaron transport, a device-specific prefactor may be
appended to the exciton hopping frequencies to investigate the effect
of e.g. impurities. These prefactors can be defined for each
individual mechanism to aid in the investigation of mechanistic
sensitivity.

Since annihilation events involve electronic processes, they are
considered athermal by default. Use of the temperature-dependent
Boltzmann factor for annihilation reactions can be enforced in the
settings.


Interlayer Processes
====================

For Dexter processes that occur across layers, a geometric average is
used to determine the transfer parameters. Förster processes,
meanwhile, require the parameters to be specified for each unique pair
of materials. To add a finer degree of control, it is possible to
provide similar definitions for the Dexter prefactors as well.

Default constructors are available for the inter-layer Förster
parameters to ease simulation setup and to aid in mechanistic
screening studies. Material characteristics are used to determine the
relevant processes and to estimate the transfer rates.


Vibronic Coupling
=================

The default expression for thermal activation assumes a ground-state
occupation for both initial and final states.

   r_{Boltzmann} = \nu\exp\left(\frac{-\Delta E}{k_b T}\right)

Marcus theory allows incorporation of a vibrational distortion of the
molecular geometry following charge transfer, expressed in terms of a
corresponding energy shift \lambda.

   r_{Marcus} = \frac{ \nu\exp\left[\frac{-\left(2\Delta E + \lambda_I
   + \lambda_F\right)^2}{8\left(\lambda_I + \lambda_F\right)k_b
   T}\right] }{ \sqrt{2\pi\left(\lambda_I + \lambda_F\right)k_b T} }

The Levich-Jortner expression is used to incorporate tunneling.

   r_{Jortner} = &\frac{\nu\exp\left(\frac{-\left(\kappa_I +
   \kappa_F\right)}{\hbar\left(\omega_I +
   \omega_F\right)}\right)}{\sqrt{2\pi\left(\lambda_I +
   \lambda_F\right)k_b T}} \sum_{i=0}^{n}
   \frac{1}{i!}\left(\frac{\kappa_I + \kappa_F}{\hbar\left(\omega_I +
   \omega_F\right)}\right)^{i}\\ &\exp\left[\frac{-\left[2\Delta E+
   \lambda_I + \lambda_F + i\hbar\left(\omega_I+\omega_F\right)\right
   ]^{2}}{8\left(\lambda_I+\lambda_F\right)k_b T}\right]

\kappa denotes the associated zero-point relaxation energy, n denotes
the number of accessible vibrational energy levels and \omega the
normal mode frequency. A single, dominant vibration is assumed to
define the vibronic coupling. A multi-modal formalism may similarly be
enabled.

   r_{multimode} =
   &\frac{\nu\left[\prod_{j}^{m}\exp\left(\frac{-\left(\kappa_{I,j} +
   \kappa_{F,j}\right)}{\hbar\left(\omega_{I,j} +
   \omega_{F,j}\right)}\right)\right]}{\sqrt{2\pi\left(\lambda_I +
   \lambda_F\right)k_b T}} \sum_{\{i\}}^{\{n\}}\left(\prod_{j}^{m}
   \frac{1}{i_j!}\left(\frac{\kappa_{I,j} +
   \kappa_{F,j}}{\hbar\left(\omega_{I,j} +
   \omega_{F,j}\right)}\right)^{i_j}\right)\\
   &\exp\left[\frac{-\left[2\Delta E+ \lambda_I + \lambda_F + \sum_{j
   }^{m}i_j\hbar\left(\omega_{I,j}+\omega_{F,j}\right)\right]^{2}}{8\
   left(\lambda_I+\lambda_F\right)k_b T}\right]

For a set of m vibrational modes, each characterized by their own
frequency and relaxation energy.

Because the vibronic coupling in the multimode form involves computing
the interactions for each combination of harmonic overtones, this
easily become a computational bottleneck in the simulation. An option
is provided to compute the multimode rates using an iterative
refinement approach, significantly reducing the overall computational
cost, though some residual errors may remain in the system.


Coulomb Interaction
===================

The voltage (V) applied to the stack creates an external field that
affects the hopping rates:

   \frac{\Delta E_{ext}}{\Delta x} = \frac{\left(V + E_{fa} -
   E_{fc}\right)}{L}

With L the stack thickness and E_{f} the Fermi levels of the anode and
cathode.

Coulomb interactions between charges are added to the external field
contribution:

   E_C = \frac{q}{4\pi \varepsilon_r \varepsilon_0}\frac{1}{R}

With q the carrier charge, \varepsilon_r the relative permittivity and
\varepsilon_0 the vacuum permittivity. The pair-wise interactions are
computed explicitly for nearby charges within a local environment,
defined by a cutoff radius. Interactions outside the cutoff radius are
accounted for by considering a layer-uniform background charge
density. To include interactions beyond the simulated region,
accounting for the decay in Coulomb interactions with distance, the
long-range contribution is truncated to a finite number of periodic
images.

To accelerate calculations, particularly for homogeneous fields, an
update interval may be specified for the long-range contribution. This
assumption is valid when the change in the inter-layer charge
distributions is slow.


Initial Charge Distribution
===========================

At the start of the simulation, no charge carriers are present in the
layers. Charge carriers must be generated through excitonic processes
or from the electrode current.

Initialization of charge carriers may be enabled to preserve carrier
densities in e.g. bulk simulations. Charge carriers are then injected
randomly throughout the stack. Carrier parameters are available to
modify the material-specific injection probability. This allows
initial charge distributions to be biased towards specific layers.

Defects may furthermore be specified to adjust the carrier densities
in individual layers. Electrons, holes, singlet and triplet excitons
can be selected.


Dipole Distributions
====================

The orientation of the molecular dipoles can affect both the internal
electric field and introduce spatial correlations of the molecular
energy levels. Molecular dipoles can be enabled to incorporate these
effects in the simulation.

Each molecule is assigned a dipole moment magnitude. A dipole
distribution can also be specified for each individual material.
Molecular dipoles at each carrier site will then be randomly generated
in accordance with this distribution. Ellipsoidal distributions are
used by default, yielding net-zero polarization, which can be modified
to describe orientational preference and alignment of the dipoles.

Transition dipole fields may similarly be specified in order to
account for orientational alignment in Förster processes. The
transition dipole vector and the molecular dipole vector are assumed
to be uncorrelated.

Spontaneous orientation polarization (SOP) can be enabled to include
field alignment of the dipoles in OLED/OFET devices. Along with the
dipole distribution, the giant surface potential (GSP) effect will
then be automatically computed and included in the simulation.


Quadrupolar Interactions
========================

Förster processes are typically mediated through dipolar coupling.
Some OLED molecules may also exhibit transition quadrupole moments
that contribute to the intermolecular interactions.

Both dipole-quadrupole and quadrupole-quadrupole contributions can be
added to individual Förster processes. These quadrupolar interactions
are characterized by effective Förster radii:

   \nu = \left(\frac{F_M}{R}\right)^M

Where M is the multiple order: 8 for dipole-quadrupole interactions,
10 for quadrupole-quadrupole interactions.


Polymeric Materials
===================

Linear, polymeric organo-electronics typically exhibit different
charge transfer characteristics in parallel or perpendicular
directions to the backbone. Polymer chains can be defined in the layer
morphology, assigning each site to a particular chain. Inter- and
intra-chain charge transfer prefactors are then appended to the
hopping frequencies:

   \nu_{\parallel} = \mu_{\parallel}\nu \qquad \wedge \qquad
   \nu_{\perp} = \mu_{\perp}\nu


QLED and QNED Devices
=====================

QLED/QNED structures can be selected when specifying a layer
morphology. Different lattice/packing models for the embedded quantum
dots are available. Defect formation is also supported. Multi-layer
core-shell quantum dots can be defined simply by selecting the
required materials.

External charge transfer into the quantum dot (or between shells)
tends to differ from the particle-interior charge transport due to the
functionalization of the capping layers. Internal and external shell-
charge transfer prefactors are therefore added to the transfer rates:

   \nu_{int} = \mu_{int}\nu \qquad \wedge \qquad \nu_{ext} =
   \mu_{ext}\nu

Optical outcoupling calculations account for the chosen morphology of
the QLED/QNED layers in the device.


Exciplexes
==========

Exciplex formation involves a charge-transfer (CT) state between donor
and acceptor molecules. Compared to an exciton, the electron-hole pair
of the exciplex is localized over 2 sites, resulting in different
behavior.

* Because the charges are spatially separated, the exciplex
  experiences Coulombic interactions, in contrast to a charge-neutral
  exciton

* An exciplex has an electric dipole, resulting in anisotropic hopping
  of charges due to the orientation of the internal field

Exciplexes can participate in the same processes described for
excitons: emission, non-radiative decay, ISC/RISC, quenching and
annihilation.

Exciplexes can hop between donor-acceptor pairs. It is also possible
for exciplex diffusion to occur through concerted hops of the electron
and hole. Both mechanisms are supported. The degree of charge
delocalization of the CT state can be tuned to account for different
magnitudes of the electric dipole.


Alternating Current
===================

By default, Bumblebee operates under a direct current. An alternating
current may also be specified in terms of a modulated signal on top of
the bias voltage. This signal is described by an amplitude, angular
frequency, phase shift and a time offset relative to the start of the
simulation. To improve sampling statistics for high-frequency signals,
an AC update interval may be specified to dilate signal updates during
the kMC sampling.


Photoluminescence
=================

Photoluminescence simulations can be performed by enabling exciton
generation following photo-absorption. For each material in the stack,
an absorption coefficient is specified as an excitation probability. A
static device exposure is defined in terms of the number of incident
photons per second per cubic nanometer of material.

Excitations can be defined to give rise to various absorption
products. Excitons, exciplexes or polarons may be specified. For
exciton generation, a singlet-triplet ratio can be used to describe
the spin statistics of the excitation.


Lifetime Simulations
====================

In order to model the effect of molecular degradation on the lifetime
device performance, sites may be deactivated following certain opto-
electronic events. These include:

* Exciton generation

* Exciton annihilation

* Polaron quenching

* Polaron hopping

* Photo-absorption

* Emission

* Non-radiative decay

Each event has a fixed probability to result in degradation of the
material. When degradation occurs, the layer composition is altered by
replacing the original material with a degraded material specified by
the user. This mechanism can also be used to specify opto-
electronically induced chemical reactions inside the stack.


OFET
====

Bumblebee can also be used to model organic field-effect transistors
(OFET). The OFET model assumes that the electrode contacts act as
gates. The source and drain contacts of the OFET are then placed
perpendicular to the gate field. A voltage can be applied between the
source and the drain.

During OFET operation, a gate field is applied between the electrodes.
Both single- and dual-gate devices can be modeled. A zero-potential
reference is introduced in order to confine the gate field to the OFET
interior. (By default, an insulator is placed at the cathode.)
Corrections to the gate field are applied to preserve this potential.
Minor field fluctuations are typically ignored by specifying a
response threshold.

It is possible to select one of two control methods for preserving the
potential:

* Modification of the external gate field

* Charge carrier injection

To inhibit rapid fluctuations in the device model due to these
potential adjustments, gate updates are typically only performed at
set intervals. The magnitude of the external field updates is
restricted to a fixed field strength stepsize. Carrier injections are
limited to one polaron per update.


Transient Responses
===================

kMC can be used to observe the change in OLED behavior upon some
external stimulus. It is possible to specify a perturbation to the
system which takes place after some initial startup time.

* Modify the voltage, temperature, fluence or the work function of the
  electrodes

* Toggle degradation events

* Match the device voltage to a target current

Multiple perturbations can be defined (at different time offsets) to
model more complex responses. Perturbations can also be combined with
the alternating current mode to study response properties under
dynamic conditions.


Pulses
------

The checkpoints module is used to perform instantaneous updates to the
device parameters. Time-resolved input signals, more closely
resembling experiment, can be achieved using the pulses module.

The pulses module allows for both single pulses and periodic
modulation of specific system parameters. Multiple signals can also be
combined to describe complex responses.

* Various built-in waveforms are provided, including sine waves,
  square waves and exponential decay

* A rectifier and/or diode can be added to filter the signal source

* Interpulse and/or intrapulse intervals can be added for repeated
  pulsing experiments

Both checkpoints and pulses can be combined to describe dynamic or
multi-variate responses.

Note:

  For impedance calculations, it is important to consider that the
  external circuit is not included in the Bumblebee simulation. When
  comparing to experimental signals, external resistances and
  reflections should be accounted for.


Simulation Volume
=================

Bumblebee assumes a 1 nm distance between sites. The thickness of the
stack is obtained as the sum of the layers. The (rectangular) surface
area of the cell can be specified through the edge lengths.

Periodic boundary conditions may be specified to model bulk material
behavior. This setting preserves the bias voltage of the electrodes,
but disables external charge exchange.


Morphology
==========

Ideal mixtures are used by default to describe layers containing
multiple materials. Self-aggregation behavior can be specified to
model non-ideal distributions.

Position-dependent morphologies can be defined in order to generate
more complex material distributions. Linear, trapezoidal and
exponential gradients are provided to specify one-dimensional material
distributions. Multiple gradients can be combined to define the final
morphology. A background material is introduced by default to assure
compositional fractions always sum to unity.

The one-dimensional profiles provide the cross-sectional average
composition along the device stack. The distribution of components
between channels remains stochastic.

Several three-dimensional morphologies are also available:

* For polymeric materials, a polymer chain network can be defined.
  Chain growth is simulated using a self-avoiding walk. To define the
  morphology, an anisotropic propagation rate is defined along with
  the desired chain length. New chains are added to the layer until
  the polymer fraction is filled. Backbone rigidity can be specified
  to limit backwalking.

* QLED/QNED morphologies allow for the embedding of quantum dots into
  the active layer. Quantum dots can be composed of multiple shells,
  each with their own thickness and material composition. Different
  packing models (including both structured and disordered lattices)
  are available to describe the distribution of particles within the
  layer. Defects can also be added to describe the effect of
  imperfections in the lattice.

* Nanoparticle morphologies allow for the inclusion of nanoparticles
  or nanocrystallites into the layer. Polydisperse solutions can also
  be defined. Similar to the quantum dots, different packing models
  can be used to describe the particle distribution inside the layer.
  These morphologies allow for modeling of hybrid OLED (HyLED)
  utilizing e.g. perovskite crystals.

When three-dimensional morphologies (PLED, QLED, QNED, HyLED) are
used, the generated structures are embedded into a background matrix
generated by the the one-dimensional profiles. Multiple morphologies
can be combined in a single layer to create, for example, a QLED
inside a polymeric matrix.


Annealing
=========

In order to simulate the voltage ramp that occurs following circuit
closure, an annealing stage can be added to the simulation.

   V\left(t\right) = V + \left(N_{anneal} -
   \frac{t}{T_{ramp}}\right)\left(\frac{V_{anneal} -
   V}{N_{anneal}}\right)

A stepped ramp function is used to increment the voltage within a
timeframe of N_{anneal} simulation steps. The voltage is updated in
T_{ramp} increments.

During OFET simulations, the voltage ramp is applied to the transistor
field instead.


Pre-equilibration
=================

Bumblebee analyses trajectory data to determine the properties of the
OLED device. These trajectories depend on the initial state of the
system. Typically, the properties of interest are obtained at
equilibrium, while the initial state may be off-equilibrium. This
introduces a bias in the trajectory statistics.

A pre-equilibration stage may be specified in the simulation settings
to allow the system to de-correlate from the initial state,
approaching the equilibrium regime. Sampling statistics will then only
be generated for the samples obtained after pre-equilibration.


Optical Outcoupling
===================

Only a part of the light that is emitted by the molecules is able to
exit the device. Optical losses, including waveguiding, re-absorption,
plasmonics and evanescence limit the external efficiency of the
device.

The companion software Firefly allows for the calculation of the
optical outcoupling efficiency for devices simulated with Bumblebee.
Optical layers, including substrates, electrodes, coatings and color
filters can be added to the OLED stack to obtain a detailed picture of
the device optics.

Wavelength-dependent complex refractive indices can be provided to
calculate the distortion of the external emission spectrum. Transition
dipole anisotropy and detailed layer morphologies are automatically
accounted for.

Bumblebee calculates the CIE color points, emission spectrum
(including color blending), luminosity and power efficiencies using
the data obtained from Firefly. Angular emission spectra and Purcell
factors are also available.

When Firefly is not enabled, Bumblebee calculations assume a 20%
outcoupling efficiency, which falls in line with typical OLED device
structures.

Note:

  Firefly is included with the Bumblebee license. AMS automatically
  installs Firefly as part of the Bumblebee module.


Simulation Settings
===================


Timescale
---------

A Bumblebee simulation typically includes a large number of processes,
spanning various timescales. To simplify the input, a reference
timescale is introduced. This allows many of the rate constants to be
re-scaled to values close to unity.

For rates described using probabilistic distributions (such as charge
hopping or exciton annihilation), prefactors are expressed in units of
the reference timestep.

Note:

  Assuming a reference timestep of 10^{-11}\,\textrm{s}, a prefactor
  of 0.9 would correspond to a process with a frequency of
  0.9\cdot{}10^{11}\,\textrm{s}^{-1}.

Fixed-rate processes may be provided in natural units, and are re-
scaled internally:

* ISC/RISC

* Radiative and non-radiative exciton decay

* Radiative and non-radiative exciplex decay

* Photoluminescent fluence

Note:

  Assuming a reference timestep of 10^{-11}\,\textrm{s}, an ISC rate
  of 10^{8}\,\textrm{s}^{-1} is converted to a probabilistic rate of
  10^{-3} occurrences per simulation step.

By adjusting the reference timeframe of the simulation, this changes
the magnitude of the rate values used internally by Bumblebee. Using
rate values closer to 1 can improve the numerical stability of the
simulation. For typical OLED devices, it is recommended to use the
default settings.

Warning:

  Changing the reference timescale requires updating all the rate
  constants in the input. (The fixed rates remain unaltered.)


Parallelization
---------------

The accuracy of the statistical estimates provided by the kMC method
depends on the number of samples. Typically, a very large number of
samples is required to obtain reliable results. In order to improve
performance, parallel sampling can be performed by selecting multiple
simulation volumes. Samples from these parallel trajectories can be
collated to obtain statistics for the aggregate dataset.

In BBinput, the number of trajectories can be selected. Multiple
parallel runs will then be started automatically. BBresults will
collect the results from the trajectories, and also provide error
estimates for the reported data.


Output
------

As the kMC simulation progresses, the statistical estimate of the
device performance is iteratively improved. In theory, the sampling
could continue until all states have been exhausted. In order to
define a practical termination condition for the simulation, a maximum
number of sampling steps is defined.

Alternative conditions are also provided:

* The simulation can terminate once a certain lifetime has been
  evaluated. (The time progression for each sample in the simulation
  is obtained as the inverse of the process frequencies)

* The simulation can terminate once a homogeneous charge density has
  been achieved. This is expressed in terms of the normalized midrange
  of the 1D current profile:

   \frac{I_{max} - I_{min}}{I_{avg}} < \tau_{conv}

* The simulation can terminate once all the excitons and/or charge
  carriers have been depleted. (As may be of interest for e.g. TRPL
  measurements)

Simulation progress is written to the output files at fixed intervals.
Default settings in the GUI automatically configure the relevant
output. Users may modify the requested properties as desired. The
simulation time increases as more output files are enabled, on account
of the time needed to write them.

The output frequency may be configured separately for output regarding
sampling statistics, or output related to device profiles. This allows
the more expensive files to be written less frequently.

Note:

  For simulations that are dominated by high-frequency events, small
  time intervals will be used by the kMC simulation. A larger number
  of significant digits for the timestamps may therefore be enabled in
  the output settings.


Acceleration
------------

Bumblebee features an accelerated kMC algorithm, which requires fewer
samples to obtain accurate device statistics. The acceleration module
is enabled by default in order to yield the best simulation
performance.

The acceleration module enhances the process sampling. Conservation of
the process distribution compared to the non-accelerated reference
system is imposed to guarantee that the kMC simulation converges to
the same estimates.


Rate Booster
~~~~~~~~~~~~

A rate booster is also available to promote anisotropic charge hopping
along the electrode current. This accelerates the sampling of charge
transfer through the stack. Note that these parameters modify the
intrinsic rate distributions specified earlier. As the bias affects
the sampling process, this invariably has an effect on the resulting
statistics. Care should be taken that the boosting method preserved
the relative rates encountered in the unbiased simulations.


Memory Usage
------------

Bumblebee has to store all polaron and exciton species in the
simulation volume. Since most OLED devices operate under sub-saturated
carrier densities, a maximum number of carrier species may be
specified to reduce memory consumption.

Dynamic memory allocation is performed whenever this limit is
exceeded, such that simulations are never terminated prematurely. As
additional allocations are quite slow compared to the normal
initialization, it is recommended to set memory limits appropriately.


Reproducibility
---------------

The stochastic sampling process occurs through the use of random
number generation. If reproducibility of the simulation trajectories
is desired, a fixed RNG seed may be specified.
