Stochastic thermodynamics of inertial-like Stuart–Landau dimer

Our starting point is the system of four coupled underdamped Langevin equations that govern the dynamics of each pair of coordinates (x_i , y_i ) of the oscillators:

$$\begin{equation}m{\ddot{x}}_{1}=(R-{x}_{1}^{2}-{y}_{1}^{2}){x}_{1}-\gamma {\omega }_{1}{y}_{1}-k\enspace ({x}_{1}-{x}_{2})-\gamma {\dot{x}}_{1}+{\xi }_{1}^{R},\end{equation} \tag{ 1 }$$

$$\begin{equation}m{\ddot{y}}_{1}=(R-{x}_{1}^{2}-{y}_{1}^{2}){y}_{1}+\gamma {\omega }_{1}{x}_{1}-k({y}_{1}-{y}_{2})-\gamma {\dot{y}}_{1}+{\xi }_{1}^{I},\end{equation} \tag{ 2 }$$

$$\begin{equation}m{\ddot{x}}_{2}=(R-{x}_{2}^{2}-{y}_{2}^{2}){x}_{2}-\gamma {\omega }_{2}{y}_{2}-k({x}_{2}-{x}_{1})-\gamma {\dot{x}}_{2}+{\xi }_{2}^{R},\end{equation} \tag{ 3 }$$

$$\begin{equation}m{\ddot{y}}_{2}=(R-{x}_{2}^{2}-{y}_{2}^{2}){y}_{2}+\gamma {\omega }_{2}{x}_{2}-k({y}_{2}-{y}_{1})-\gamma {\dot{y}}_{2}+{\xi }_{2}^{I}.\end{equation} \tag{ 4 }$$

Both oscillators have the same mass m. The parameter R represents the strength of the external potential, and acts as a control parameter such that in the overdamped-uncoupled-noiseless limit (γ ≫ m, k → 0, and ${\xi }_{i}^{R(I)}(t)=0$) the oscillators settle in their own limit cycle with radius $\sqrt{R}$ for R > 0 and to a fixed point if R < 0. Each two-dimensional oscillator is governed by a frequency ωi (scaled by the damping coefficient) and the oscillators are coupled via a spring with strength k. The parameter γ represents the phenomenological Stokesian damping coefficient. The random noise ${\xi }_{i}^{R,I}$, i = 1, 2 are Gaussian with zero mean and $\langle {\xi }_{i}^{\alpha }(t)\;\;{\ast}\;\;{\xi }_{j}^{\beta }({t}^{\prime })\rangle =2{\gamma }_{j}{k}_{\text{B}}{T}_{j}\delta (t-{t}^{\prime }){\delta }_{\alpha ,\beta }{\delta }_{i,j}$, where ⟨⋅⟩ represents the ensemble average over the noises. The strength of the noise is chosen such that the fluctuation dissipation relation holds. Throughout this work we will set R = 1 as the unit scale in which other parameters are measured, ω₁ = 1 and ω₂ = ω₁ + Δω, unless mentioned otherwise. The temperatures T₁ and T₂ could be different, but we will focus on the case where they are the same, so that the nonequilibrium driving comes only from the deterministic part of the dynamics.

A potentially useful and equivalent way to rewrite these equations is to use the complex notation z₁ = x₁ + iy₁, z₂ = x₂ + iy₂, which yields the following compact expressions:

$$\begin{align}m{\ddot{z}}_{1}& =\left(R+\text{i}\gamma {\omega }_{1}-\vert {z}_{1}{\vert }^{2}\right){z}_{1}-\gamma {\dot{z}}_{1}+k({z}_{2}-{z}_{1})+{\xi }_{1},\\ m{\ddot{z}}_{2}& =\left(R+\text{i}\gamma {\omega }_{2}-\vert {z}_{2}{\vert }^{2}\right){z}_{2}-\gamma {\dot{z}}_{2}+k({z}_{1}-{z}_{2})+{\xi }_{2}.\end{align} \tag{ 5 }$$

Above the complex noise ${\xi }_{i}={\xi }_{i}^{R}+\text{i}{\xi }_{i}^{I}$. The above set of equations with complex variable z_i converges to the standard Stuart–Landau dimer in the overdamped limit γ ≫ m and the underdamped version will hereon be referred to as the inertial-like Stuart–Landau dimer (ISLD). Note that we consider the coordinates (x_i , y_i ) of the oscillator to be independent of each other and to describe position coordinates, due to which the presence of the second-order term refers to an inertia-like term. The independence between x and y is not necessary [43], but it is true in several systems of interest, e.g. neuron firing experiments wherein for a single SLO there are two independent variables namely the amplitude of the pulse and the duration between two firings [25] or nano-electro mechanical systems wherein if one goes beyond the slow variation of the amplitude approximation then a second order derivative for the complex variable z appears naturally [44, 45]. In the case of x and y being canonically conjugate (in the sense of Hamiltonian dynamics) the standard Stuart–Landau model itself would include an inertial term and the complex variable equations can be dealt with the use of Wirtinger derivatives. Throughout this work, we do not consider this case and refer the interested readers to references [46, 47]. Since even the standard Stuart–Landau equation could include effects of inertia, we use the term inertia-like to avoid conceptual pitfalls.

To the best of our knowledge, the noisy ISLD has never been explored before. In the remainder of this section we will fully explore the phase space and bifurcation diagram of this model in the zero temperature limit to illustrate the existence of bi-stability that arises only in the underdamped regime.

The physical steady-state phases of the ISLD are classified in figures 1(a) and (b) for the noiseless limit with small and large frequency difference Δω, respectively. We first focus on the noiseless case that exhibits a complex phase diagram and then explain the effect of noise. For simplicity, we have set γ = 1.

Zoom In Zoom Out Reset image size

Noiseless without an inertia-like term.- In the case without an inertia-like term (m = 0), there are three distinct phases, an ILC, amplitude death (AD), and an in-phase coherent limit cycle (CLCi). For γΔω < 1 the AD regime vanishes and we have an exceptional point that demarcate the ILC to CLCi transition. For γΔω > 1 the AD region emerges such that its area in the k–γΔω plane increases [48, 49]. In this overdamped case, the boundaries of the AD regime can be analytically predicted using linear stability analysis and for γΔω < 1 the critical coupling k = k_c = γΔω/2 [49] (exceptional point) segregates the ILC and CLCi regions.

The physical picture in this regime is quite intuitive, for small k each oscillator vibrates with its own frequency and hence there is no notion of synchronization giving an ILC due to the absence of noise (T_i = 0), figure 1(c). Note here that the second oscillator (red solid line) has a very small oscillation in its amplitude not visible due to the scale implying that it is indeed an ILC oscillator. As k increases, the frequency mismatch between the oscillators decreases and the oscillators motion destructively interfere to give rise to a fixed point solution known as AD, figure 1(f). For large values of k, since there is a strong coupling between the oscillators a self-organizing mechanism kicks in and the oscillators synchronize with their frequencies and amplitudes becoming equal with a finite phase difference as shown in figure 1(d). This is known as the CLCi regime since the phase differences are smaller than π/2. As the k increases further, the phase difference decreases and in the limit k → ∞ the oscillators are perfectly in phase.

Noiseless with an inertia-like term: small frequency difference γΔω < 1.- Adding an inertia-like term to the system leads to rich limit cycle dynamics. If m is sufficiently large with γΔω small, incoherent (ILC) and coherent (CLCi) limit cycles can coexist as seen in figure 1(a). Thus, the ISLD exhibits bistability, similar to that of Kuramoto model with inertia [50]. Such bistable behavior is also known to exist in dimers exhibiting explosive synchronization for large complex networks [51, 52]. The boundaries of various phases are no longer analytically tractable unlike the m = 0 case. The boundaries that we can obtain from linear stability is one that separates the AD region (orange line in figure 1(b)) and the red line in figures 1(a) and (b) that corresponds to an exceptional 'line' (see appendix A for more details).

Figures 2(a) and (b) show the bifurcation diagrams of amplitudes, |z_i |, of two oscillators when m = 2.0 and γΔω = 0.5 as a function of k. In order to obtain the bifurcation diagram, we initialize the system in a state given by the asymptotic limit of the previous coupling k − δk plus a small perturbation (forward direction, black lines in figure 2). In the ILC regime ($k< {k}_{\text{c}}^{\text{f}}$), the amplitude oscillates as seen by the spread of |z_i | (black region hidden under red in figures 2(a) and (b). At the critical value of the coupling ${k}_{\text{c}}^{\text{f}}$ a small fluctuation of initial conditions triggers a discontinuous jump and the oscillators synchronize in phase (black line, hidden under red, in figures 2(a) and (b)). In the reverse direction [53], (red lines in figure 2) for $k > {k}_{\text{c}}^{\text{b}}$ the oscillators remain synchronized in phase (CLCi) and we observe a single red line in figures 2(a) and (b). At the backward critical value of the coupling ${k}_{\text{c}}^{\text{b}}$ the system transits from CLCi to an ILC which is maintained for all $k< {k}_{\text{c}}^{\text{b}}$. Thus, the bifurcation diagrams not only allow us to obtain all the phase boundaries numerically, e.g. the boundary (black solid line) separating the bistable I and CLCi phase in figure 1(a), but also clearly shows the hysteric synchrony of incoherent and coherent states in the region between ${k}_{\text{c}}^{\text{b}}$ and ${k}_{\text{c}}^{\text{f}}$ [29, 30] that crucially depends on the choice of initial conditions indicating bistability.

Zoom In Zoom Out Reset image size

Noiseless with an inertia-like term: large frequency difference γΔω > 1.- When the difference of frequencies of the two oscillators is large the dynamics becomes highly complex in the presence of an inertia-like term as compared to the small γΔω regime. In this parameter regime, there exists an out of phase coherent limit cycle (CLCo) in which both oscillators have the same frequency but their amplitudes are different and they differ by a phase larger than π/2 (see figure 1(e)). At small mass m, unlike the low frequency difference case, a new boundary emerges (solid black line in figure 1(b)) that separates the ILC and CLCo regime. For large m, the boundary at the critical coupling k = k_c = γΔω/2 re-emerges and separates either the CLCo and bistable or ILC and bistable regime. The large frequency difference regime also produces a new bistable regime in which two different types of coherent states can coexist, CLCi and CLCo (see blue region in figure 1(b)). This type-II bistable region is distinctly different that the type-I, encountered in the small frequency difference regime that permits the coexistence of ILC and CLCi. At large m, both type-I and type-II bistable regions exist at large frequency differences, giving a rich complex phase diagram.

Figures 2(c) and (d) show the bifurcation diagrams of amplitudes, |z_i |, of two oscillators when m = 0.5 and γΔω = 4.0 (large frequency difference) as a function of k. As k increases (black line hidden under red), the amplitude |z_i | changes from the ILC to CLCo at a critical value ${k}_{\text{c}}^{0}$ and hence we see the spread reduce to a single line around k ≈ 1.3. As k increases further, the oscillators are synchronized (CLCo) up to a $k={k}_{\text{c}}^{\text{f}}$. At this point, a small fluctuation of initial conditions triggers a discontinuous jump to the CLCi regime. In the reverse direction when k is reduced (red lines in figures 2(c) and (d)), the oscillators maintain an in-phase synchronous behaviour (CLCi) up to a different $k={k}_{\text{c}}^{\text{b}}$ at which the system transits to being out-of-phase (CLCo). At ${k}_{\text{c}}^{0}$ the system returns to the ILC phase which shows a spread in the bifurcation diagram. In this case again the bifurcation diagram shows a hysteresis and allows us to obtain the boundaries separating the CLCo phase (black lines in figure 1(b)). In the bistable regime, CLCi and CLCo can be distinguished based on the amplitude values of the oscillators. For CLCi, both oscillators have the same amplitude (compare red lines in figures 2(c) and (d)) whereas for CLCo they have different amplitudes (compare black lines in figures 2(c) and (d)).

Noiseless CLC solution.- For the system with an inertia-like term, we can find the semi-analytic solution in the CLC regime assuming that the two oscillators have the same frequency but different amplitude and phase. Thus, using ${z}_{1}={r}_{1}{\text{e}}^{\text{i}({\Omega}t+{\phi }_{1})}$ and ${z}_{2}={r}_{2}{\text{e}}^{\text{i}({\Omega}t+{\phi }_{2})}$ in equation (5) with ξ₁(t) = ξ₂(t) = 0 and solving for r_i , Ω, and Δϕ = ϕ₁ − ϕ₂ we obtain,

$$\begin{equation}{\Omega}=\frac{{\omega }_{1}{r}_{1}^{2}+{\omega }_{2}{r}_{2}^{2}}{{r}_{1}^{2}+{r}_{2}^{2}},\hspace{25.0pt}\text{sin}({\Delta}\phi )=-\frac{\gamma {r}_{1}{r}_{2}{\Delta}\omega }{k({r}_{1}^{2}+{r}_{2}^{2})},\end{equation} \tag{ 6 }$$

$$\begin{equation}R-{r}_{1}^{2}+m{{\Omega}}^{2}+k\left(\frac{{r}_{2}}{{r}_{1}}\enspace \mathrm{cos}({\Delta}\phi )-1\right)=0,\end{equation} \tag{ 7 }$$

$$\begin{equation}R-{r}_{2}^{2}+m{{\Omega}}^{2}+k\left(\frac{{r}_{1}}{{r}_{2}}\enspace \mathrm{cos}({\Delta}\phi )-1\right)=0.\end{equation} \tag{ 8 }$$

Equations (7) and (8) should be solved simultaneously to obtain r_i . For the CLCo regime it is non-trivial to obtain the analytic solution and reduces to finding the roots of a quartic equation but for the CLCi regime since r₁ = r₂ = r we obtain r² = R + mΩ² + k[cos(Δϕ) − 1] with Ω = (ω₁ + ω₂)/2 and $\mathrm{cos}({\Delta}\phi )=\sqrt{1-{\gamma }^{2}{\Delta}{\omega }^{2}/4{k}^{2}}$. In the CLCi regime, the presence of an inertia-like term modifies only the radius (r > 0 since 4k² > γ²Δω² in the CLCi regime) of the limit cycle and not the frequency or phase difference. On the other hand in the CLCo regime the presence of an inertia-like term affects all parameters.

In the large coupling k limit, assuming that the radii of the two limit cycle oscillators are finite, equations (7) and (8) which need to be satisfied simultaneously reduce to

$$\begin{align*}\frac{{r}_{2}}{{r}_{1}}\enspace \mathrm{cos}({\Delta}\phi )-1& =0,\\ \frac{{r}_{1}}{{r}_{2}}\enspace \mathrm{cos}({\Delta}\phi )-1& =0,\end{align*}$$

which only permits one solution, i.e. the CLCi solution with r₁ = r₂, implying that the system would always end up coherently synchronized in phase for large couplings.

Noisy regime.- In presence of noise, i.e. at a finite temperature, the noise averaged asymptotic solution decays to a fixed point as seen in figures 1(g)–(j) (black and red solid lines). In this case, at the level of averages, we do not have a complex phase diagram and noise wipes out all signatures of the limit-cycle behavior and bistability in the asymptotic limit. The individual trajectories (grey lines in figures 1(g)–(j)) carry huge fluctuations and may carry important correlations which can be wiped out in the noise averaged quantities like ⟨z_i ⟩. Hence, these observables are not appropriate quantities to understand the behavior of the oscillators in different phases. In other words, the noisy ISLD does not exhibit different phases if we characterize its phase using the observables ⟨z_i ⟩ as one would normally do in the deterministic case.

It is worth noting here that unlike the equilibrium scenario wherein phases can be fully characterized by the partition function, in nonequilibrium such a universal function does not exist [54]. Hence, one set of observables showing a lackluster behavior does not necessarily mean that all observables will behave similarly. In the next section we will see how thermodynamic observables like heat and work retain information of the different phases and hence show universal behaviors depending on the phase of the oscillators.

In this section, we show how a single noisy inertial-like SLO can be obtained from the dynamics of a charged particle subject to a quartic potential and a constant magnetic field, and connected to a heat bath. This allows us to justify the form of the parameters entering the SLO equation, and inform our thermodynamic analysis in the next section.

Consider an electric point particle with charge q and mass m, with coordinates X and Y, velocities $\dot{X}$ and $\dot{Y}$, subject to an external force deriving from a potential U which we assume to be central in the plane (X, Y) (e.g. the horizontal components of a potential electric field), and a Lorentz force due to a constant magnetic field B normal to the plane (X, Y). We also assume the particle is in contact with a heat bath (e.g. a viscous fluid) of strength γ and temperature T, resulting in both a drag force and a random force. The equations of motion of this system are then

$$\begin{equation}m\ddot{X}=-{\partial }_{X}U+qB\dot{Y}-\gamma \dot{X}+{\xi }_{X},\end{equation} \tag{ 9 }$$

$$\begin{equation}m\ddot{Y}=-{\partial }_{Y}U-qB\dot{X}-\gamma \dot{Y}+{\xi }_{Y},\end{equation} \tag{ 10 }$$

where ξ_X (t) and ξ_Y (t) are Gaussian white noises with variance ⟨ξ_X (t)ξ_X (t')⟩ = ⟨ξ_Y (t)ξ_Y (t')⟩ = 2k_B Tγδ(t − t'). Note that the potential U is arbitrary for now, to be as general as possible, but we will need to specify it to a quartic function later in order to recover the Stuart–Landau equations.

As is well known, the magnetic force does no work, though it is not quite a potential force. Let us define half the cyclotron frequency of the particle: ω = qB/2m. The conservative part of the dynamics, obtained for γ = 0, has then two conserved quantities: (i) the energy

$$\begin{equation}h=\frac{m}{2}({\dot{X}}^{2}+{\dot{Y}}^{2})+U({X}^{2}+{Y}^{2}),\end{equation} \tag{ 11 }$$

and (ii) the third component of the angular momentum

$$\begin{equation}{L}_{3}={\overrightarrow{e}}_{3}\cdot \left[\overrightarrow{r}\times (m\dot{\overrightarrow{r}}+q\overrightarrow{B}\times \overrightarrow{r})\right]=m\omega ({X}^{2}+{Y}^{2})+m(X\dot{Y}-Y\dot{X}).\end{equation} \tag{ 12 }$$

Given the form of the coupling with the bath, which satisfies the Einstein relation with respect to the energy (and does not involve the angular momentum), the stationary distribution of this system is simply

$$\begin{equation}\rho (X,Y,\dot{X},\dot{Y})={\mathcal{N}}^{-1}\enspace \mathrm{exp}[-\beta h],\end{equation} \tag{ 13 }$$

with the natural normalisation $\mathcal{N}=\int \mathrm{d}X\int \mathrm{d}Y\int \mathrm{d}\dot{X}\int \mathrm{d}\dot{Y}\enspace \mathrm{exp}[-\beta h]$. Note that, even though L₃ is not relevant to our following analysis we give the relevant equations for the sake of completeness.

Defining Z = X + iY and ξ(t) = ξ^R (t) + iξ^I (t), the equations of motion (9) and (10) become

$$\begin{equation}m\ddot{Z}=-{\partial }_{\bar{Z}}U(\vert Z{\vert }^{2})-(\text{i}2m\omega +\gamma )\dot{Z}+\xi ,\end{equation} \tag{ 14 }$$

with $\bar{Z}$ being the complex conjugate of Z. Let us define new coordinates z = x + iy in the rotating frame through

$$\begin{equation}z={\mathrm{e}}^{\text{i}\omega t}Z.\end{equation} \tag{ 15 }$$

Placing ourselves in the rotating frame of angular velocity ω simplifies the velocity-dependent part of the dynamics, at the cost of an extra rotational force. Moreover, given that an isotopic Gaussian white noise is invariant under rotation, we deduce an equation for $\ddot{z}$:

$$\begin{equation}m\ddot{z}=-{\partial }_{\bar{z}}\left(U(\vert z{\vert }^{2})+m{\omega }^{2}\vert z{\vert }^{2}/2\right)+\text{i}\gamma \omega z-\gamma \dot{z}+\xi .\end{equation} \tag{ 16 }$$

If the potential U is chosen to be quartic in |z|, we see that this becomes the noisy Stuart–Landau equation for a single oscillator. Importantly, the rotational term is proportional to γ, i.e. it is a virtual force resulting from friction in a rotating frame, and vanishes for γ = 0, unlike the centrifugal force −mω² z. Note that, by construction, this system is in fact at equilibrium in a rotating frame, which does not affect its thermodynamics but gives its dynamics the appearance of non-equilibrium. All systems of this type are described by the so-called GENERIC formalism [37–39].

Let us therefore fix the value of the potential U from here on. For the sake of simplicity, we will choose it such that it cancels the centrifugal term mentioned above:

$$\begin{equation}U(\vert z{\vert }^{2})=\frac{\vert z{\vert }^{4}}{4}-(R+m{\omega }^{2})\frac{\vert z{\vert }^{2}}{2}.\end{equation} \tag{ 17 }$$

We thus recover the uncoupled inertial-like SLO equation (5) for k = 0.

Under these new coordinates, the energy becomes

$$\begin{equation}h=\frac{m}{2}\left[{(\dot{x}+\omega y)}^{2}+{(\dot{y}-\omega x)}^{2}\right]+U({x}^{2}+{y}^{2}),\end{equation} \tag{ 18 }$$

giving us the corresponding stationary distribution, equation (13), and the third component of the angular momentum simplifies to

$$\begin{equation}{L}_{3}=m(x\dot{y}-y\dot{x}).\end{equation} \tag{ 19 }$$

An important remark is that, although we have identified an energy (18) for the stationary distribution of a single SLO, this energy is not a Hamiltonian for the conservative part of the inertial-like Stuart–Landau equations, i.e. the equations of motion for γ → 0 do not derive from h in the usual way. The reason for this is that the coordinate transformation from Z to z is not canonical, and this explains why the single SLO does not appear to be detail-balanced even though it is.

As we mentioned earlier, the Einstein relation is verified between the noise and the dissipation, relative to an inverse temperature β and the Shannon entropy for independent particles. The heat exchanged with the reservoir is then unambiguously defined as the rate at which energy is lost by the system. Using the dynamics (14) of a single oscillator i to compute this rate, and allowing for different values of the frequencies by denoting them ω_i , we find

$$\begin{equation}\delta {Q}_{i}=-{\dot{h}}_{i}{\vert }_{k=0}=\gamma ({\dot{X}}_{i}^{2}+{\dot{Y}}_{i}^{2})-({\dot{X}}_{i}{\xi }_{i}^{R}+{\dot{Y}}_{i}{\xi }_{i}^{I}),\end{equation} \tag{ 20 }$$

where the second part averages out to 0 due to the noises. Since we dealt with only a single oscillator in the previous sub-section the quantities like energy h or phase space variables $\left\{\dot{x},x,\dot{y},y\right\}$ did not have a sub-index but from now on the sub-index i is the oscillator index. Note that both parts of the heat vanish for γ → 0, as expected. In terms of the rotating coordinates z_i , this becomes

$$\begin{equation*}\delta {Q}_{i}=\gamma \left[{({\dot{x}}_{i}+{\omega }_{i}{y}_{i})}^{2}+{({\dot{y}}_{i}-{\omega }_{i}{x}_{i})}^{2}\right]-\left[({\dot{x}}_{i}+{\omega }_{i}{y}_{i}){\xi }_{i}^{R}+({\dot{y}}_{i}-{\omega }_{i}{x}_{i}){\xi }_{i}^{I}\right].\end{equation*}$$

As a side note, the angular momentum exchanged by oscillator i with the reservoir is given by

$$\begin{equation}-{\dot{L}}_{3}=m\gamma ({X}_{i}{\dot{Y}}_{i}-{Y}_{i}{\dot{X}}_{i})+({X}_{i}{\xi }_{i}^{I}-{Y}_{i}{\xi }_{i}^{R})\end{equation} \tag{ 21 }$$

$$\begin{equation}=m\gamma ({x}_{i}{\dot{y}}_{i}-{y}_{i}{\dot{x}}_{i})-m\gamma \omega ({x}_{i}^{2}+{y}_{i}^{2})+({x}_{i}{\xi }_{i}^{I}-y{\xi }_{i}^{R}).\end{equation} \tag{ 22 }$$

where once again the noisy part averages out to 0 and both terms disappear for γ = 0.

Defining work is slightly less straightforward, as it requires us to make a meaningful choice of a reference potential meant to describe the system without external driving. Without a proper physical interpretation of the system and of the driving, any potential could be a candidate. In the present case, we have two interesting possibilities:

Case I: the first one is to consider the potential defined above for each oscillator separately, and define work as the energy variation in the system due to the inclusion of the interaction k. In this case, the undriven system is simply the two oscillators, each following equation (14) with the appropriate frequency ω_i . The driving consists in adding the terms proportional to k in equation (5). The work rate can then be computed by applying this driven dynamics to the reference potential, in order to estimate the rate of energy change, and extracting the term which is not heat, which by construction is the term proportional to k. Hence, we get

$$\begin{align}\delta {W}^{(0)}& ={\dot{h}}_{1}+{\dot{h}}_{2}+\delta {Q}_{1}+\delta {Q}_{2}\\ & =-k\left[({x}_{1}-{x}_{2})({\dot{x}}_{1}-{\dot{x}}_{2})+({y}_{1}-{y}_{2})({\dot{y}}_{1}-{\dot{y}}_{2})+{\Delta}\omega ({x}_{1}{y}_{2}-{y}_{1}{x}_{2})\right].\end{align} \tag{ 23 }$$

Case II: alternatively, we can consider the interaction term in equation (5) from a different reference system which is that of two charged particles coupled by a spring interaction k in their original coordinate frame. This means adding a spring potential of strength k in the original reference frame, so that the system has a total energy

$$\begin{equation}H={h}_{1}+{h}_{2}+\frac{k}{2}\left[{({X}_{2}-{X}_{1})}^{2}+{({Y}_{2}-{Y}_{1})}^{2}\right].\end{equation} \tag{ 24 }$$

In the rotating frames, this extra term becomes more complicated because of the different frequencies:

$$\begin{align}H& ={h}_{1}+{h}_{2}+\frac{k}{2}\left({x}_{1}^{2}+{y}_{1}^{2}+{x}_{2}^{2}+{y}_{2}^{2}-2({x}_{1}{x}_{2}+{y}_{1}{y}_{2})\mathrm{cos}\left({\Delta}\omega t\right)\right.\\ & \quad -\left.2({x}_{1}{y}_{2}-{y}_{1}{x}_{2})\mathrm{sin}\left({\Delta}\omega t\right)\right).\end{align} \tag{ 25 }$$

For the same reason, the interaction term in the reference equations of motion for z_i is time-dependent in the rotating frame if the two frequencies ω_i are different, so that replacing it by a time-independent k as seen in equation (5) amounts to having a work-performing external driving. In this case, when computing work in the same way as case I above, the first term of equation (23) disappears, and we are left with

$$\begin{equation}\delta {W}^{(k)}=\dot{H}+\delta {Q}_{1}+\delta {Q}_{2}=-k{\Delta}\omega ({x}_{1}{y}_{2}-{y}_{1}{x}_{2})\end{equation} \tag{ 26 }$$

which is surprisingly simple. Note that this case has the added advantage of producing no work in more situations (δW^(k) = 0 for ω₁ = ω₂ even if k ≠ 0), whereas in case I this would lead to work that averages out to 0 along any periodic trajectory. The definition of work we use is however ultimately a matter of choice, unless we have a physical reason to prefer one reference energy over the other. That being said, the difference between the two is a total time derivative, so that the choice does not affect the time-averaged output of the system.