Developing a time-domain method for simulating statistical behaviour of many-emitter systems in the presence of electromagnetic field

TL;DR

This paper introduces a novel time-domain numerical method that incorporates probabilistic molecular transitions into classical electrodynamics to simulate the statistical behavior of many-emitter systems interacting with electromagnetic fields.

Contribution

The paper develops a new computational approach that integrates quantum-like transition probabilities into classical Maxwell's equations for accurate statistical emission modeling.

Findings

Successfully validated against experimental data

Accurately predicts emission directionality and coherence

Demonstrates applicability to complex nano-optical systems

Abstract

In this paper, one of the major shortcomings of the conventional numerical approaches is alleviated by introducing the probabilistic nature of molecular transitions into the framework of classical computational electrodynamics. The main aim is to develop a numerical method, which is capable of capturing the statistical attributes caused by the interactions between a group of spontaneous as well as stimulated emitters and the surrounding electromagnetic field. The electromagnetic field is governed by classical Maxwell's equations, while energy is absorbed from and emitted to the (surrounding) field according to the transitions occurring for the emitters, which are governed by time-dependent probability functions. These probabilities are principally consistent with quantum mechanics. In order to validate the proposed method, it is applied to three different test-cases; directionality of fluorescent emission in a corrugated single-hole gold nano-disk, spatial and temporal coherence of fluorescent emission in a hybrid photonic-plasmonic crystal, and stimulated emission of a core-shell SPASER. The results are shown to be closely comparable to the experimental results reported in the literature.

Tables

Table 1. Table 1: Optical properties of emitters used in test-cases.

	excitation wavelength	emission wavelength	quantum yield	excited state life-time	${\tau }_{10r}$	${\tau }_{10nr}$	${\tau }_{21nr}$	${\tau }_{20nr}$
Alexa Fluor 647	650 nm	672 nm	0.3	1.04 ns	3.5 ns	1.7 ns	0.01 ns	11.4 ns
Rhodamine 6G	525 nm	550 nm	0.95	4.08 ns	4.3 ns	92.7 ns	0.01 ns	824.1 ns
Sulfohodamine 101	575 nm	591 nm	0.8	4.2 ns	5.2 ns	23.8 ns	0.01 ns	257.5 ns
Oregon Green 488	494 nm	524 nm	0.9	4.3 ns	4.8 ns	48.8 ns	0.01 ns	265.9 ns

Equations

$$\ddot{p}+\gamma \dot{p}+{\omega }_a^{2}p=(e^2/m)\mathbf{E}(t)\cdot \hat{\mathbf{p}}.$$

$$P_{02}(t_n)=\exp \left(-{\left(p(t_n)-p_{max}\right)}^2/2{\sigma }^2\right),$$

$$p(\omega )=\frac{(e^2/m)|E_0|}{\sqrt{{\left({\omega }_a^2-{\omega }^2\right)}^2+{\left(\gamma \omega \right)}^2}}.$$

$$\Delta p=p_{max}-p({\omega }_a-\frac{\Delta {\omega }_a}{2}).$$

$$\epsilon \dot{\mathbf{E}}=\nabla \times \mathbf{H}-\sum _{i=1}^N{\dot{\mathbf{p}}}_i.$$

$$P_{21nr}(\tilde{t})=$$

$$P_{20nr}(\tilde{t})=$$

$$P_{10nr}(\tilde{t})=$$

$$P_{10r}(\tilde{t})=$$

$$\mathbf{F}(\tilde{t})={\mathbf{F}}_0\exp \left(-\frac{(\tilde{t}-t_0{)}^2}{{\sigma }_e^2}\right)\sin ({\omega }_e(\tilde{t}-t_0)),$$

$$A={\tau }_{21nr}/({\tau }_{21nr}+{\tau }_{20nr}),$$

$$\Gamma (\tau )=⟨u(t+\tau )u^{\ast }(t)⟩,$$

$${\tau }_c={\int }_{-\infty }^{\infty }|\gamma (\tau ){|}^{2}d\tau .$$

$${\mathcal{V}}_{ij}=\frac{2\sqrt{u_i^2(t)u_j^2(t)}}{u_i^2(t)+u_j^2(t)}{M}_{ij}.$$

$$M_{ij}=\frac{⟨u_i(t)u_j^{\ast }(t)⟩}{\sqrt{⟨u_i^2(t)⟩⟨u_j^2(t)⟩}}$$

$$P_{10stim}(t_n)=\exp \left(-{\left(p_{stim}(t_n)-p_{stim,max}\right)}^2/2{\sigma }_{stim}^2\right),$$

$${\ddot{p}}_{stim}+{\gamma }_{stim}{\dot{p}}_{stim}+{\omega }_e^2{p}_{stim}=(e^2/m)|\mathbf{E}(t)|.$$

$$\begin{array}{rl}{\mathbf{F}}_{stim}(\tilde{t})={\mathbf{F}}_{stim,0}& \exp \left(-\frac{(\tilde{t}-t_{stim,0}{)}^2}{{\sigma }_{stim,e}^2}\right)\\ \times & \sin ({\omega }_e(\tilde{t}-t_{stim,0}))\end{array}$$

Full text

I Introduction

Since early developments in the modern optics, light-matter interaction Frisk Kockum et al. (2019) has been the central issue in numerous applications, ranging from lasers to photonic computers. In this field, a vast number of researches has been triggered by Purcell’s statement Purcell (1995); Pelton (2015) that the behaviour of a photon emitter is highly sensitive to the surrounding electromagnetic (EM) field; both the decay rate  Geddes and Lakowicz (2005) and the emission power Russell et al. (2012) of the emitter can be effectively controlled by adjusting the EM field. Such adjustment can be performed using photonic crystal cavities Englund et al. (2005); Koenderink et al. (2005); Zhang et al. (2018), metallic surfaces and nano-particles Shimizu et al. (2002); Ribeiro et al. (2017), plasmonic structures Kinkhabwala et al. (2009); Shi et al. (2014a) and hybrid photonic-plasmonic structures (HPPSs) Shi et al. (2014b). Another fact to consider is that the emitted light itself can also affect the neighbouring luminescent molecules either directly, as a short-range interaction, or by exciting photonic and/or plasmonic guided modes, which leads to a long range interaction. In this sense, one needs to deal with complicated time-varying internal and external interactions in a many-body problem Richter et al. (2015); Agarwal (1971). Therefore, further investigations are demanded to develop efficient means to control the response and enhance the emission to a higher degree for applied purposes Ribeiro et al. (2017).

Considering the practical limitations of the delicate experimental setups, a rather low-cost numerical method can be effectively employed to provide guidelines for experiments Bradford and Shen (2014); Lee and Mycek (2018) while it also advances our fundamental understanding of physical phenomena Cao et al. (2000); Kivisaari et al. (2018). However, the presence of numerous emitters in the vicinity of a nano-structure results in special collective attributes, e.g. coherence, which requires developing an approach to numerically capture the statistical physics of the phenomena. It is necessary to ensure that the probabilistic nature of molecular transitions is taken into account and consequently, one can utilize a statistical approach to infer the collective attributes. In other word, any two emitters of the ensemble that are in the same environmental condition do not deterministically behave in the same way. In addition, the numerical method should be capable of handling the interaction between emitters and the external EM field as well as the emitter-photon-emitter interactions.

Emitters can be numerically simulated using two different class of methods; those developed based on a macroscopic viewpoint, i.e. using statistically averaged quantities, like effective optical parameters, e.g. permittivity, conductivity, and wave number Cao et al. (2000); Hagness et al. (1996), or population densities for molecular energy levels Chang and Taflove (2004); Nagra and York (1998); Ziolkowski et al. (1995). In the second class of methods, a microscopic viewpoint is adopted, i.e. one focuses on the behaviour of a single molecule using quantum mechanical methods, like solving the atomic master equation, density function equation, or Schrödinger’s equation Dzsotjan et al. (2010); Dum et al. (1992); Hong et al. (2015); Chen et al. (2010), or modeling the molecular dipole moment of emitters in a classical manner using a damped driven harmonic oscillator differential equation Englund et al. (2005); Genevet et al. (2010); Wang et al. (2011); Bauch and Dostalek (2013); Hoang et al. (2015); Javadi et al. (2018). The latter approach has been frequently used in the literature due to its simplicity and capability to be implemented within the framework of conventional numerical methods developed for EM wave propagation, e.g. Finite-Difference Time-Domain (FDTD) Taflove and Hagness (2005).

The macroscopic viewpoint is well qualified for simulating a bulk of emitting matters; as an example, Chang and Taflove Chang and Taflove (2004) successfully simulated laser gain material by developing a semi-classical approach. Nevertheless, in this viewpoint, the behaviour of numerous microscopic emitters, that are subject to identical conditions, is represented by the averaged macroscopic quantities. Therefore, despite its promising performance for a bulk of active matter, this viewpoint is not a suitable choice for cases of many emitters with individually different dipole orientations and local (especially time-dependent) heterogeneity at the molecular-scale. Moreover, in this class of approaches, the probabilistic nature of problems is lost as a result of the deterministic governing equations.

On the other hand, numerical methods are developed adopting microscopic viewpoint with the quantum mechanical approach Chen et al. (2010). In these approaches, due to a prohibitive computational cost, the governing equation, e.g. Schrödinger’s equation, including the terms corresponding to the external field-emitter and emitter-emitter interactions, can be solved only for a rather small number of emitters. This shortcoming can be resolved by considering all emitters as identical to each other Richter et al. (2015). However, this requires the environmental conditions to be completely the same for all the emitters. Moreover, acquiring the microscopic viewpoint with a classical approach for modeling the dipoles, one cannot make any distinction between the excitation and emission frequency in cases of three level emitters like fluorescent molecules Witthaut and Sørensen (2010); Roy (2011). More importantly, representing the dipole moment of distinct emitters using a single deterministic equation for harmonic oscillator, the statistical nature is lost.

The subject of the present work is to show how the probabilistic nature of the molecular transitions can be microscopically taken into account while the statistical attributes of a rather large set of emitters are macroscopically calculated. To this end, a semi-classical approach is proposed, in which transition probabilities are introduced in order to handle the behaviour of emitters, which are considered as single molecules with particular behaviour and distributed in the computational cells. To the best of our knowledge, this is the first time that such approach has been acquired. It is worth to note that since the presented method is based on generating random samples from probability distributions in a repeated manner, it can be categorized within the broad class of Monte Carlo methods. This class of methods has shown a promising performance in handling the probabilistic nature of many-body problems Vishwanath et al. (2002); Lee and Mycek (2018). Here, the introduced algorithm is implemented within the framework of the FDTD method and validated against previously reported experimental results for spontaneous and stimulated emissions; directionality and coherence of spontaneous emission and stimulated emission of a SPASER. It is worth noting that in the present work, the implemented algorithm is developed for a three-level fluorescent molecule as the emitter. However, it is straightforward to generalize the algorithm in order to include two- and four-level emitters as well.

II Numerical modeling

An emitter with three energy bands, e.g. a fluorescent molecule can be modeled as a three (energy) level system, for which absorption and both radiative and non-radiative transitions occur Jablonski (1933); Albani (2008) (Fig. 1). From the quantum mechanical viewpoint, a time dependent probability function can be associated to each of these possible transitions.

In the following, the implementation of these probability functions within the framework of a time-domain method, which is originally developed for simulating EM field propagation, is briefly described for a many-emitter system. For more convenience in this paper, the subscripts are corresponding to the energy levels as numbered in Fig. 1. Moreover, the energy levels of each molecule are identified using three occupation numbers $n_{0}$, $n_{1}$, and $n_{2}$, where at each time-step, only one of them is 1 while the other two are 0. This is due to the fact that at any instance of time, only one of the energy levels can be occupied by the molecule. Here, $t$ is the total time of the simulation, while $\tilde{t}$ is the time passed from the last transition occurred for the molecule, $P_{ij}$ is the time-dependent probability of transition from level $i$ to level $j$, and the subscripts $r$ and $nr$ refer to the radiative and non-radiative transitions, respectively.

At each time-step, the state is checked for each molecule in the system according to these criteria:

• •

For a molecule in the ground state ($n_{0}=1$), it is only possible to transit into the second excited state ($0\overset{absorb}{\rightarrow}2$) by absorbing an enough amount of energy. Computationally, this occurs if the randomly generated number $0<r<1$ is less than (or equal to) the corresponding probability $P_{02}(t)$.

• •

For a molecule in the second excited state ($n_{2}=1$), the non-radiative decay is possible to either the first state ($2\overset{nr}{\rightarrow}1$) or the ground state ($2\overset{nr}{\rightarrow}0$). If random number $r$ is less than (or equal to) $P_{21nr}(\tilde{t})$ transition $2\overset{nr}{\rightarrow}1$ occurs, and on the other hand, the molecule experiences $2\overset{nr}{\rightarrow}0$ if $P_{21nr}(\tilde{t})<r\leq(P_{21nr}(\tilde{t})+P_{20nr}(\tilde{t}))$.

• •

For a molecule in the first excited state ($n_{1}=1$), both the radiative and non-radiative decays to the ground state are possible. Transition $1\overset{nr}{\rightarrow}0$ occurs if $r\leq P_{10nr}(\tilde{t})$. Else, if $P_{10nr}(\tilde{t})<r\leq(P_{10nr}(\tilde{t})+P_{10r}(\tilde{t}))$, the radiative transition of $1\overset{r}{\rightarrow}0$ takes place and consequently, a wave packet with the central frequency equal to emission frequency of the fluorophore $\omega_{e}$ is emitted from a point source at the position of the emitting molecule.

Once a transition occurs, the occupation numbers $n_{0}$, $n_{1}$, and $n_{2}$ are correspondingly reset, e.g., $0\overset{absorb}{\rightarrow}2$ is associated with resetting numbers as $n_{0}=0$ and $n_{2}=1$. During emission, the re-excitation of the corresponding molecule is prevented by setting its emission flag to on, which means the state of the molecule is ignored and the emission continues until the emitted wave packet vanishes. At this moment, the emission flag is off and the molecule is at its ground state.

The probabilities associated with the above-mentioned transitions between these levels are estimated as described in the following. It is evident that these probabilities can be modified into more sophisticated functions that are obtained by quantum mechanical analysis of molecules. However, the following formulas show the simplest probability functions that satisfy the physical requirements and work successfully for the test-cases solved in this work.

Absorption: For a molecule with the frequency of the maximum absorption, $\omega_{a}$, the absorption of energy from the field and excitation of the molecule becomes more probable if the electric field imposed at the position of molecule, $\mathbf{E}(t)$ incorporates a frequency component tending to $\omega_{a}$. In order to examine this possibility, the occurrence of the resonance of a harmonic oscillator driven by $\mathbf{E}(t)$ is checked. This oscillator resembles the electric dipole moment of the molecule, $\mathbf{p}$. Within the context of the time-domain method used in the present work, i.e. FDTD, the equation that governs the response of the harmonic oscillator is considered as an auxiliary differential equation (ADE) Kashiwa and Fukai (1990). It must be highlighted that using the aforementioned ADE is the most efficient way to check the frequency components of $\mathbf{E}(t)$ against $\omega_{a}$. The implemented ADE is

[TABLE]

Here, $\gamma$ is a damping factor and $e$ and $m$ are the charge and mass of electron, respectively. Vector $\hat{\mathbf{p}}$ is a unit vector representing the orientation of dipole moment of the emitter, i.e. $\mathbf{p}=p\hat{\mathbf{p}}$. The magnitude of the dipole moment is initially zero and updated (at time step $n+1$) using a second-order central time-marching scheme as done for Maxwell’s equations in the adopted FDTD method Taflove and Hagness (2005). At the onset of resonance, the amplitude of $\mathbf{p}$ tends to its maximum value $p_{max}$ and consequently, the transition from level 0 to 2 becomes more probable. Therefore, the corresponding probability is estimated as

[TABLE]

where $t_{n}$ represents time at $n$th time-step, i.e $t_{n}=n\Delta t$, and $\sigma$ is determined in terms of the absorption bandwidth of the fluorophore, $\Delta\omega_{a}$. In order to derive an equation for $\sigma$, one can consider the external field in its simplest form $\mathbf{E}(t)=\mathbf{E_{0}}\sin(\omega t)$ and analytically solve Eq. (1) for the amplitude of $\mathbf{p}$, which depends on $\omega$ as

[TABLE]

As illustrated in Fig. (2), for an absorption bandwidth of $\Delta\omega_{a}$, an amplitude interval of $\Delta p$ is defined as

[TABLE]

In this way, the standard deviation $\sigma$ is considered equal to $\Delta p$.

It must be noted that since $p_{max}\propto\frac{1}{\gamma\omega_{a}}$, both $p_{max}$ and damping factor $\gamma$ are arbitrary factors which should be determined correspondingly.

Upon the occurrence of resonance, energy is also absorbed from the external EM field, which is governed by Ampere-Maxwell’s relation

[TABLE]

Here, $N$ denotes the number of fluorophores associated with the computational grid cell, which is proportional to the density of the fluorescent material at the same point within the FDTD discretized domain. It is possible to have various fluorescent densities at different locations of the structure. If referring to Eq. 2 transition $0\overset{absorb}{\rightarrow}2$ occurs, $p$ and its first temporal derivative are set to zero. On the other hand, if $0\overset{absorb}{\rightarrow}2$ does not occur, the dipole moment of the molecule is updated via equation 1.

Non-radiative transitions: For the present model, three non-radiative transitions are considered; $2\overset{nr}{\rightarrow}1$, $2\overset{nr}{\rightarrow}0$ and $1\overset{nr}{\rightarrow}0$. The corresponding probabilities are estimated as

[TABLE]

Here, $\tau_{ijnr}$ represents the time-constant for a non-radiative decay between levels $i$ and $j$. The asymptotic behavior of these functions guaranties a definite decay at an infinitely long time.

Radiative transition: Radiative transition is only considered as a decay from level 1 to level 0, for which the corresponding probability is estimated as

[TABLE]

where $\tau_{10r}$ is the time constant of the radiative decay. Once the transition occurs, a wave-packet is emitted with a central-frequency equal to the emission frequency of the fluorophore, $\omega_{e}$. This is implemented as a soft point source, which is a vector function, $\mathbf{F}(\tilde{t})$, that is superimposed to the electric (and/or magnetic) field at the position of the molecule. This function is formulated as

[TABLE]

where $\mathbf{F}_{0}=F_{0}\hat{\mathbf{n}}$. In this equation, $\hat{\mathbf{n}}$ is a randomly oriented unit vector and $F_{0}$ is determined in a manner that the total energy of the wave-packet is equal to the energy of an emitted photon. Here, $\sigma_{e}$ is associated with the band-width of the emission and $t_{0}\approx 3\sigma_{e}$ is a time-offset, which guarantees that the at $\tilde{t}=0$, $\mathbf{F}$ approaches zero.

The normalization factors $A,\,B,\,C,\,D$ are calculated based on two physical concepts; first, it is impossible for a molecule to permanently stay at an excited state and consequently, $A+B=1$ and $C+D=1$. Second, the probability ratio of the transitions is inversely related to the corresponding decay times, $\tau_{ij}$. Therefore, one can obtain

[TABLE]

The flowcharts illustrated in Figures 3 and 4, present the implemented algorithm in more details.

In this work, the propagation of EM fields in a three-dimensional domain and successive time-steps is simulated using an FDTD package, which is developed in C++ language. Convolutional perfectly matched layer (CPML) Roden and Gedney (2000); Berenger (2002) boundary condition is used for reducing the reflection from exterior boundaries of the simulation domain. Whenever needed, the total-field/scattered-field technique and the Drude-Lorentz model are implemented to handle the incident field and dispersive materials, respectively Taflove and Hagness (2005). In all simulations, the electric and magnetic fields are initially set to zero. It is worth noting that the implementation of the same algorithm is also possible for other numerical methods, e.g. the finite-element time-domain method.

III Applications for many spontaneous emitters

In the following, the proposed method is applied to two phenomena observed for many-emitter systems; directionality of fluorescence in the presence of a plasmonic nano-structure Aouani et al. (2011) and fluorescence coherence obtained by utilizing an HPPS Shi et al. (2014b). The present method is validated by comparing the numerical results with the corresponding experimental results reported in the literature.

Here, the spatial discretization of the solution domain is set considering the criteria for minimizing the dispersion error of the FDTD method Taflove and Hagness (2005) while resolving all structural details. Governing equations are solved for a three-dimensional Cartesian mesh with Yee cells of $\Delta x=\Delta y=\Delta z$. The time-step is set according to the Courant condition. The physical properties of the specific emitters simulated in the present work are shown in Table 1.

III.1 Directional spontaneous emission

In this section, the central hole of a gold nano-disk with two concentric grooves is filled by excited fluorescent molecules as shown in Fig. 5. The capability of this plasmonic system in producing directional emission has been experimentally studied by Aouani et al. Aouani et al. (2011). The geometrical parameters (see Fig. 5), as well as the physical properties of the base structure (gold) and the fluorescent (Alexa Fluor 647 and Rhodamine 6G) molecules, are set the same as those reported in Aouani et al. (2011).

The properties of Alexa Fluor 647 and Rhodamine 6G are set according to Buschmann et al. (2003) and Magde et al. (2002), respectively (see Table 1). In this test-case, the caution is that the resulted emission distribution is hardly distinguishable due to the masking effect of the incident pump beam. Here, in order to eliminate the need for a post-processing procedure, the fluorophore molecules are considered to be initially-excited ($n_{2}=1$ for all molecules). In this way, there is no need to impose a pump beam.

Here, two different simulations are separately done one for Alexa Fluor 647 and one for Rhodamine 6G molecules as the fluorescent molecules, while the output emission is detected at two perpendicular arc-ports on the hemisphere encircling the disk as illustrated in Fig. 5. The long-time average, as well as the ensemble average of the intensity are presented in Fig. 6. The results are in a good agreement with those reported in Aouani et al. (2011) (see Fig. 6, dashed lines), i.e., for Alexa Fluor 647 the peak intensity is observed at a polar angle of around 27 degrees while for Rhodamine 6G the emission become concentrated at the zero polar angle.

III.2 Fluorescence coherence

The main aim of developing the proposed method is to capture the statistical attributes of a many-emitter system. In this sense, the method is utilized to simulate the response of a set of fluorescent molecules positioned adjacent to a photonic crystal (PC) constructed by triangular arrangement of polystyrene spheres of 500 nm diameter. This PC is placed on top of a 200 nm thick silver slab to form an HPPS as shown in Fig. 7. The incident pump ($\lambda=532$ nm) works continuously in $y$-direction while the output emission is detected in $z$-direction. The capability of this structure in producing a coherent light from fluorescent spontaneous emissions was previously reported by Shi et al. Shi et al. (2014b). Coherence of a fluorescent light is one of the statistical phenomena, which to the best of authors’ knowledge has not yet been numerically addressed. In this case, the fluorescent material, fluorophore-doped polyvinil alcohol (PVA), forms a 50 nm thick layer on top of the HPPS and also fills the vacancy between spheres (see the cross-sectional view in Fig. 7). In order to keep the numerical test-case substantially similar to the reported experimental setup, properties of the fluorophores, i.e. excitation and emission wavelength are set according to the physical properties of Sulforhodamine 101 (S101) as $\lambda_{ex}=$ 575 nm and $\lambda_{em}=$ 591 nm, respectively (Table 1). A $y$-directed 532 nm continuous-wave is used to pump the structure from the center of $x$-$z$ plane.

In Fig. 8, the time-evolution of the field distribution is plotted on a surface passing through the center of the polystyrene spheres in the $x$-$y$ plane. Since the low-intensity fluorescent emissions are masked by the pump intensity that would result in a non-clear field distribution, here, the results are plotted merely for a pulse-train. It is worth noting that for all other simulations of the current test-case, the continuous wave is used. It is observed that the first pump pulse passes through the domain without any fluorescent molecule emission (Fig. 8b). The fluorescent emission is seen in Figs. 8c-h since the molecules have gained the excitation energy required for emission. The honey-comb like pattern (seen more clearly in Fig. 8g and h) is caused by the fluorescent emission of molecules filling the vacancy between the polystyrene spheres.

In order to estimate the degree of temporal coherence, two different approaches has been employed; in one approach, the temporal coherence function (TCF) is utilized, which is the auto correlation of the signal,

[TABLE]

where angle brackets and superscript $*$ denote the time averaging and complex conjugate, respectively Goodman (2015). Here, $u(t)$ is the electric field. The degree of coherence is conclusively determined as $\gamma(\tau)={\Gamma(\tau)}/{\Gamma(0)}$ and thus, the coherence time is

[TABLE]

Using this approach for the present test-case, the coherence time is $\tau_{c}=1.14\times 10^{-13}$ sec, this is approximately equivalent to the wavelength bandwidth of $|\Delta\lambda|=6.9$ nm. In the other approach, the desirable wavelength spectrum is obtained using the Fourier transform of the time-varying electric field. For the present test-case, this spectrum is shown in Fig.9. By fitting a Gaussian function into the figure, the emission bandwidth of the peak is calculated to be approximately $|\Delta\lambda|=7.1$ nm.

It must be noted that in the present work, a perfectly ordered structure is modeled while any experimental test is subject to fabrication defects. These structural defects lead to a wider bandwidth and therefore, a lower temporal coherence can be detected for the experimental setup. This is due to the fact that any defect/disorder perturbs the Hamiltonian of the structure, which in turn broadens the band of energy Mukherjee and Gordon (2012). Considering this issue, there is a good agreement between the result obtained using the proposed numerical method and the $|\Delta\lambda|=9$ nm bandwidth reported in the literature Shi et al. (2014b).

On the other hand, the most reliable and extremely practical approach to the estimation of spatial coherence is the Young’s double slits technique Martienssen and Spiller (1964) for which, the interference fringe visibility is analytically calculated as a function of the separation distance between slits as Goodman (2015)

[TABLE]

In this equation,

[TABLE]

is the mutual coherence function between electric fields $u_{i}(t)$ and $u_{j}(t)$ detected at two distinct points $i$ and $j$ placed on a plane perpendicular to the direction of detection, which represent the positions of the slits. For the present test-case, the fringe visibility is plotted in Fig. 10. It is evident that visibility significantly decreases as the slits separation distance increases beyond $4\>\mu$m and practically, no interference pattern can be observed for a separation distance of larger than $8\mu$m. This is in close agreement with the result reported in the reference Shi et al. (2014b); clearly visible interference fringes were observed for the slit-separation distance of $4\mu$m and less clear interference fringes were recognized for the slit-separation distance of $7\mu$m. Nevertheless, it is noteworthy that using numerical simulation, it can be observed (Fig.10) how the fringe visibility varies for a wide range of separation distances. In Fig. 10 (inset), the maxima and minima of the visibility calculated using the proposed numerical method are plotted. In this figure, the intensity of modes as a function of propagation length as reported by Shi et al. Shi et al. (2014b) is also illustrated for comparison.

IV Extending the proposed method

One of the main advantages of the proposed method is its capability to be modified for other applications beyond many-spontaneous-emitters by modifying the probability functions. In this section, the proposed method is modified to capture the statistical behaviour of a many-emitter system, which also exhibits stimulated transitions besides the spontaneous transitions. To this end, it is needed to modify the proposed algorithm in order to also include the probability function $P_{10stim}$ required to model the stimulated emission transition. This probability should be based on transition-field coupling and thus can be defined analogous to Eq. 2 as

[TABLE]

where $p_{stim}$ is obtained using the following harmonic oscillator equation

[TABLE]

Here, the subscript $stim$ corresponds to the stimulated transition and $\omega_{e}$ is the frequency of the stimulating photons. It must be noted that upon transition a wave-packet of the form

[TABLE]

is emitted. Since the polarization of a photon released in stimulated emission must be aligned with the polarization of the stimulating photon, the amplitude of the wave packet is $\mathbf{F}_{stim,0}=F_{0}\hat{\mathbf{E}}(t)$, where $\hat{\mathbf{E}}(t)$ represents the unit vector aligned with the electric field at the position of the molecule and the instance of stimulated transition. It is worth noting that the previously proposed semi-classical FDTD approaches (for example see Chang and Taflove (2004) and Ziolkowski et al. (1995)) incorporate the statistically averaged quantities from a deterministic viewpoint. Therefore, those methods are only capable of modeling a bulk of emitting material, in contrast to many-emitter systems with localized sources that are successfully simulated using the present method.

In order to verify the performance of the proposed method, a core-shell type surface plasmon amplification by stimulated emission of radiation (SPASER) is simulated and the numerical results are compared to the experimental results reported by Noginov et al. Noginov et al. (2009). The structure is schematically shown in Fig. 11; the core is made of gold with the diameter of 14nm, which is enclosed by an Oregon Green 488 (see Table 1) doped silica shell of 44nm diameter. Here, the goal is only to investigate the performance of the method in simulating stimulated emission and hence, the pump mechanism is considered to be out of scope. In this sense, an initial energy-level population inversion is considered for the dye-molecules, i.e. all Oregon Green 488 molecules are initially at level 2. The emission spectrum is shown in Fig. 12, which shows two peaks around $\lambda_{spon}=520$ nm corresponding to the spontaneous emission and $\lambda_{stim}=540$ nm corresponding to the stimulated emission. Here, the peak wavelengths are calculated by matching two Lorentzian functions to the data as illustrated in Fig. 12. The agreement between these results and those reported in Noginov et al. (2009), i.e. $\lambda_{spon}=520\pm 20$ nm and $\lambda_{stim}=531$ nm, approves the validity of the proposed method for many-spontaneous/stimulated-emitters (see the inset of Fig. 12). However, the slight deviation in $\lambda_{stim}$ is caused by the so-called stair-case error generally occurs in the FDTD method. In the experimental case, the stimulated emission is resulted from the feedback provided by the surface plasmons of the core inside a cavity of a spherical shape. In the FDTD simulation, on the other hand, the spheres cannot exactly be fitted by the Cartesian grid. Here, mapping of the faces onto the computational cells leads to the appearance of some stairs at the surface of the spheres. Thus, the resonance modes of the numerically modeled simulated sphere are slightly different from those of a real spherical cavity and consequently, the peak of the stimulated emission obtained in the simulation is slightly different than that of the experiment. In order to alleviate this error one needs to either use a highly refined domain discretization that needs prohibitively intense computations or utilize a body-fitted grid, which is out of the scope of the present work.

V Conclusion

A numerical method was developed, which is capable of capturing the statistical nature of the interactions between a group of emitters and the surrounding EM field. The proposed algorithm was devised in a manner that is basically consistent with physical principles. The method has been validated for three different many-emitter systems; in the first test-case, the directionality of the fluorescent emission is captured. In the second test-case, the spatial and temporal coherence are simulated for an ensemble of spontaneous emitters. This important statistical attribute can be numerically captured merely by utilizing the proposed probabilistic approach. In order to show the capability of the proposed method beyond the spontaneous emitters, in the last test-case, the stimulated emission of a SPASER is simulated.

It is evident that the applications of the proposed method are not limited to the test-cases simulated in this work; the method can be utilized to capture other statistical attributes, e.g. transition rates. In addition, by using different probability functions and/or including two or four energy-levels, the method can be further developed for numerous cases with various types of emitters. Acquiring intensity-dependent probability functions is also a means to obtain more realistic physics. Moreover, nonlinear effects, e.g. two-photon absorption or up-conversion, are possible to be captured by adding virtual states and considering corresponding transitions in the algorithm.

It is worth noting that using the proposed method, the run-times are increased by less than $5\%$ comparing to a sole electromagnetic solver; however, the computational costs (memory and run-time) depend on the number of emitters. Considering the complexities and costs associated with experiments, the proposed method is a promising means to facilitate new designs for optical structures and advance the fundamental understanding of the statistical phenomena in the field of light-matter interaction.