Lorentz Oscillators#

The optical interactions between light and matter at the nanometer scale are important phenomena for a variety of research fields, and a rigorous understanding of these interactions requires the use of quantum electrodynamics (QED) theory. However, nanometer-scale structures are often too complex to be solved rigorously using only QED; in these cases, a classical approach that invokes QED results in a phenomenological way can be applied.

PyCharge uses the Lorentz oscillator (LO) model, an approximation derived from quantum theory (e.g., from the time-dependent Schrödinger), to simulate the interaction of a bound charge (such as an electron) with light.

In the classical model, an oscillating dipole produces electromagnetic radiation that dissipates energy and modifies the self-consistent dipole moment. The recoil force, \(\mathbf{F}_\mathrm{r}\), acting on the accelerating point charges in the dipole is called the radiation reaction or radiation damping force. The equation of motion for an undriven LO (e.g., in vacuum) that includes the radiation reaction force is given by

(1)#\[m \mathbf{\ddot r_{\rm dip}}(t) + \omega_0^2 m \mathbf{r_\mathrm{dip}}(t) = \mathbf{F}_\mathrm{r}(t),\]

where \(\mathbf{r}_{\rm dip}\) is the displacement from the LO’s negative charge to its positive charge, \(\mathbf{\ddot r_{\rm dip}}\) is its second derivative with respect to time, \(\omega_0\) is its natural angular frequency, and \(m\) is the effective mass (also called reduced masss) of the LO, defined as

(2)#\[m=\frac{m_1 m_2}{m_1+m_2},\]

where \(m_1\) and \(m_2\) are the masses of the two point charges in the dipole.

The radiation reaction force \(\mathbf{F}_\mathrm{r}\) acting on the accelerating point charges is described by the Abraham–Lorentz formula for non-relativistic velocities:

(3)#\[\mathbf{F}_\mathrm{r}(t)=\frac{q^2}{6\pi\epsilon_0 c^3}\,\mathbf{\dddot r_\mathrm{dip}}(t),\]

where \(\mathbf{\dddot r_\mathrm{dip}}\) is the third derivative of the displacement between the two charges. We can approximate \(\mathbf{\dddot r}_\mathrm{dip}\approx -\omega_0^2\mathbf{\dot r}_\mathrm{dip}\) in the above equation if the damping introduced by radiation reaction is negligible (i.e., \(|\mathbf{F}_\mathrm{r}| \ll \omega_0^2 m |\mathbf{r}_\mathrm{dip}|\)), satisfying

(4)#\[\frac{q^2 \omega_0}{m} \ll 6\pi \epsilon_0 c^3.\]

In an inhomogeneous environment, an oscillating electric dipole experiences an external electric field \(\mathbf{E}_\mathrm{d}\) as a driving force: the component of the total electric field along the polarization direction at the dipole’s origin (center of mass) \(\mathbf{R}\), generated by other sources in the system and its own scattered field.

If the above condition is satisfied, the equation of motion for a driven LO is

(5)#\[\mathbf{\ddot d}(t) +\gamma_{0} \mathbf{\dot d}(t) +\omega_{0}^{2} {\bf d}(t) = \frac{q^2}{m} \mathbf{E}_\mathrm{d}(\mathbf{R}, t),\]

where \({\bf d}=q{\bf r_{\rm dip}}\) is the dipole moment, \(\mathbf{\dot d}\) and \(\mathbf{\ddot d}\) are its first and second derivatives, \(\gamma_0\) is the free-space energy decay rate given by

(6)#\[\gamma_0 = \frac{q^2\omega_0^2}{6\pi\epsilon_0c^3m}.\]

This equation of motion corresponds to a Lorentzian atom model with transition frequency \(\omega_0\) and linewidth \(\gamma_0\) (where \(\gamma_0 \ll \omega_0\)), valid for non-relativistic velocities since relativistic mass effects are neglected.

PyCharge Implementation#

The driven Lorentz Oscillator equation of motion (Eq. (5)) is the core of dynamic dipole simulations in PyCharge. This ODE is solved numerically by the simulate() function.

You can create a Lorentz Oscillator source using the dipole_source() factory function. This function automatically constructs the necessary ODE function that simulate() will use to time-step the dipole’s state.

Quantum Correspondence#

It is useful to discuss how the decay dynamics of LOs relate to those of a quantum two-level system (TLS) in certain limits. In the weak-excitation (linear response) regime, we can connect the quantum mechanical equations of motion for a TLS to the classical LO equations by replacing \(q^2/m\) with \(q^2 f/m\), where \(f\) is the oscillator strength, defined as

(7)#\[f = \frac{2m \omega_0 d_0^2}{\hbar q^2},\]

where \(d_0 = |\mathbf{d}(t=0)|\). We then recover the standard expression for the spontaneous emission (SE) rate \(\gamma_{0,\mathrm{TLS}}\) from an excited TLS:

(8)#\[\gamma_{0,\mathrm{TLS}} = \frac{\omega_0^3 d_0^2}{3\pi\epsilon_0 \hbar c^3}.\]

An alternative argument relating dipole moment to radiative decay rate is to equate the total mean energy of the LO with the ground-state energy of a quantized harmonic oscillator:

\[\frac{m\omega_0^2 d_0^2}{q^2} = \frac{\hbar\omega_0}{2},\]

yielding \(q^2/m = 2\omega_0 d_0^2/\hbar\), as expected. Using Fermi’s golden rule and the interaction Hamiltonian \(H_\mathrm{int} = -{\bf d} \cdot \hat{\bf E}\), we obtain rate equations for the populations of an isolated TLS in vacuum:

(9)#\[\dot n_\mathrm{e}(t) = -\gamma_0 n_\mathrm{e}(t),\]
(10)#\[\dot n_\mathrm{g}(t) = \gamma_0 n_\mathrm{e}(t),\]

where \(n_g\) and \(n_e\) are the populations of the ground and excited states, respectively (\(n_g+n_e=1\)). Here, \(\gamma_0\) is identical to the Einstein A coefficient, so the energy decay rate equals the population decay rate. The analogy between LO dynamics and TLS populations holds only in the weak-excitation limit.

Energy of a Dipole#

The total energy \(\mathcal{E}\) of a dipole, which is the sum of its kinetic and potential energies, can be calculated as

(11)#\[\mathcal{E}(t) = \frac{m \omega_0^2}{2q^2} d^2(t) + \frac{m}{2q^2} \dot d^2(t),\]

where \(\dot d = |\mathbf{\dot d}|\).

In a simulation, the dipole moment \(\mathbf{d}(t)\) and its derivative \(\mathbf{\dot d}(t)\) are calculated from the state array returned by the simulate() function. The state array contains the position and velocity of each charge in the dipole at each time step.

Since the total dipole energy \(\mathcal{E}\) is proportional to the excited-state population, the normalized population can be determined from

(12)#\[n_e(t)=\frac{\mathcal{E}(t)}{\max (\mathcal{E})}.\]

References

    1. Novotny and B. Hecht, Principles of Nano-Optics, Ch. 8

    1. Milonni and J. Eberly, Lasers Physics, Ch. 3

  • D.J. Griffiths, Introduction to Electrodynamics, Ch. 11