Coupled Lorentz Oscillators

Note: See the appendix of the PyCharge paper for the quantum electrodynamics (QED) expressions for the radiative decay rates and coupling parameters of point dipoles.

It is well known that an atom's surrounding environment modifies its radiative properties. In the classical model, the modification of the SE rate is generated by the scattering of the atomic field (as the LO is driven by the electric field at its origin position), while in QED theory the SE rate is stimulated by vacuum field fluctuations or radiation reaction, which partly depends on the ordering of the quantum field operators. Regardless, in the weak coupling regime (where the atom-field coupling constant is much less than the photon decay rate inside the cavity), the interactions can be treated perturbatively such that QED and classical theory yield the same results for the modification of the SE rate. An exception is when the surrounding medium contains gain. The modification of radiative properties for two coupled LOs in close vicinity is given by invoking QED theory and using the dyadic Green's function for a dipole.

The classical analogs of the superradiant and subradiant states of two coupled TLSs (where the dipoles are quantized) occur when they are polarized along the same axis and begin either in phase (direction of the two dipole moments are equal) or out of phase (direction of the two dipole moments are reversed), respectively. PyCharge can calculate the frequency shift \(\delta_{12}\) and SE rate \(\gamma^{\pm}\) of two coupled LOs in either collective state by curve fitting the discretized kinetic energy (KE) values, which are calculated by PyCharge at each time step, to the expected harmonic equation

\[\begin{equation} \label{eq:KE} \mathrm{KE}= Ae^{-(\gamma^{\pm} t)}\sin\left((\omega_0\pm\delta_{12}) t + \phi\right)^2, \end{equation}\]

where \(A\) and \(\phi\) are constants necessary to accurately fit the function and are dependent on the initial conditions of the simulation. The curve fit should be performed using the kinetic energy values after a number of time steps have elapsed in the simulation to allow the scattered fields to propagate back to the LO's origin position. When the two coupled LOs are in the superradiant or subradiant states, the population of their excited state and their total energy \(\mathcal{E}\) are exponentially decaying functions with a decay rate of \(\gamma^{+}\) or \(\gamma^{-}\), respectively.

It is also useful to note that the total EM power radiated by an accelerating point charge in a vacuum (at non-relativistic speeds) can be calculated using the Larmor formula:

\[\begin{equation} P(t) = \frac{q^2 a^2(t)}{6 \pi \epsilon_0 c^3}. \end{equation}\]

The power radiated by a dipole can also be calculated using the above equation by replacing \(q^2 a^2\) with \(|\mathbf{\ddot d}|^2\). Assuming that the dipoles begin oscillating at \(t=0\) s, the radiated energy at time \(t'\) can be calculated by integrating the radiated power from \(t=0\) s to \(t=t'\) (which can be approximated with PyCharge using a discrete integration). As well, if there are two or more dipoles in a system that interact, then each dipole will `absorb' a certain amount of energy \(W_\mathrm{abs}\) radiated from the other dipoles. The total (constant) energy of a system that contains \(N\) dipoles is the sum of the energy gains and losses of all the dipoles, given by

\[\begin{equation} W_\mathrm{total} = \sum_{\mathrm{i}=1}^N \left( \mathcal{E}_i(t') - W_{\mathrm{abs},\,i}(t') + \int_0^{t'} P_i(t)\,dt \right), \end{equation}\]

where \(\mathcal{E}_i\) is the total energy (sum of the kinetic and potential energies) of the \(i\)th dipole in the system.