Self-consistent numerical modeling of radiatively damped Lorentz oscillators

Recent progress towards realizing quantum emitters (QEs) suitable for integration in quantum information networks stimulates the demand for a self-consistent numerical approach to describe scattering of radiatively limited QEs in complex dielectric environments. As the QE response has to be characterized without the use of phenomenological damping parameters, the divergent nature of the pointlike emitter's in-phase self-field has to be carefully dealt with to avoid unphysical frequency shifts. Here we provide a solution to this problem and show two ways to obtain accurate results in the weak excitation limit using finite-difference time-domain algorithms. One of these approaches lays important groundwork needed for future simulations of nonlinear QE networks. The solution for dealing with the frequency shift reveals that dynamical contributions to the resonant depolarization field of arbitrarily small dielectric objects make crucial contributions to the net dipole moment induced by an external field when radiative scattering is the only source of damping.

Figure 1

(a) Dielectric function for a sphere of LO material given by Eq. (13), with ${\omega }_0=3.55\times {10}^{15}$ rad/s, $\mu =690$ D, $\mathrm{\Gamma }=1\times {10}^{12}$ rad/s, $R=10$ nm, and $V=4\pi R^3/3$. The inset shows the dielectric function near ${\omega }_0$, marked with a dashed line. (b) Mie scattering intensity spectra generated by sphere LO materials described by the dielectric function in (a), but for sphere radii ranging between 5 and 30 nm, as labeled on the plot. The broadest peaks in each spectra are the dipole-resonance peaks. (c) Magnitudes of the resonant frequency shifts normalized by $\mathrm{\Gamma }$, versus the sphere radii in (b). The data connected by solid, dashed, and dash-dotted lines are calculated for spheres with LO material permittivities given by Eq. (16), where $\zeta =0,1$ and $\sqrt{12/15}$, respectively. Triangle and circles markers indicate negative and positive shifts, respectively. Note the $\mathrm{\Gamma }$ normalization applied to the shift is appropriate because the absolute frequency shifts and the linewidths are derived from the second and third lines, respectively, of Eq. (6), using a common scale factor proportional to the assumed dipole transition moment squared [see Eq. (17)].

Figure 2

Finite-difference time-domain simulations of the electric field intensity spectra generated by the sphere of LO material (${\omega }_0=3.55\times {10}^{15}$ rad/s, $\mu =690$ D, $\mathrm{\Gamma }=1\times {10}^{12}$ rad/s) with a 10-nm radius and permittivity given by Eq. (16) with $V=4\pi R^3/3$, undergoing free decay in vacuum, with mesh sizes $\mathrm{\Delta }x$ labeled on the plot. The electric field is monitored 27 nm away from the sphere center. The asterisks mark the dipole-resonance peaks. The bottom spectrum is for a single-mesh cell LO with $\mathrm{\Delta }x=20$ nm and $V=\mathrm{\Delta }{x}^3$, with its electric field monitored at the LO position. Each spectrum is calculated by taking the Fourier transform of the electric field after the excitation pulse has passed.

Figure 3

FDTD simulation results for an LO material sphere freely decaying in vacuum with $R=10$ nm, ${\omega }_0=3.55\times {10}^{15}$ rad/s, $\mu =690$ D, and permittivity given by Eq. (16) with $\zeta =\sqrt{12/15}$, for different mesh sizes $\mathrm{\Delta }x$. The simulated (a) frequency shift $\mathrm{\Delta }{\omega }_{\mathrm{sim}}$ and (b) decay rate ${\mathrm{\Gamma }}_{\mathrm{sim}}$ are both normalized to $\mathrm{\Gamma }$. The volume $V$ in Eq. (16) is set to $4\pi R^3/3$. The smallest five mesh sizes (rightmost data points in each plot) correspond to the resonant peaks marked by the asterisks in the upper five spectra in Fig. 2.

Figure 4

(a) Layout used for three-dimensional FDTD simulations with background dielectric permittivity ${\epsilon }_{\mathrm{B}}$ (white), plotted in the $x-y$ plane (left) and $x-z$ plane (right). The edges of the Yee cells are indicated by the light-gray grid lines. The oscillator material (circle) is placed on a Yee-cell $E_z$ field location, which is offset by half a mesh step in $\widehat{\mathbf{z}}$ from the Yee-cell corner. The field at the $E_z$ location is directly monitored by a detector (red cross). A total-field scattered-field (TFSF) source generates plane waves polarized in $\widehat{\mathbf{z}}$ (double arrow in right schematic) traveling in the direction indicated by the black arrow, and only allows scattered fields to exist outside of it. The right schematic shows a gold half-space that spans the $x-y$ plane that is introduced in some inhomogeneous environment simulations. (b) FDTD simulation results for the free decay of a single Yee-cell LO material with $V=\mathrm{\Delta }{x}^3$ (open empty markers) and LO plugin material (closed colored markers) in vacuum with resonance frequency ${\omega }_0=3.55\times {10}^{15}$ rad/s and dipole strength $\mu =30.8$ D ($\mathrm{\Gamma }=2\times {10}^9$ rad/s from Fermi's golden rule). The simulation layout in (a) is implemented, without the gold half-space. On the left axis, the downward and upward triangles plot the magnitude of the simulated frequency shift relative to ${\omega }_0,\mathrm{\Delta }{\omega }_{\mathrm{sim}}$, normalized by $\mathrm{\Gamma }$, for LOs with $\zeta =\sqrt{12/15}$ and $\zeta =1.15$, respectively. They are plotted as a function of $\beta \left(\mathrm{\Delta }x\right)$ (bottom axis), which is the mesh size expressed as a fraction of the wavelength in ${\epsilon }_{\mathrm{B}}$ and as a function of parts-per-wavelength (PPW), $1/\beta \left(\mathrm{\Delta }x\right)$ (top axis). The simulated decay rates normalized by $\mathrm{\Gamma }$ are also plotted (circles, right axis), for the original LO material and the plugin material, using $\zeta =\sqrt{12/15}$.

Figure 5

The real part of normalized harmonic self-field found by FDTD simulations multiplied by the mesh cell volume $\mathrm{\Delta }{x}^3,{\tilde{\mathbf{E}}}_{\mathrm{CW}}({\mathbf{r}}_{\mathbf{0}},\beta )/p_0\mathrm{\Delta }{x}^3$, is plotted for simulations with ${\epsilon }_{\mathrm{B}}=1$ (circles) and ${\epsilon }_{\mathrm{B}}=1.96$ (diamonds) as a function of $\beta \left(\mathrm{\Delta }x\right)$. The real part of the normalized harmonic self-field predicted with our model ${\tilde{\mathbf{E}}}_{\mathrm{LO}}^{\mathrm{div}}\left(\beta \right)/\tilde{\mathbf{p}}\left({\omega }_0\right)\mathrm{\Delta }{x}^3$, with the best-fit values ${\zeta }_1=1.1493\pm 0.0003$ (black line) and ${\zeta }_2=1.1521\pm 0.0006$ (light line), is plotted for calculations with ${\epsilon }_{\mathrm{B}}=1$ and ${\epsilon }_{\mathrm{B}}=1.96$ (refractive index $n=1.4$), respectively.

Figure 6

(a) Relative error magnitude between our model, ${\tilde{\mathbf{E}}}_{\mathrm{LO}}^{\mathrm{div}}\left(\beta \right)/\tilde{\mathbf{p}}\left({\omega }_0\right)\mathrm{\Delta }{x}^3$, and the FDTD results, $Re[{\tilde{\mathbf{E}}}_{\mathrm{CW}}({\mathbf{r}}_{\mathbf{0}},\beta )/p_0]$, for the simulations in Fig. 5, plotted as a function of $\beta ={\omega }_0\sqrt{{\epsilon }_{\mathrm{B}}}\mathrm{\Delta }x/\left(2\pi c\right)$. Best-fit values ${\zeta }_1=1.1493$ and ${\zeta }_2=1.1498$, applied to $f\left(\mathrm{\Delta }x\right)$ [Eq. (20)] for simulations with ${\epsilon }_{\mathrm{B}}=1$ (circles, solid black line) and ${\epsilon }_{\mathrm{B}}=1.96$ (diamonds, solid colored line), respectively. The red dashed line shows the relative error between the real part of the field amplitude averaged over one mesh cell ($\zeta =1$) and $Re[{\tilde{\mathbf{E}}}_{\mathrm{CW}}({\mathbf{r}}_{\mathbf{0}},\beta )/p_0]$ for vacuum simulations (${\epsilon }_{\mathrm{B}}=1.96$ results are not shown, as they are approximately the same). The shaded region indicates where the relative error causes frequency shifts of less than a linewidth in LO simulations. (b) Relative error magnitude between the predicted imaginary contribution, ${\left(k_{\mathrm{B}}\right)}^3/\left(6\pi {\epsilon }_0{\epsilon }_{\mathrm{B}}\right)$, and the FDTD results, $Im[{\tilde{\mathbf{E}}}_{\mathrm{CW}}({\mathbf{r}}_{\mathbf{0}},\beta )/p_0]$, for simulations with ${\epsilon }_{\mathrm{B}}=1$ (circles, solid black line) and ${\epsilon }_{\mathrm{B}}=1.96$ (diamonds, solid colored line).

Figure 7

The electric field intensity spectrum of a Lorentz oscillator (LO) in free decay, located 9 nm off the surface of a gold half-space is plotted (shaded region). The LO is simulated with FDTD methods using a plugin material with resonance frequency ${\omega }_0=2.51\times {10}^{15}$ rad/s (${\lambda }_0=750$ nm) and dipole strength $\mu =51.9$ D ($\mathrm{\Gamma }=2\times {10}^9$ rad/s) and $\mathrm{\Delta }x=$ 2 nm, that applies ${\mathbf{E}}_{\mathrm{LO}}^{\mathrm{div}}\left(t\right)$ to correct for the divergent self-field contribution. The simulated spectrum is fit by a Lorentzian line shape (dotted line).