Direct space-time modeling of mechanically dressed dipole-dipole interactions with electromagnetically-coupled oscillating dipoles

The spontaneous emission (SE) process can be modelled from the viewpoint of vacuum fluctuations in quantum field theory or from radiation reaction [1], with the latter mapping on to classical and semiclassical field-matter theories, whereby the radiative dynamics can be captured through polarization dipoles and Lorentz oscillators (LOs) [2] coupled to classical fields.

In both classical and quantum field pictures, the SE rate is governed by the (projected) photonic local density of states (LDOS) [3, 4, 5, 6], which is related to the imaginary part of the system’s electric field Green function $\mathbf{G}$, through

for an emitter at position ${\bf r}_{0}$, with dipole moment ${\bf d}$ (assumed real) and resonance frequency $\omega_{0}$. The classical Green function of the medium is defined from

where $\epsilon({\bf r},\omega)$ is the dielectric constant of the medium, $c$ is the speed of light in vacuum, and ${\mathbf{I}}$ is the unit tensor (dyad), and we assume a nonmagnetic medium.

For example, for a dipole in free space, as depicted in Fig. 1(a), we have ${\rm Im}{G}_{ii}=\omega^{3}/(6\pi c^{3})$, and

is the free space SE decay rate, for an emitter with resonance frequency $\omega_{0}$ with $d=\lVert{\bf d}\rVert$. In a dielectric with $\epsilon_{\rm B}=n_{\rm B}^{2}$, this rate is simply modified by $n_{\rm B}$, the background refractive index of the medium. Moreover, for any inhomogeneous structure, the SE rate can be controlled by the LDOS (or projected LDOS, as the dipole direction also matters), which is modified by the electromagnetic environment surrounding the emitters. The pioneering work of Edward Purcell first pointed out that the SE rate of quantum emitters could be altered by placing them inside a resonant cavity [7], causing enhancement or inhibition of the SE rate. Subsequent research has extensively explored this process across various physical systems, such as emitters within photonic crystals [8, 9, 10, 11, 12, 13, 14].

An alternative approach for modifying the SE dynamics involves coupling to other atoms, or to a collection of atoms, typically in the near-field regime (i.e., much less than their emission wavelength). This results in collective excitations that decay either faster (superradiance) or slower (subradiance) than the individual emitters, which has been demonstrated in a number of systems [15, 16, 17, 18, 19]. These coupled dipole studies reveal that the modified SE rate is dependent on interatomic distances, which can significantly impact the amplitude and width of the emission spectrum [17, 20].

For two atoms, with their center-of-mass (COM) at positions ${\bf R}_{n}$ and with resonant energies $\omega_{n}$, with $n=1,2$, then we also have an incoherent photon transfer rate [21]¹¹1Here we also allow for complex dipole moments

and a coherent coupling rate

where a Markov approximation has been applied (e.g., the rates are assumed to be constant over the frequency range of interest, so all rates are instantaneous in time).

Note also that the above expressions can be derived classically or quantum mechanically, with a full classical-quantum correspondence (discussed in more detail in Appendix B). For example, in the quantum case, when the dipole bosonic operators $b_{n}$ and $b_{n}^{\dagger}$ are introduced for each atom, the interaction Hamiltonian for dipole-dipole coupling is $H_{dd}=\hbar g_{12}(b_{1}+b^{\dagger}_{1})(b_{2}+b^{\dagger}_{2})$, and within a rotating-wave approximation, $H_{dd}^{\rm RWA}=\hbar g_{12}(b_{1}^{\dagger}b_{2}+b^{\dagger}_{2}b_{1})$. With the appropriate initial condition, this can lead to the creation of superradiant states (Dicke states [22], specialized to two dipoles), which decay with the rate $\gamma^{+}=\gamma_{0}+\gamma_{12}$, or subradiant states, which decay at the rate $\gamma^{-}=\gamma_{0}-\gamma_{12}$. For simplicity, here we assume equal emitter strengths and resonance frequencies, so that $\gamma_{0}=\gamma_{11}=\gamma_{22}$, $\gamma_{12}=\gamma_{21}$, and $g_{12}=g_{21}\equiv g$ [see Fig. 1(b)].

Recently, the direct numerical modelling of dipole-dipole interactions, with no intrinsic approximations (such as the dipole approximation, rotating wave approximation or Markov approximation), was demonstrated by directly solving the self-consistent electrodynamics simulations of LOs and oscillating point charges [23]. This direct modelling approach can have enormous benefits, as it connects classical dynamical simulations to complex SE dynamics, often in ways that are intractable or difficult to account for in standard quantum theory (or standard classical theory). For example, one does not have to make a rotating wave approximation, and one can also go beyond a dipole approximation, two common approximations in quantized light-matter interaction theories.

A particularly intriguing way of modifying the SE dynamics is to move the COM of the emission dipole itself, e.g., modeled as a moving atom [24]. At relativistic velocities, this can yield very small corrections to $\gamma_{0}$ [25, 26]. A more practical and significant way to modify the SE rate is through mechanical perturbations of the photonic environment, resulting in an optomechanical coupling, e.g., where the SE can be altered by a vibrating cavity wall [27, 28, 29]. As is well known, if the atom is close to a mirror (e.g., a half-open cavity), the emitted photon can be reflected by the mirror, which (depending on the angle of emission, frequency, and precise distance from the mirror) can either create a constructive or destructive interference at the position of the atom [4, 30, 16, 31, 32].

Indeed, it is now well known that the SE rate and the effective level (frequency) shift of a quantum emitter, such as a two-level atom, can be modified by its interaction with an oscillating mirror. This interaction manifests in a positional-dependent SE rate and effective level shift, scaling with $\cos(2k_{0}R_{0})$ and $\sin(2k_{0}R_{0})$, respectively, where $k_{0}$ is the wave vector and $R_{0}$ is the atom-mirror separation [33, 34]. In Ref. [33], using a rotating wave approximation and one-dimensional quantized fields, the SE of a two-level system near an oscillating plate was investigated. The study focused on field modes propagating within a limited solid angle, $\theta,$ towards the mirror and reflected back to the atom, influenced by the mirror oscillation amplitude. In a Markov approximation (instantaneous coupling), when $\gamma\tau\ll 1$, $v_{\rm mirror}=R_{M}\omega_{M}\ll c$ and $\omega_{M}\tau\ll 1$ (with $\tau=2R_{0}/c$ the round-trip time, and $R_{M}$ and $\omega_{M}$ are the (mechanically) oscillating amplitude and frequency of the mirror), the modified SE rate was found to be

and similarly for the atom frequency, $\omega_{0}^{\prime}=\omega_{0}-\theta\frac{\gamma}{2}{J}_{0}(2k_{0}R_{M})\sin(2k_{0}R_{0})$, where $\theta$ is less than 1 and depends on the scattering geometry, and ${J}_{0}$ is the zeroth-order Bessel function of the first kind, which represents the mirror’s oscillation effect. Modulating the mirror with the dependence $R_{M}\sin(\omega_{M}t)$, yields the following amplitude for an initially excited atom:

where $|C(t)|^{2}$ is the time-dependent population. The excited atom population decays exponentially at a modified rate, $\gamma^{\prime}$, which includes an oscillating modulation, where the amplitude of this modulation decreases when the oscillation frequency of the mirror increases. In the resolved sideband regime, where $\gamma<\omega_{M}$, this modification is small.

In Ref. [34], the authors considered a more general model, where a full three-dimensional electromagnetic field is quantized with the time-dependent boundary conditions determined by the (adiabatically) oscillating mirror. A general orientation of the atomic dipole moment was considered, but the oscillating frequency is only a few GHz compared with the transition frequency of an atom ($\approx 1$ THz for that work), and the oscillating amplitude was $R_{M}/R_{0}=0.1$. Another work investigated the SE rate and emission spectrum of a TLS positioned within a dynamic photonic crystal, which showed how SE modifications relate to the time-dependent photonic density states [35]. These findings suggested that the dynamical control of the boundary conditions causes a further modification of the SE process [33, 35, 34]. This is a nontrivial regime as the SE process is now non-stationary because of the time-dependent boundary conditions.

Apart from the interest in fundamental optics, the investigation of vibrational atoms and molecules is also important in the context of Raman spectroscopy and Surface-Enhanced Raman Spectroscopy (SERS), which forms a model of molecular optomechanics [36, 37, 38]. Molecular optomechanics provides insights into the fundamental principles governing the interaction between light and matter, as well as the development of technologies or applications, such as ultrasensitive detectors that can measure tiny forces and displacements at atomic scale [38, 39, 40].

In this paper, we introduce a direct space-time model for coupling mechanical oscillations to radiatively coupled emitters, where we directly move (oscillate) the COM of the dipoles mechanically, e.g., instead of vibrating mirrors which cause phonons to couple with photons [33, 34]. The main numerical simulations are carried out using PyCharge, a recently developed Python package that can calculate the electromagnetic fields and (vector and scalar) potentials generated by moving point charges and can self-consistently simulate dipoles as LOs [23]. The electric dipoles (LOs) are formed by a coupled oscillating point charges with opposite signs, equal magnitude and a separated displacement [see Fig. 1(a)]. The charges oscillate around the COM of the electric dipole along the direction of polarization. The origin of the LOs can be set to stationary [see Fig. 1(b)] or move as a function of time [see Fig. 2(a)]. The LOs are naturally damped within the self-consistent solver, with radiative decay. The motion of the LOs is restricted to nonrelativistic velocities; however, PyCharge can model a wide range of effects in electromagnetism, including “Braking radiation” (Bremsstrahlung) [41].

There are several benefits to using PyCharge to model the dynamics of LOs. First, it has already been shown to be very accurate for obtaining the known rates of dipole-dipole coupling from QED theory (typical numbers are $<0.2\%$ relative error, even for very close dipole separations, deep sub-wavelength) [23]. Secondly, PyCharge can explore effects that cannot be solved analytically or via the usual simulation software for electromagnetism. For instance, the numerical simulation of PyCharge does not rely on a Markovian approximation, dipole approximation, nor any rotating wave approximation, which are often required in other approaches, such as quantum master equations.

Specifically, in this work we use PyCharge to explore the modified SE decay dynamics from a pair of atoms (treated as $s$-polarized dipoles) that are coupled through photon-coupling effects such as Förster coupling [42], when one of the atoms are subject to a mechanical oscillation on the COM. In general, an atom subject to a periodic oscillation will form Floquet states [43, 44, 45], which can then dress the usual dipole-dipole-induced superradiant and subradiant states to have multiple resonances, separated (in the perturbative limit) by $\pm l\omega_{M}$ ($l$ integer). Thus, we also develop a rigorous Floquet analysis of two-coupled atoms (harmonic oscillators) with mechanical perturbations, to help with the interpretation and prediction of our direct time-dependent solution using PyCharge. In PyCharge, the dipoles are modeled as classical LOs [4], by directly solving the equation of motion for the LOs including the scattered fields responsible for modified SE.

The rest of our article is organized as follows: In Sec. II we present the main theoretical background for the direct electromagnetic simulations in PyCharge. In Sec. III, we then describe an analytical model of two coupled quantum harmonic oscillators (HOs), with periodic modulation the dipole-dipole coupling, and show how this can be solved using Floquet theory. In Sec. IV, we present several example results showing the electric dipole dynamics of two coupled LOs with and without mechanical oscillations. We also show results for the emitted spectra, for different driving parameters, and each example regime is corroborated using the analytical Floquet theory. Finally, in Sec. V, we give a summary and conclusions. In addition, we present two Appendices. In Appendix A, we present the field solutions from oscillating point dipole in a homogeneous medium. While in Appendix B, we present a quantum and classical field theory for dipole-dipole coupling, which we use in the main text to derive a Floquet solution with periodic oscillations in the dipole-dipole coupling.

To numerically model the radiative decay through radiation reaction, we start with one initially excited dipole. Figure 2 shows a schematic of the two coupled LOs with mechanical oscillations governed by Eq. (8), as well as the potential energy levels. The LOs are formed by two opposite point charges, $\pm q_{n}$ and masses $m_{n\pm}$ for each dipole [having the effective (reduced) mass $m_{{\rm eff},n}=m_{n+}m_{n-}/(m_{n+}+m_{n-})$], with a separate displacement, which oscillates around the origin along the $y$-axis ($s$-polarized dipoles). We can then set an initial displacement for the dipoles, from $y_{n}(t=0)$. For example, we can set $y_{1}(t=0)=y_{0}$ and $y_{2}(t=0)\approx{0}$, where $y_{0}=1\,$nm. Thus, we consider Dipole 1 is initially excited, and Dipole 2 will become excited through interactions from the electric field generated by Dipole 1. As they oscillate, the LOs will radiate (dissipate) energy that causes the dipole moment to decrease [4], but potentially oscillate when coupled to nearby LOs.

Considering two spatially separated LOs with the same effective mass, $m_{{\rm eff},1}=m_{{\rm eff},2}=m_{{\rm eff}}$, the same charge, $q_{1}=q_{2}=q$, and the same transition frequency, $\omega_{1}=\omega_{2}\equiv\omega_{0}$, in free space, the coupled dipole equations of motion are

where ${\bf E}_{d}$ is the component of the total electric field in the direction of polarization generated by all sources in the system (including scattered fields), which is the sum of Coulomb and radiative electric field contributions, given in Eqs. (25) and (26) of Appendix. A, respectively. In the coupled dipole equation of motion, $\mathbf{R}_{1}(t)$ is the position vector of the COM of Dipole 1 which is in periodically mechanical motion, ${\bf R}_{2}$ (assumed fixed) is the position vector of the COM of Dipole 2. Since any relevant physics here depends on the relative distance, we define ${\bf R}(t)={\bf R}_{1}(t)-{\bf R}_{2}$. Thus, the mechanical motion is defined from

where we let ${\bf R}_{2}={\bf 0}$, ${\bf R}_{0}=R_{0}{\bf\hat{x}}$ and ${\bf R}_{M}=R_{M}{\bf\hat{x}}$, so that $R_{M}$ is the mechanical oscillating amplitude of Dipole 1, $R_{0}$ is the initial static distance between the COMs of the two dipoles along $x$-axis, and $\omega_{M}=2\pi/T$ is the mechanical oscillating frequency. Therefore, the relative time-dependent distance between the two LOs is simply $R(t)\equiv\lVert{\bf R}(t)\rVert=R_{0}+R_{M}\sin(\omega_{M}t)$ [see Fig. 2(a)].

For our study, we will assume that the average mechanical velocity of the dipole is $v_{M}\approx R_{M}\omega_{M}\ll c$. In addition, we constrain the movement of dipoles and charges to non-relativistic velocities, so that the dissipation of charges due to radiation reaction damping is negligible. Thus, we set a reasonable maximum velocity of charge to be $c/100$. Thus, the mechanical oscillation parameters $R_{M}$ and $\omega_{M}$ play the significant role in our study, and we primarily focus on their respective dimensionless ratios, denoted as ${R_{M}}/{R_{0}}$ which is constrained to values less than $0.8$ (this limitation is imposed to ensure that the dipoles remain spatially separate and do not come in direct contact with each other) and ${\omega_{M}}/{g}$.

In the SE rate regime, the energy decay rate of an emitted photon is equivalent to the population decay rate, where the free-space decay rate, $\gamma_{0}$, can also be modelled as an Einstein A coefficient [46]. However, it is important to emphasize that the relationship between the dynamics of LOs and the population of two-level systems’ states can only be connected in a weak excitation regime (since the LO model does not account for a two-level system saturation effects).

To connect to the LO equivalent of population decay, one can also calculate the self-consistent dipole moment of the LOs at each time step by solving Eq. (8). The total energy $\mathcal{E}$ of an oscillating dipole, is then given by

which is the sum of the potential energy and kinetic energy of the dipole, respectively, and $d_{n}(t)$ is the time-dependent dipole moment of the $n$th dipole. We assume that the total energy $\mathcal{E}_{n}$ is proportional to the population of the excited state (mapping on to a two-level system model), due to the equivalence of the energy decay rate and the population decay rate. Subsequently, we can write the scaled population of the excited states as the normalized total energy [23],

for the $n$th dipole.

As described in Appendix. B, the system Hamiltonian for two coupled quantum HOs, when considering the coupling between two electric dipoles in free space, with equal oscillator frequencies $\omega_{0}$, takes on the form

where $b_{n}$ ($b_{n}^{\dagger}$) is the bosonic annihilation (creation) operator of the $n$th dipole and we assume the dipoles are coupled through a near-field electromagnetic interaction, with $g=q_{1}q_{2}/(8\pi\hbar\epsilon_{0}m_{\rm eff}\omega_{0}R^{3})$, where $R=\lVert{\bf R}_{1}-{\bf R}_{2}\rVert.$ More generally, this interaction strength can be obtained from Eq. (5), which depends on the real part of the Green function, ${\bf G}({\bf R}_{1},{\bf R}_{2})$. A correspondence between the classical dipole model can easily be made by connecting $q_{n}^{2}/m_{{\rm eff},n}=2\omega_{n}d_{n}^{2}/\hbar$ (see Appendix. B, and also Ref. [23]), and in terms a quantum dipole moment, one can also show that $g=g_{12}$, introduced earlier (see also Refs. [47, 48]).

Below we will use the Hamiltonian for the quantum HO problem, but we stress that exactly the same resonances occur for a classical Hamiltonian system, which is easy to recover by simply replacing the harmonic oscillator operators by (classical) $c$-numbers. This classical correspondence is possible since the system obeys a linear response, where the resonances yield the polariton frequencies. Interestingly, note that the specific form of the system Hamiltonian is similar to a quantum Hopfield-like model without the diamagnetic term [49, 50], which is needed for describing the coupling between a quantum dipole and a single quantized cavity mode, otherwise a phase transition (often ill-defined) can occur when $\eta=g/\omega_{0}>0.5$. In the present case, such a limit cannot be reached, since the dipole approximation is only valid when $R>4a$, where $a$ is the size (length) of the dipole (for example, the dipole moment $d=ea$).

In order to understand the dynamics of the system and calculate the normal modes or resonances (in the absence of dissipation), we write the Heisenberg equation of motion as $\mathrm{d}\bm{B}/\mathrm{d}t=-\mathrm{i}\mathbf{H}\,\bm{B}$, with $\bm{B}=[b_{1},b_{2},b_{1}^{\dagger},b_{2}^{\dagger}]^{T}$, where

and the normal modes of the hybrid system are $\omega_{\pm}=\sqrt{\omega_{0}^{2}\pm 2g\omega_{0}}$, which are net positive since $g<\omega_{0}/2$ for dipole-dipole interactions within the dipole approximation.

We next present an analysis of the two coupled HOs with relative COM movement where the system is Floquet engineered through the time-dependent dipole-dipole coupling rate $g\to g(t)$. This time-dependent character of the system is through the time-dependent factor to the interaction Hamiltonian $H_{\rm dd}={d}_{1}{d}_{2}/(4\pi\epsilon_{0}R^{3})$ via $R=R_{0}\to R(t)=R_{0}+R_{M}\sin\omega_{M}t$, which gives $g\to g(t)$, where

Due to the periodicity of the time-dependent coupling, with period $T$, the Hamiltonian is also periodic: ${H}(t)={H}(t+T)$, and so is the dynamical matrix of the equation of motion: $\mathbf{H}(t)=\mathbf{H}(t+T)$.

We can thus expand the time-dependent Hamiltonian as a Fourier series $\mathbf{H}(t)=\sum_{m\in\mathbb{Z}}\mathbf{H}_{m}\,\mathrm{e}^{\mathrm{i}m\omega_{M}t}$, with $\mathbf{H}_{m}=({1}/{T})\int_{T}\mathrm{d}t\,\mathrm{e}^{-\mathrm{i}m\omega_{M}t}\,\mathbf{H}(t)$, yielding

where

Notably, we see that the static coupling rate, $g_{0}$, is renormalized by the average of the time-dependent term [51]:

and thus the mechanical drive renormalizes the dc Hamiltonian through an effective increase in $g$, which is the undriven coupling term. Later, we will see such an effect in our simulations when $R_{M}$ is sufficiently large.

Note that $\mathbf{H}_{m}$ separates into a time-independent part for $m=0$ and $m\neq 0$ (for finite $R_{0}$, including a shift due to $R_{M}\neq 0$), and a time-dependent interaction. Thus, while $R_{M}$ is related to the time-dependent interaction, there is a static contribution from that term as well. Solving the time-dependent equation of motion (which is similar to the time-dependent Schrödinger equation with the Floquet theory), yields $\bm{B}_{\alpha}(t)=\mathrm{e}^{-\mathrm{i}\varepsilon_{\alpha}t}\bm{F}_{\alpha}(t)$, where $\varepsilon_{\alpha}$ is the Floquet quasienergy [52], and the Floquet mode $\bm{F}_{\alpha}(t)$ is $T$-periodic [53, 54]. The Floquet solutions, $\{\bm{B}_{\alpha}(t)\}$, form a complete basis for any value of $t$, thus $\bm{B}(t)=\sum_{\alpha}c_{\alpha}\bm{B}_{\alpha}(t)$, where $c_{\alpha}=\bm{F}_{\alpha}^{T}\bm{B}_{\alpha}(0)$, with $\bm{F}_{\alpha}\equiv\bm{F}_{\alpha}(0)$. Resonances occur at differences between Floquet quasienergies [55].

To compute the Floquet modes [51], we use a Fourier series expansion of $\bm{F}_{\alpha}(t)=\sum_{l\in\mathbb{Z}}\mathrm{e}^{\mathrm{i}l\omega_{\rm M}t}\bm{F}_{\alpha l}$, where the Fourier coefficient vectors $\bm{F}_{\alpha l}$ are Floquet sidebands. Although the dynamical matrix is $4\times 4$ so that one might expect nominally $4$ eigenenergies and thus 4 quasienergies confined within a $\omega_{M}$ energy range, we only see 2 lines corresponding to the two modes (the determinant has only two real roots, the polariton resonances) as the Floquet theory is applied to the equation of motion (not the Hamiltonian) representing the effective transitions–from the dressed polariton states to the dressed ground state, within the system [51]. Indeed, the coupling between these quasienergy lines is drive-dependent and for small enough $\eta$, $R_{M}$ and $\omega_{M}$, the drive might not have enough push to couple all the manifolds of the energy states [see the right-most part of Fig. 2(b)], and thus, only a limited number of them come into a Brillouin zone to make a coupling and transition possible. This will also reveal the dynamics of the two dressed polariton modes, that is equivalent in classical and quantum pictures [50].

Although we have made several key approximation in the theory above, we will see later that the resonances are in very good agreement with the full electrodynamic simulations. To be specific, PyCharge employs a full wave approach without any approximation based on nonperturbative, linear electromagnetic theory (derived from Maxwell’s equations), so it can capture the nondipolar, and non-Markovian aspects of the coupling, as well as the damping. However, the Floquet analysis above properly accounts for the main effect of the modulated dipole-dipole interaction including the intrinsic asymmetry stepping from the $1/R^{3}(t)$ term, while giving a solid theoretical background to the main effects that we expect to see when we modulate the dipole-dipole coupling with a periodic drive. Moreover, the analysis can also be used to look at quantum effects for coupled two-level systems, which yields different spectral signatures than the coupled harmonic oscillator problem (primarily due to saturation effects and biexcitonic coupling). This will be the subject of future work.

In this section, we show the results of the PyCharge simulations in time and frequency domains, and show its agreement with the Floquet theory. We let the first dipole be initially excited (internal charges are within their maximum relative distance) and the second one initially in the ground state (internal charges are within their minimum relative distance). In Figs. 3, 4 and 5, we first investigate the time evolution of the scaled population of the excited states of the two coupled LOs ($s$-dipoles) in their panels (a), then show the Fourier transformation of the time-dependent dipole moments that shows the spectral peaks of the transitions in their panels (b), and finally show their corresponding Floquet quasienergy diagrams in panels (c-d) to explain that the transitions occur in the differences of the quasienergies. We use a fixed value of the (static) coupling between the dipoles of $g/\omega_{0}=0.00136=g$, which is large enough to easily resolve the dipole-dipole splitting and coherent oscillations, since $g\gg\gamma_{0}$. In addition, we also choose various values for the mechanical drive’s frequency and dynamical coupling. Specifically, we first use $R_{M}=0.1R_{0}$ and a rather low drive frequency of $\omega_{M}=g$, which is already large enough to show clear new resonances, as shown in Fig. 3; in addition, we investigate $R_{M}=0.1R_{0}$ and a higher drive frequency of $\omega_{M}=5g$ in Fig. 4, and we increase the dynamical coupling by $R_{M}=0.35R_{0}$ and a rather high drive frequency of $\omega_{M}=5g$ in Fig. 5. Below we discuss each of these three regimes.

Figure 3(a) shows the PyCharge simulation result of the (scaled) population of the excited states of the two coupled dipoles with (Dipole 1^′ and $2^{\prime}$, solid lines) and without (Dipole 1 and 2, dashed lines) mechanical oscillations, plotted as a function of time with the following mechanical parameters: $R_{M}=0.1R_{0}$ and $\omega_{M}\approx g$ ($T=0.079\gamma_{0}^{-1}$), where $T$ refers a cycle of mechanical oscillation period. The time is scaled by the free-space decay rate ($\gamma_{0}$ = 21.48 GHz). Dipole 1 and 1^′ [$y_{1}(0)=1\,$nm] are initially excited, while Dipole 2 and 2^′ are not [$y_{2}(0)\approx 0$], corresponding to the schematic diagram in Fig. 2(a). All of the LOs have a natural angular frequency $\omega_{0}$ of 200 THz. To ensure very good accuracy, the numerical simulation runs over 10 million time iterations (125,663 periods of $\omega_{0}$) with a time step $\Delta t=10^{-17}\,$s [$(\omega_{0}\Delta t)^{-1}\approx 80$]. The dipole-dipole initial separation is $R_{0}=50\,$nm ($y_{0}/R_{0}=0.02$) along the $x$-axis, with charge magnitude $q=10e$, effective mass $m_{\rm eff}=\hbar/2\omega_{0}y_{0}^{2}$ and static coupling strength $g\equiv g_{12}=g_{21}=79.7\gamma_{0}$ ($\tau_{dd}=\pi g^{-1}=0.039\gamma_{0}^{-1}$, where $\tau_{dd}$ represents the period of a cycle in which a photon is emitted and subsequently re-absorbed by the same dipole, or simply the round trip time in a classical picture). The time domain simulations show an oscillatory decay (caused by dipole-dipole interactions) of the populations to a steady state, where, with the onset of the mechanical oscillation, the phase of the oscillatory decay changes dynamically; thus, the time interval between the same peaks of the moving and stationary dipoles is sometimes shorter and sometimes longer. However, this oscillatory decay behavior for the mechanically oscillating dipoles is always faster than the stationary dipoles; this indicates the emergence of higher frequency peaks in the spectral analysis [cf. spectral peaks in panel (b) and also panel (c)].

Figure 3(b) displays the spectral peaks of the dipoles in the frequency domain (shown in a logarithmic scale) for the corresponding cases in Figs 3(a), centered at the dipoles’ transition frequency $\omega_{0}$ and normalized by $\omega_{M}$. The spectra are obtained by a fast Fourier transform (FFT) of the dipole moment $d_{n}(t)$ with a Hamming window function applied to enhance peak clarity and smoothness. We see that all of the resonance peaks from the field spectrum plots are separated by an integer number of the driving frequency $\pm\omega_{M}$ ($n$ integer) when $\omega_{M}\leq g$, as shown in Fig. 3(c). The blue and red dot-dashed vertical lines in panel (b) are the same as the blue and red thin solid lines representing the Floquet quasienergies, calculated from Eq. (15), in panel (c), for the desired values of the parameter space pointed out by the vertical transparent gray line, and with a magnified view in panel (d). This spectral behavior contrasts significantly with that of two coupled LOs with a fixed COM, which produce only two frequency peaks (the well known superradiant and subradiant states, separated by $2g$, which also connect to Dicke states [22]). As anticipated, mechanical oscillation can modify these states to a manifold of dressed (Floquet) states, shown in Fig. 2, where energy levels are separated by integer multiples of $\omega_{M}$.

Figure 3(c) shows the Floquet quasienergies near $\omega_{0}$, scaled with $\omega_{M}$, where $\omega_{M}/\omega_{0}=\eta=0.00136$, for two coupled HOs as a function of normalized oscillating amplitude ($R_{M}/R_{0}$) for the same regime as that of panels (a) and (b). Here, we include a truncated number of Brillouin zones (separated by the grey dotted lines). Construction of the Brillouin zones, which is dependent on the strength of the external drive, assists the coupling of the energy states by the external quanta; thus, with a limited (not enough) number of Brillouin zones, due to the small choices of drive frequency and amplitude, the complete intermixing of all the system’s eigenenergy states might not be possible within a Brillouin zone. Moreover, within the small range of the drive amplitude, we do not yet observe the creation of effective anticrossings in the dynamical picture between the quasienergy levels, which would only occur for sufficiently larger values of $R_{M}/R_{0}$.

In coupling scenarios where $\omega_{M}\leq g$, the spectral lines are initially separated by the integers of the driving frequency, as demonstrated in Fig. 3(c), where we normally expect a single peak structure within the range of one $\omega_{M}$. We also see that at $\omega_{M}\sim g$ (and, in general, increasing the coupling strength via $g$ and $R_{M}$) we are well at the verge of entering into the regime of splitting the peak, as seen in panel (b), where also the initial splitting of the Floquet quasienergy lines almost begin to occur near the vertical grey line in panel (c); which is clearer in the magnified region in panel (d). However, a split in the spectral lines with an internal separation of $2\eta$ is observed when $\omega_{M}>g$ (and, more effective for $\omega_{M}\gg g$), because the splitting of the quasienergy lines is now more effective and also they may undergo through an ac Stark shift. Additionally, the center of splits is still separated by an integer multiple of $\omega_{M}$ because still we encounter transition only between the self repetition manifolds (going back and forth in different Brillouin zones from a state) of the polariton states only.

Figure 4 is similar to Fig. 3, but with a larger driving frequency, $\omega_{M}\approx 5g$ ($T=0.016\gamma_{0}^{-1}$). Comparing the results with (primed dipoles) and without mechanical oscillation (unprimed dipoles), i.e., comparing panel (a)’s in Figs. 3 and 4, the driving frequency clearly alters the phase of the population oscillations in a different way due to time-dependent fluctuation in the dipole-dipole coupling strength. As shown in Fig. 4(a), we observe an overall faster exchange of energy between Dipole 1^′ and 2^′ (mechanical) compared with Dipole 1 and 2 (non-mechanical), when the driving frequency is greater than the dipole-dipole frequency shift, $\omega_{M}=5g>g$, indicating an increase of the frequency shift. In other words, for $\omega_{M}=g$ in Fig. 3(a), we see that the populations of mechanical dipoles (1’ and 2’) oscillate sometimes faster, sometimes slower than the non-mechanical cases (1 and 2). However, once we have a larger driving frequency, the population exchange of mechanically moving dipoles are always faster than the non-mechanical case. This can imply that for $\omega_{M}\gg g$, one can expect an overall faster exchange of energy between Dipole 1^′ and 2^′ (mechanical).

Furthermore, with $\omega_{M}=5g>g$, we observe that the resonance peaks tend to exhibit a splitting into superradiant and subradiant states with internal separation around 2$g$, as illustrated in Fig. 4(b). The inset of panel (b) provides an enlarged view of the peaks associated with the dressed polariton states. Additionally, the center of the splitting peaks is still separated by an integer of driving frequency. The vertical blue and red dot-dashed lines show these phenomena are aligned with our Floquet analysis, showing the excellent agreement between Floquet analysis and PyCharge simulations.

We next increase the oscillating amplitude that can bring the dipoles much closer to each other to explore the role of the oscillating strength (of the mechanical drive), again with $\omega_{M}=5g$, as shown in Fig. 5, where we now use a larger driving amplitude, $R_{M}=0.35R_{0}$. The change in the phase of population oscillations has not been significant, but the frequency of the population oscillations increases and decreases dramatically, as the effective $\tau_{dd}$ is significantly reduced when the dipoles are brought closer together, and significantly increases when the dipoles are farther apart. This causes the resonances to be modified dynamically, as shown in Fig. 5(b-d).

We also find that the internal separation between resonance peaks got further split when $R_{M}/R_{0}$ is large enough. Our study indicates that the oscillating amplitude $R_{M}$ is able to influence the positioning of resonant frequencies, causing a splitting of the original peaks, as well as a change in the separation between the split peaks (when the amplitude is sufficiently large). The ensuing dynamic control is also corroborated by Floquet analysis, see panels (c) and (d) of Fig. 5, where we see the renormalization of the dc eigenenergies (dashed lines) and the ac Stark shift of quasienergies (thin solid lines) are becoming largely impactful and are well distinguishable from the dc eigenenergies (thick solid lines).

In Figs. 3, 4 and 5, Floquet quasienergies and states are obtained by analyzing the system’s Floquet dynamical matrix corresponding to the equation of motion, which was constructed in an extended Hilbert space incorporating time periodicity, similar to the approach in Ref. [51], and their correspondence to spectral features of the dipole moment $|{d}_{n}(\omega)|$. Hence, this approach confirms the presence of well-defined Floquet states, with spectral peaks appearing at integer multiples of the driving frequency. Additionally, the theoretical Floquet model was validated by the numerical simulations, providing insights into the effects of mechanical motion and dipole-dipole interactions on the system’s energy structure.

Additionally, a closer examination of the results in the Floquet analysis reveals that, when the parameter $R_{M}/R_{0}$ reaches sufficiently high values, the quasienergy states exhibit crossing/anticrossing behavior. The region of crossing and avoided crossing occurs more rapidly with increasing oscillating amplitude, resulting in a dispersion and merging of the resonances. This phenomenon is better illustrated in Fig. 5(c). The region of crossing/anticrossing shifts toward larger values of normalized driving amplitude as the frequency of the drive increases. Thus, the amplitude modulation may decreases as driving frequency increases, similar to the findings in [33], and the amplitude of this modulation decreases when the oscillation frequency of the mirror increases.

Finally, we remark that we have also investigated the influence of the mechanical oscillation phase, along with the consequences of altering the signs of sine and cosine functions. We found that these effects are minimal when the mechanical frequency is sufficiently high. However, they can become more pronounced at lower driving frequencies, which was not the domain of interest for this study.

We have presented a direct numerical simulation of moving point charge potentials, using PyCharge, and studied how the (scaled) population dynamics and dipole spectra of two two-dipole coupled LOs can be significantly influenced by periodic mechanical oscillations, in a non-trivial way. We also investigated how the spectral characteristics and resonance peaks vary with the mechanical perturbation parameters, $R_{M}/R_{0}$ and $\omega_{M}/g$, and showed how the resonance peaks positioning, spectral splitting of the peaks, as well as the separation between the split peaks change with an increased oscillating amplitude $R_{M}$.

These results were fully corroborated with an intuitive and rigorous Floquet theory, where we added a periodic oscillation to the usual dipole-dipole coupling term (which can be derived from classical or quantum Maxwell equations), which couples two HOs through a time-dependent coherent coupling term. Both the direct numerical and Floquet analysis demonstrated a reciprocal relationship between $R_{M}/R_{0}$ and $g$, with resonances showing a separation by an integer of driving frequency, namely $\pm l\omega_{M}$, and these resonances further split into superradiant and subradiant states when drive frequency becomes larger than the static dipole-dipole coupling. Significantly, in scenarios with sufficiently large values of both $R_{M}$ and $\omega_{M}$, there was a reduction in the radiative decay rate of the system, implying that mechanical effects can effectively enhance the coupling strength and extend the lifetime of the coupling system.

Related effects have been predicted in studies of an atom in front of an oscillating mirror (phonon-photon coupling) [33, 34], and time-modulated photonic crystals (time-modulated band-gap) [35], where the sidebands in the emitted spectrum are separated from the atomic transition frequency by the modulation frequency of the plate or photonic lattice, $\omega_{M}$; however, the lateral peaks here were strongly asymmetric due to the different density of states at the edges of the photonic band gap [35]. In contrast, for our case, which deals with finite-size dipoles in free space, the lateral peaks are symmetric because the photonic density of states (DOS) is essentially the same at frequencies of the peaks.

The advantage of PyCharge is that we do not need to make any of the usual approximations (dipole approximation and rotating wave approximation) and although everything is classical in nature, we are still able to show clear analogies with quantum optical systems (quantum emitters and quantized fields) [33, 34]. Additionally, we have demonstrated with a stronger oscillating strength or faster driving frequency, the resonance peaks can be further modified. We have gone beyond usual perturbative optomechanical studies, where they have linearized Taylor series expansion [56, 57, 58].

Overall, our findings provide a deeper understanding of the interplay between light and matter with mechanical motion, using coupled atoms described by LOs and the self-consistent fields. This work lays the foundation for exploring more complex optomechanical systems, potentially leading to novel applications in fundamental nano-optics, quantum information processing, and sensing and communication protocols (e.g., using coupled atoms as nano antennas). Moreover, it showcases the opportunities of using numerical solutions of retarded potentials and moving point charges and dipoles (using and extending the open-source PyCharge software), in contrast to the usual numerical Maxwell solvers with the physical fields (such as through FDTD and finite-element solvers).

Finally, we remark these classical studies can also be used to help develop more advanced quantum models, since a classical-quantum correspondence is only expected in some limit, e.g., linear response with no saturation effects (HO models). Thus, future work, for example, could look at the dynamical coupling between dipoles treated as two-level systems, where the eigenstates also contain a two quanta excitation, and saturation and nonlinear effects can become important.