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:

	$\displaystyle m\ddot{x}_{1}$	$\displaystyle=$	$\displaystyle(R-x_{1}^{2}-y_{1}^{2})~{}x_{1}-\gamma\omega_{1}y_{1}-k~{}(x_{1}-x_{2})-\gamma\dot{x}_{1}+\xi_{1}^{R},$		(1)
	$\displaystyle m\ddot{y}_{1}$	$\displaystyle=$	$\displaystyle(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},$		(2)
	$\displaystyle m\ddot{x}_{2}$	$\displaystyle=$	$\displaystyle(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},$		(3)
	$\displaystyle m\ddot{y}_{2}$	$\displaystyle=$	$\displaystyle(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}.$		(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 ($\gamma>>m$, $k\rightarrow 0$, and $\xi_{l}^{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 $\omega_{l}$ (scaled by the damping coefficient) and the oscillators are coupled via a spring with strength $k$. The parameter $\gamma$ represents the phenomenological Stokesian damping coefficient. The random noise $\xi_{j}^{R,I}$, $j=1,2$ are Gaussian with zero mean and $\langle\xi_{l}^{\alpha}(t)*\xi_{j}^{\beta}(t^{\prime})\rangle=2\gamma_{j}k_{B}T_{j}\delta(t-t^{\prime})\delta_{\alpha,\beta}\delta_{l,j}$, where $\langle\cdot\rangle$ 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, $\omega_{1}=1$ and $\omega_{2}=\omega_{1}+\Delta\omega$, unless mentioned otherwise. The temperatures $T_{1}$ and $T_{2}$ 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_{1}=x_{1}+iy_{1},z_{2}=x_{2}+iy_{2}$, which yields the following compact expressions:

	$\displaystyle m\ddot{z}_{1}$	$\displaystyle=$	$\displaystyle\left(R+i\gamma\omega_{1}-|z_{1}|^{2}\right)z_{1}-\gamma\dot{z}_{1}+k(z_{2}-z_{1})+\xi_{1},$
	$\displaystyle m\ddot{z}_{2}$	$\displaystyle=$	$\displaystyle\left(R+i\gamma\omega_{2}-|z_{2}|^{2}\right)z_{2}-\gamma\dot{z}_{2}+k(z_{1}-z_{2})+\xi_{2}.$		(5)

Above the complex noise $\xi_{l}=\xi_{l}^{R}+i\xi^{I}_{l}$. The above set of equations with complex variable $z_{l}$ converges to the standard Stuart-Landau dimer in the overdamped limit $\gamma>>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 (NEMS) 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 Refs. [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 inertial-like Stuart-Landau dimer 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 Figs. 1 a and b for the noiseless limit with small and large frequency difference $\Delta\omega$, 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 $\gamma=1$.
 
Noiseless without an inertia-like term.– In the case without an inertia-like term ($m=0$), there are three distinct phases, an incoherent limit cycle (ILC), amplitude death (AD), and an in-phase coherent limit cycle (CLCi). For $\gamma\Delta\omega<1$ the AD regime vanishes and we have an exceptional point that demarcate the ILC to CLCi transition. For $\gamma\Delta\omega>1$ the AD region emerges such that its area in the $k$-$\gamma\Delta\omega$ plane increases [48, 49]. In this overdamped case, the boundaries of the AD regime can be analytically predicted using linear stability analysis and for $\gamma\Delta\omega<1$ the critical coupling $k=k_{c}=\gamma\Delta\omega/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_{l}=0$), Fig. 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, Fig. 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 Fig. 1 d. This is known as the in-phase coherent limit cycle (CLCi) regime since the phase differences are smaller than $\pi/2$. As the $k$ increases further, the phase difference decreases and in the limit $k\rightarrow\infty$ the oscillators are perfectly in phase.
 
Noiseless with an inertia-like term: small frequency difference $\gamma\Delta\omega<1$.– Adding an inertia-like term to the system leads to rich limit cycle dynamics. If $m$ is sufficiently large with $\gamma\Delta\omega$ small, incoherent (ILC) and coherent (CLCi) limit cycles can coexist as seen in Fig. 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 Fig. 1 b) and the red line in Figs. 1 a and b that corresponds to an exceptional ‘line’ (see 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 $\gamma\Delta\omega=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-\delta k$ plus a small perturbation (forward direction, black lines in Fig. 2). In the ILC regime ($k<k_{c}^{f}$), the amplitude oscillates as seen by the spread of $|z_{i}|$ (black region hidden under red in Figs. 2 a and b). At the critical value of the coupling $k_{c}^{f}$ a small fluctuation of initial conditions triggers a discontinuous jump and the oscillators synchronize in phase (black line, hidden under red, in Figs. 2 a and b). In the reverse direction [53], (red lines in Fig. 2) for $k>k_{c}^{b}$ the oscillators remain sychronized in phase (CLCi) and we observe a single red line in Figs. 2 a and b. At the backward critical value of the coupling $k_{c}^{b}$ the system transits from CLCi to an ILC which is maintained for all $k<k_{c}^{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 Fig. 1 a, but also clearly shows the hysteric synchrony of incoherent and coherent states in the region between $k_{c}^{b}$ and $k_{c}^{f}$ [29, 30] that crucially depends on the choice of initial conditions indicating bistability.
 
Noiseless with an inertia-like term: large frequency difference $\gamma\Delta\omega>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 $\gamma\Delta\omega$ 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 $\pi/2$ (see Fig. 1 e). At small mass $m$, unlike the low frequency difference case, a new boundary emerges (solid black line in Fig. 1 b) that separates the ILC and CLCo regime. For large $m$, the boundary at the critical coupling $k=k_{c}=\gamma\Delta\omega/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 Fig. 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 $\gamma\Delta\omega=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_{c}^{0}$ and hence we see the spread reduce to a single line around $k\approx 1.3$. As $k$ increases further, the oscillators are synchronized (CLCo) up to a $k=k_{c}^{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 Figs. 2 c and d), the oscillators maintain an in-phase synchronous behaviour (CLCi) up to a different $k=k_{c}^{b}$ at which the system transits to being out-of-phase (CLCo). At $k_{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 Fig. 1 b). In the bistable regime, CLCi and CLCo can be distinguised based on the amplitude values of the oscillators. For CLCi, both oscillators have the same amplitude (compare red lines in Figs. 2 c and d) whereas for CLCo they have different amplitudes (compare black lines in Figs. 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}e^{i(\Omega t+\phi_{1})}$ and $z_{2}=r_{2}e^{i(\Omega t+\phi_{2})}$ in Eq. (2.1) with $\xi_{1}(t)=\xi_{2}(t)=0$ and solving for $r_{i}$, $\Omega$, and $\Delta\phi=\phi_{1}-\phi_{2}$ we obtain,

	$\displaystyle\Omega$	$\displaystyle=$	$\displaystyle\frac{\omega_{1}r_{1}^{2}+\omega_{2}r_{2}^{2}}{r_{1}^{2}+r_{2}^{2}},\quad\quad\sin(\Delta\phi)=-\frac{\gamma r_{1}r_{2}\Delta\omega}{k(r_{1}^{2}+r_{2}^{2})},$		(8)
			$\displaystyle R-r_{1}^{2}+m\Omega^{2}+k\left(\frac{r_{2}}{r_{1}}\cos(\Delta\phi)-1\right)=0,$
			$\displaystyle R-r_{2}^{2}+m\Omega^{2}+k\left(\frac{r_{1}}{r_{2}}\cos(\Delta\phi)-1\right)=0.$

Equations (8)-(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_{1}=r_{2}=r$ we obtain $r^{2}=R+m\Omega^{2}+k[\cos(\Delta\phi)-1]$ with $\Omega=(\omega_{1}+\omega_{2})/2$ and $\cos(\Delta\phi)=\sqrt{1-\gamma^{2}\Delta\omega^{2}/4k^{2}}$. In the CLCi regime, the presence of an inertia-like term modifies only the radius ($r>0$ since $4k^{2}>\gamma^{2}\Delta\omega^{2}$ 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 effects all parameters.

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

	$\displaystyle\frac{r_{2}}{r_{1}}\cos(\Delta\phi)-1$	$\displaystyle=$	$\displaystyle 0,$
	$\displaystyle\frac{r_{1}}{r_{2}}\cos(\Delta\phi)-1$	$\displaystyle=$	$\displaystyle 0,$

which only permits one solution, i.e., the CLCi solution with $r_{1}=r_{2}$, 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 Figs. 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 Figs. 1 g - j) carry huge fluctuations and may carry important correlations which can be wiped out in the noise averaged quantities like $\langle z_{i}\rangle$. 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 $\langle z_{i}\rangle$ 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 $\gamma$ and temperature $T$, resulting in both a drag force and a random force. The equations of motion of this system are then

	$\displaystyle m\ddot{X}$	$\displaystyle=$	$\displaystyle-\partial_{X}U+qB\dot{Y}-\gamma\dot{X}+\xi_{X},$		(9)
	$\displaystyle m\ddot{Y}$	$\displaystyle=$	$\displaystyle-\partial_{Y}U-qB\dot{X}-\gamma\dot{Y}+\xi_{Y},$		(10)

where $\xi_{X}(t)$ and $\xi_{Y}(t)$ are Gaussian white noises with variance $\langle\xi_{X}(t)\xi_{X}(t^{\prime})\rangle=\langle\xi_{Y}(t)\xi_{Y}(t^{\prime})\rangle=2k_{B}T\gamma\delta(t-t^{\prime})$. 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: $\omega=qB/2m$. The conservative part of the dynamics, obtained for $\gamma=0$, has then two conserved quantities: i) the energy

$$h=\frac{m}{2}(\dot{X}^{2}+\dot{Y}^{2})+U(X^{2}+Y^{2}),$$

(11)

and ii) the third component of the angular momentum

$$L_{3}=\vec{e}_{3}\cdot\left[\vec{r}\times(m\dot{\vec{r}}+q\vec{B}\times\vec{r})\right]=m\omega(X^{2}+Y^{2})+m(X\dot{Y}-Y\dot{X}).$$

(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

$$\rho(X,Y,\dot{X},\dot{Y})=\mathcal{N}^{-1}\exp[-\beta h],$$

(13)

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

Defining $Z=X+iY$ and $\xi(t)=\xi^{R}(t)+i\xi^{I}(t)$, the equations of motion (9) and (10) become

$$m\ddot{Z}=-\partial_{\bar{Z}}U(|Z|^{2})-(i2m\omega+\gamma)\dot{Z}+\xi,$$

(14)

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

$$z=\mathrm{e}^{i\omega t}Z.$$

(15)

Placing ourselves in the rotating frame of angular velocity $\omega$ 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}$:

$$m\ddot{z}=-\partial_{\bar{z}}\left(U(|z|^{2})+m\omega^{2}|z|^{2}/2\right)+i\gamma\omega z-\gamma\dot{z}+\xi.$$

(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 $\gamma$, i.e., it is a virtual force resulting from friction in a rotating frame, and vanishes for $\gamma=0$, unlike the centrifugal force $-m\omega^{2}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, 38, 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:

$$U(|z|^{2})=\frac{|z|^{4}}{4}-(R+m\omega^{2})\frac{|z|^{2}}{2}.$$

(17)

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

Under these new coordinates, the energy becomes

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

(18)

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

$$L_{3}=m(x\dot{y}-y\dot{x}).$$

(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 $\gamma\rightarrow 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 $\beta$ 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 $\omega_{i}$, we find

$$\delta Q_{i}=-\dot{h}_{i}|_{k=0}=\gamma(\dot{X}_{i}^{2}+\dot{Y}_{i}^{2})-(\dot{X}_{i}\xi_{i}^{R}+\dot{Y}_{i}\xi_{i}^{I})$$

(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 $\{\dot{x},x,\dot{y},y\}$ 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 $\gamma\rightarrow 0$, as expected. In terms of the rotating coordinates $z_{i}$, this becomes

$\displaystyle\delta Q_{i}$

$\displaystyle=$

$\displaystyle\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].$

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

	$\displaystyle-\dot{L}_{3}$	$\displaystyle=$	$\displaystyle m\gamma(X_{i}\dot{Y}_{i}-Y_{i}\dot{X}_{i})+(X_{i}\xi_{i}^{I}-Y_{i}\xi_{i}^{R})$		(21)
		$\displaystyle=$	$\displaystyle 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}).$		(22)

where once again the noisy part averages out to $0$ and both terms disappear for $\gamma=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 Eq. (14) with the appropriate frequency $\omega_{i}$. The driving consists in adding the terms proportional to $k$ in Eq. (2.1). 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

	$\displaystyle\delta W^{(0)}$	$\displaystyle=$	$\displaystyle\dot{h}_{1}+\dot{h}_{2}+\delta Q_{1}+\delta Q_{2}$
		$\displaystyle=$	$\displaystyle-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].$

Case II: Alternatively, we can consider the interaction term in Eq. (2.1) 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

$\displaystyle H$

$\displaystyle=$

$\displaystyle h_{1}+h_{2}+\frac{k}{2}\left[(X_{2}-X_{1})^{2}+(Y_{2}-Y_{1})^{2}\right].$

(24)

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

	$\displaystyle H$	$\displaystyle=$	$\displaystyle h_{1}+h_{2}+\frac{k}{2}\Big{(}x_{1}^{2}+y_{1}^{2}+x_{2}^{2}+y_{2}^{2}-2(x_{1}x_{2}+y_{1}y_{2})\cos\left(\Delta\omega t\right)$		(25)
			$\displaystyle-2(x_{1}y_{2}-y_{1}x_{2})\sin\left(\Delta\omega t\right)\Big{)}.$

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 $\omega_{i}$ are different, so that replacing it by a time-independent $k$ as seen in Eq. (2.1) 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 Eq. (3.2) disappears, and we are left with

$$\delta W^{(k)}=\dot{H}+\delta Q_{1}+\delta Q_{2}=-k\Delta\omega(x_{1}y_{2}-y_{1}x_{2})$$

(26)

which is surprisingly simple. Note that this case has the added advantage of producing no work in more situations ($\delta W^{(k)}=0$ for $\omega_{1}=\omega_{2}$ even if $k\neq 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.