Fluorescence enhancement

The enhancement of the emission rate of a quantum emitter due to the interaction with an optical antenna structure can be simulated in FDTD. In this example we show how to calculate the normalized radiative and loss decay rates, as well as the quantum efficiency. Results are validated against an analytic solution for a metal nanosphere antenna with radius much smaller than the excitation wavelength.

Overview

Understand the simulation workflow and key results

Fluorescence occurs when quantum emitters (e.g. molecules or quantum dots) are excited by absorption of a photon and then later decay by spontaneous emission of another photon (typically of longer wavelength) and by non-radiative relaxation mechanisms (such as phonon excitation). Nearby structures can enhance the fields that cause the excitation as well as enhance the decay rate due to spontaneous emission. Both types of enhancement can be calculated with FDTD. The present example explains how to calculate the decay rate enhancement with respect to spontaneous emission in the absence of the structure.

Key objects in the simulation are:

• Dipole source located at the emitter position. In the FDTD simulation this dipole source models the radiation characteristics of the emitter.
• Nearby physical objects (acting as antenna), such as the nanosphere shown here.
• Box monitors around the dipole and the entire system.

The box monitors measure the amount of power emitted by the dipole, the power absorbed by the nearby objects and the power emitted to the far field. These quantities are used to calculate the radiative decay and absorption rates (normalized by the decay rate in the absence of any structures), which in turn determine the quantum efficiency.For information on calculating the excitation field enhancement, see the Mie scattering example. Calculation of the absolute decay rates, as well as the effect of the intrinsic non-radiative decay rates, are beyond the scope of this example. Such information is typically obtained by experiment.

Run and results

Instructions for running the model and discussion of key results

Decay rates and quantum efficiency

The analysis group “quantum_efficiency”, added from the Object Library , determines the total power radiated by the dipole and by the whole structure (dipole and nanosphere antenna) in order to calculate the normalized radiative and loss decay rates, as well as the quantum efficiency in the Appendix. The results below for a dipole at a distance of d =15nm from the sphere surface show that the total decay rate (radiative + loss decay rates) is enhanced by up to ~22 times. The quantum efficiency is less than 1 due to the absorption from the nanosphere.

The script fluorescence_enhancement_comparison_theory.lsf  generates the plots below for the radiative and loss decay rates and the quantum efficiency calculated by the analysis group “quantum_efficiency” (see above), together with the analytical results for a distance d±dz, where dz=5nm is the mesh step. Note that the results are very sensitive to the position of the dipole and the nanosphere, which can only be accurately defined to within about dz in the FDTD simulation.

The simulation results are within the boundaries given by the analytical calculation; however, the current mesh step is quite coarse, which leads to a large uncertainty in the results. For more accuracy it is important to refine the mesh, as discussed in Convergence .

Important model settings

Description of important objects and settings used in this model

Mesh override region

Mesh settings have a significant impact on the simulation results due to:

• The sensitivity to the dipole-antenna distance (see Run and results ).
• The effect of the mesh size on the accuracy of the power radiated by the dipole (as explained below in Box monitor vs. dipolepower).
• “Lightning rod” effects caused by the staircasing when meshing with a rectangular grid in the Appendix).
• [Advanced] The precise behavior of the dipole source depending on its position relative to the mesh cell: For highest accuracy, the dipole must be placed halfway between two mesh points along the direction of the dipole orientation (for a more detailed discussion see Green’s Function and LDOS ). This can be achieved by adjusting the location and/or span of the mesh override region in that direction depending on the chosen mesh step. The advantage of placing the dipole correctly is that the results for coarse mesh will be more accurate and they will converge faster as the mesh is refined. The z span of the mesh override region is 105 nm so that the dipole, located at z = 35 nm, lands half way between two mesh points along the z direction.

Box monitor vs. dipolepower

You can choose between two methods to calculate the power radiated by the dipole in the “quantum efficiency” analysis group, depending on the value of the analysis variable “use dipole box”:

• With the script command dipolepower (“use dipole box” = 0): Use this is option for coarse meshes when there are not enough mesh cells between the source area of the dipole (grey area in the GUI object) and the transmission box monitor enclosing the dipole (for a more detail discussion see this example ).
• With a small transmission box monitor (“use dipole box” = 1): Use this option for fine meshes, especially if the mesh step is of the order of \( \lambda/1000 \); in this regime the dipolepower script command might return incorrect results.

Simulation auto shutoff

The simulation auto-shutoff is used to terminate simulations when the estimated electromagnetic energy in the system has decayed to a small fraction of its peak value. However, the default value of this fraction (1e-5) is not enough in this example because we have a broadband simulation and a structure where light can be potentially trapped due to resonant behavior. This early termination leads to ripples in the results from frequency domain monitors at certain wavelengths. To avoid this problem the auto-shutoff limit is decreased from its default 1e-5 to 1e-8.

Symmetry

Depending on the symmetry of the antenna structure and the configuration of the source (position and orientation of the dipole moment), symmetric/antisymmetric boundary conditions can be used to reduce simulation time and memory. In this example the dipole and the center of the nanosphere antenna are located along the z axis, and so we can use antisymmetric boundaries both along the x and y directions for a dipole oriented in the z direction.

Updating the model with your parameters

Instructions for updating the model based on your device parameters

Taking the model further

Information and tips for users that want to further customize the model

Isotropic emission

In this example we have assumed the orientation of the dipole moment of the quantum emitter is well-defined and so a specific orientation of the dipole source is used; however, in many situations the emission process is isotropic (e.g. in spherical quantum dots). This type of emission requires running three separate FDTD simulations of the same dipole oriented along the X, Y and Z axes and averaging the decay rates calculated for each orientation. This is the approach used to model incoherent unpolarized dipoles .

Optimization

There are several structure parameters that can be changed to maximize the fluorescence enhancement, for example: antenna size, material, shape, antenna-dipole distance, etc. The optimization process might be difficult to do using the built-in optimization tool , depending on the complexity of the figure of merit. There could be some advantage to use the a 3rd party driven optimization workflow , such as MATLAB.

Convergence

Tips for ensuring that your model is giving accurate results

The initial mesh settings in this example lead to results that are not very accurate. This is a reasonable first step because the simulations will take less than a minute to complete on most computers. The next step is to test the convergence of the results as the mesh is refined. This process should be done carefully due to the sensitivity of the calculation to the mesh.

For the geometry discussed in this example, the analytical calculation in the script file fluorescence_enhancement_comparison_theory.lsf provides a method to verify that the simulation results are within the uncertainty in the dipole-antenna distance as resolved by the FDTD mesh. As an example, the figure below shows the results for a finer mesh step of 2.5 nm in the mesh override region. Note that the uncertainty estimated from the analytical calculation becomes smaller but the FDTD result is still within the analytical boundaries.

Even for other more complex geometries, the general simulation setup can be optimized first by using this simple example and then replacing the nanosphere by the actual antenna structure.

As you reduce the mesh step, bear in mind that the dipole power calculation based on dipolepower (“use dipole box” = 0) will become unreliable for fine meshes and so it is necessary enable the use of the transmission box monitor (“use dipole box” = 1). In the example above there is still no significant difference between the two methods but for finer meshes, it will be important to set “use dipole box” = 1.

Additional resources

Additional documentation, examples and training material

Related publications

1. L. Novotny and B. Hecht, "Principles of Nano-Optics", Cambridge University Press, Cambridge, (2006), Chapter 8.
2. P. Bharadwaj and L. Novotny, "Spectral dependence of single molecule fluorescence enhancement," Opt. Express 15, 14266-14274 (2007).
3. J. Gersten and A. Nitzan, "Spectroscopic properties of molecules interacting with small dielectric particles", J. Chem. Phys. 75, 1139 (1981).
4. Chowdhury, et al, “Systematic Computational Study of the Effect of Silver Nanoparticle Dimers on the Coupled Emission from Nearby Fluorophores”, J. Phys. Chem. C., 112(30), 11236-11249 (2008).

Acknowledgments

We gratefully acknowledge collaborative development of this example with Joseph Lakowicz and Mustafa Chowdhury at the Center for Fluorescence Spectroscopy, Medical Biotechnology Center, University of Maryland School of Medicine, Stephen Gray at Argonne National Labs, and discussions with David Ginger at the University of Washington.

Related Ansys Innovation Courses

Appendix

Additional background information and theory

Simulating fluorescence with FDTD

While fluorescence is a quantum mechanical effect, the radiation characteristics of a quantum emitter can be considered using classical electromagnetism, assuming the quantum emitter acts as a point-like dipole with a time-varying dipole moment (see Novotny, et al.). The dipole source in FDTD can be used to study these radiation characteristics in a homogeneous medium, or in a non-uniform environment where the emission is affected by the scattering and absorption of light by structures near the dipole. The latter is possible since the dipole source is implemented as a soft source .

This section presents a discussion of the theory of fluorescence enhancement focusing on the information that can be extracted from FDTD simulations. The main results relevant for this example are summarized in the subsection Quantum efficiency.

Fluorescence enhancement and decay rates

The emission rate of a quantum emitter \( \gamma_{\text{emi}} \) in the presence of an antenna structure is determined by the excitation rate \( \gamma_{\text{exc}} \) and the quantum efficiency (or yield) of the process \( QE \). Assuming these two effects to be uncorrelated, we can write (see Bharadwaj, et al.):

$$ \gamma_{\text{emi}} = \gamma_{\text{exc}} \cdot QE $$

Using the superscript ‘0’ to denote corresponding free-space (homogeneous medium) quantity, the enhancement of the fluorescence due to the interaction with the antenna structure becomes:

$$ \frac{\gamma_{\text{emi}}}{ \gamma^{0}_{\text{emi}}} = \frac{\gamma_{\text{exc}}}{ \gamma^{0}_{\text{exc}}} \frac{QE}{QE^{0}} $$

The enhancement of the excitation rate is given by the ratio of the electric fields with and without antenna under the excitation conditions:

$$ \frac{\gamma_{\text{exc}}}{ \gamma^{0}_{\text{exc}}} = \frac{|\bf{E}|^2}{|\bf{E}^{0}|^2} $$

The Mie scattering example shows how to calculate this field enhancement.

The quantum efficiency is the ratio between the decay rate due to radiative transitions and the total decay rate. Therefore, we have:

$$ QE^0 = \frac{ \gamma^{0}_{\text{rad}}}{ \gamma^{0}_{\text{rad}} + \gamma^{0}_{\text{nr}}} \, \text{ and } \, QE = \frac{ \gamma_{\text{rad}} }{ \gamma_{\text{rad}} + \gamma_{\text{nr}} + \gamma_{\text{loss}}} $$

where:

• \( \gamma_{\text{rad}} \) and \( \gamma^0_{\text{rad}} \) are the decay rates of radiative transitions in the system with and without antenna structure, respectively. They typically correspond to the rate at which photons are collected in the far field experimentally.
• \( \gamma_{\text{nr}} \) and \( \gamma^0_{\text{nr}} \) are the decay rates of nonradiative transitions in the emitter with and without antenna structure, respectively. Nonradiative transitions include all the decay processes in the emitter that are not electromagnetic, for example, phonon excitations (heat). The presence of the antenna generally does not affect these processes and so we will assume a single intrinsic non-radiative decay rate \( \gamma_{\text{nr}} = \gamma^0_{\text{nr}} \). This decay rate cannot be calculated with FDTD.
• \( \gamma_{\text{loss}} \) is an additional nonradiative decay rate of the transitions in the emitter leading to photons that are absorbed (or otherwise lost) in the antenna when present. Typically, it can be determined by the rate of photons that cannot be collected in the far field experimentally.

Normalized decay rates

The absolute decay rates \( \gamma_{\text{rad}} \), \( \gamma^{0}_{\text{rad}} \), \( \gamma_{\text{loss}} \) cannot be calculated with FDTD. However, it can be shown that for dipole transitions that only occur through radiation, the quantum mechanical decay rate in an inhomogeneous environment \( \gamma_{\text{dip}} \)  and the classical power radiated by the dipole in the same environment \( P_{\text{dip}} \) (see Novotny, et al.) are related by:

$$ \frac{\gamma_{\text{dip}}}{\gamma^{0}_{\text{rad}}} = \frac{P_{\text{dip}}}{P^{0}_{\text{rad}}}$$

where \( \gamma^{0}_{\text{rad}} \) and \( P^{0}_{\text{rad}} \) are the decay rate and radiated power of the dipole in a homogeneous environment. \( P^{0}_{\text{rad}} \) can be calculated for a dipole source in FDTD using sourcepower . The emission rate enhancement \( \gamma_{\text{dip}} / \gamma^{0}_{\text{rad}} \) is known as the Purcell factor .

We can write analogous expressions for the normalized decay rates:

$$ \frac{ \gamma_{\text{rad}} }{ \gamma^{0}_{\text{rad}} } = \frac{ P_{\text{rad}} }{ P^0_{\text{rad}} } \, \text{ and } \, \frac{ \gamma_{\text{loss}} }{ \gamma^{0}_{\text{rad}} } = \frac{ P_{\text{abs}} }{ P^0_{\text{rad}} } $$

where:

• \( P_{\text{rad}} \) is the power radiated out of the system (emitter and antenna) that is collected in the far field. This can be calculated in FDTD simulations using a transmission box monitor enclosing the entire system.
• \( P_{\text{loss}} \) is the power absorbed by the antenna. In FDTD simulations it can be calculated directly with a transmission box monitor enclosing the antenna only. However, it is more convenient to calculate the power \( P_{\text{dip}} \) radiated by the emitter in the presence of the antenna and then use \( P_{\text{loss}} = P_{\text{dip}} - P_{\text{rad}} \). The power \( P_{\text{dip}} \) can be found using either the script command dipolepower or a transmission box monitor enclosing the dipole alone.

Quantum efficiency

We now assume that the emitter has a high intrinsic quantum efficiency and thus we can set \( QE^{0} = 1 \) and \( \gamma^0_{\text{nr}} = 0 \). Using the relationships between power and decay rates above, the quantum efficiency ratio with respect to \( QE^{0} \)  can be written as:

$$ \frac{QE}{QE^{0}} = QE = \frac{ \gamma_{\text{rad}} }{ \gamma_{\text{rad}} + \gamma_{\text{loss}} } = \frac{ P_{\text{rad}} }{ P_{\text{rad}} + P _{\text{loss}} }  $$

The analysis group “quantum_efficiency” used in this example follows this assumption and returns the quantum efficiency ratio above, as well as the following normalized decay rates:

• gamma.gamma_radiative: \( \gamma_{\text{rad}} / \gamma^{0}_{\text{rad}} = P_{\text{rad}} / P^0_{\text{rad}} \)
• gamma.gamma_loss: \( \gamma_{\text{loss}} / \gamma^{0}_{\text{rad}} = P _{\text{loss}} / P^0_{\text{rad}} \)
• gamma.gamma_dipole: \( \gamma_{\text{dip}} / \gamma^{0}_{\text{rad}} = P _{\text{dip}} / P^0_{\text{rad}} \)

Lightning rod effects

In this example, we are studying a dipole that is very close to a curved metal surface. Since the curved surface is meshed in a rectangular mesh, the sphere will be represented as a group of rectangular regions. When the structure is metal, this can lead to strong local fields near sharp points (or corners) in the meshed structure. This simulation artifact might affect the accuracy of the results if the dipole is located near one of these sharp points. As a consequence, there are some unfortunate geometries where the results will converge slowly as the mesh is refined.

In the sphere case, this situation will most likely to occur when the radius, \( a \), is an integer multiple of half the mesh step. For instance, a sharp point will arise if the mesh cell at \( (0,0,a) \) is counted as inside the metal sphere, while neighboring cells, for example \( (dx,0,a) \), are in air. To alleviate this problem, reduce the radius of the sphere by a small amount, as done in this example, where the intended radius (20 nm) was reduced by 0.01 nm.

If you are interested in testing your own structures, use an index monitor and disable the spatial interpolation under the 'Advanced' tab. You can then look at all three components of the index (x, y and z) and identify any potential sharp points near the dipole location, as illustrated below for spheres with slightly different radius. Near the dipole location, the one on the left has a sharp point while the one on the right has a flatter surface after meshing, and so the results for the latter are expected to converge faster with smaller mesh step.