General

Introduction

Bumblebee utilizes kinetic Monte Carlo simulations to model the opto-electronic properties of OLED stacks. The transport of the electrons through the device is described stochastically in order to determine the time-evolution and space-dependence of excitonic processes.

Input parameters can be obtained from

  • the built-in materials database, or

  • through the AMS OLED workflow, or

  • through manual input.

Bumblebee may be used

  • from the command line,

  • from Python, or

  • through a dedicated web interface.

Features

Bumblebee is capable of simulating OLED stacks containing numerous layers comprised of both pure and composite materials.

  • An extensive array of excitonic interactions can be included to describe the opto-electronic character

  • Device degradation can be modeled to estimate stack lifetimes

  • Photoluminescent absorption allows replication of dark-field and PV studies

  • Optical outcoupling simulations predict emission spectra at the device level

  • Small-signal analysis with AC modulation

Theory

Kinetic Monte Carlo 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

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. These crossing events are currently implemented with a static rate (independent of temperature or site environment).

Exciton 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

In order to compute the effect of dipole orientations on the electric field, molecular dipoles may be enabled. Random orientations of the dipoles are generated at each carrier site. A material-specific dipole magnitude is specified. Dipole orientations are randomly drawn from an ellipsoidal distribution, which can be modified to describe orientational preference.

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.

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\]

Alternating Voltage

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

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

The OFET model assumes that the electrode contacts act as gates. The source and drain are 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 of electrode Fermi levels

  • 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.

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

When the layer composition is made up of a mixture of different materials, a random distribution of components in the layer is assumed. A clustering modification may be defined to mimic self-aggregation.

Advanced compositions can be defined in order to generate more complex distributions. Linear and trapezoidal 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 linear profile provides the average composition along the device stack. The distribution of components between channels remains random.

For polymeric materials, a polymer chain morphology 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 lock backwalking until a minimum number of growth steps has occurred.

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

Kinetic Monte Carlo 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.

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

  • 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.

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. open circuits)

Simulation progress is written to the output files at fixed intervals. Default settings in the web interface 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.

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

Acceleration

A rate booster is 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.