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.

\[\begin{split}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]\end{split}\]

\(\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.

\[\begin{split}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]\end{split}\]

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.