Stochastic thermodynamics of inertial-like Stuart-Landau dimer

Stuart-Landau limit-cycle oscillators are a paradigm in the study of coherent and incoherent limit cycles. In this work, we generalize the standard Stuart-Landau dimer model to include effects due to an inertia-like term and noise and study its dynamics and stochastic thermodynamics. In the absence of noise (zero-temperature limit), the dynamics show the emergence of a new bistable phase where coherent and incoherent limit cycles coexist. At finite temperatures, we develop a stochastic thermodynamic framework based on the dynamics of a charged particle in a magnetic field to identify physically meaningful heat and work. The stochastic system no longer exhibits the bistable phase but the thermodynamic observables, such as work, exhibit bistability in the temporally metastable regime. We demonstrate that the inertial-like Stuart-Landau dimer operates like a machine, reliably outputting the most work when the oscillators coherently synchronize and unreliable with minimum work output when the oscillators are incoherent. Overall, our results show the importance of coherent synchronization within the working substance in the operation of a thermal machine.


Introduction
Macroscopic heat engines that convert wasteful heat into useful mechanical work have been a cornerstone in establishing the laws of thermodynamics [1]. Due to the technological advancements over the past several decades, there is a renewed interest in meso-and micro-scopic engines [2,3,4,5]. In these nanoscale machines, which are connected to macroscopic heat baths, thermal fluctuations that randomize the dynamics of the small working substance are no longer negligible [6]. Such small systems thus, can not be described by traditional thermodynamics and has led to the advent of stochastic thermodynamics [7,8,9] that provides a rigorous mathematical framework to define not only the average thermodynamic quantities like heat and work, but also their distributions.
Nanoscale heat engines have been extensively studied exploring the roles of external time dependent driving [10,11], active dissipation [12,13,14,15], particleparticle interaction [16], network topology [17], collective phenomenon [18,19,20], etc., to improve the performance metrics of the engine. We focus on the fascinating phenomenon of synchronization [21] and explore whether the working substance comprising of synchronizing particles can make better machines? We consider the paradigmatic model, Stuart-Landau [22,23] dimer, used for understanding collective behaviours such as synchronisation and oscillation suppression in coupled nonlinear oscillators [24,25,26,27,28], and extend this overdamped-deterministic model to the underdamped-stochastic regime, incorporating the effects of inertia and temperature. Surprisingly, even though the inertial-like stochastic Stuart-Landau dimer has never been explored, its famous counterpart the Kuramoto model has been extensively studied. The inertial stochastic Kuramoto model is used to describe synchronisation of fireflies [29,30], disordered Josephson Junction arrays [31], power grids [32], and has also been used to explore nonequilibrium phase transitions [33,34].
In this work, our first main aim is to understand the effects of the inertial-like term on the phase diagram for the nonlinear oscillator model and explore how synchronisation affects the work and heat exchanges which characterise it as a thermodynamic machine. Similar to the inertial Kuramoto model [35,36], we find that due to the inertia-like term the zero-temperature phase diagram develops bistability leading to a coexistence of coherent and incoherent limit cycles. Secondly, we build a stochastic thermodynamics framework and show that the thermodynamic observables (heat and work) depend strongly on the physical interpretation of the dynamical equations being analyzed. In other words, given the dynamical equations there does not exist a unique mathematical framework to obtain physically meaningful heat and work. Specifically, the definitions of thermodynamic quantities depend on how one segregates the various parameters (system or reservoir parameters) appearing in the dynamical equations and if the dynamics has the correct time-reversal symmetry [37,38,39]. The physically meaningful thermodynamics depends also on whether the system is described by underdamped or overdamped dynamics [40], whether in a Newtonian or accelerating frame [41], and whether any of the forces are pseudovectors (such as magnetic fields) [42].
The paper is organised as follows: In Sec. 2, we start with defining the stochastic inertial-like Stuart-Landau dimer, and we explore its deterministic phase diagram. This will guide us while studying the stochastic version of the model by the introducing the thermal noise. In Sec. 3, we show that, though the inertial-like Stuart-Landau oscillator might appear to be driven by non-conservative environmental forces (involving work production even without dimer interaction), it can in fact be seen as an equilibrium charged particle in a magnetic field, described in a rotating reference frame; this gives us a well-motivated physical frame for defining work and heat in the interacting case. In Sec. 4, we numerically study the behavior of our model in terms of average work in different bistable regimes, for this we take a cue from the deterministic case. We also study the full probability distributions and the reliability of our system. We conclude in Sec. 5.

Stuart-Landau Dimer
The Stuart-Landau oscillator (SLO) is a prototypical system exhibiting Hopf bifurcation and limit cycle oscillations that can reveal universal features of many practical systems [22,23,24,25]. In this section, we introduce and study the inertial-like Stuart-Landau dimer, comprising of two coupled SLOs, and investigate the effects of the inertial-like term on its phase diagram. The specific parameters entering the equations will be physically justified in the next section.

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: 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 overdampeduncoupled-noiseless limit (γ >> m, k → 0, and ξ R(I) l (t) = 0) the oscillators settle in their own limit cycle with radius √ R for R > 0 and to a fixed point if R < 0. Each twodimensional oscillator is governed by a frequency ω l (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 ξ R,I j , j = 1, 2 are Gaussian with zero mean and ξ α l (t) * ξ β j (t ) = 2γ j k B T j δ(t − t )δ α,β δ l,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 = 1 and ω 2 = ω 1 + ∆ω, 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: Above the complex noise ξ l = ξ R l +iξ I l . The above set of equations with complex variable z l 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 (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.

Phase diagram and temporal behavior
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 ∆ω, 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.
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 γ∆ω < 1 the AD regime vanishes and we have an exceptional point that demarcate the ILC to CLCi transition. For γ∆ω > 1  ). ILC, CLCi, CLCo, and AD represent incoherent limit cycle, in-phase coherent limit cycle, out of phase coherent limit cycle, and amplitude death. The coexistence of incoherent and coherent limit cycles is shaded in green (bistable I) whereas the coexistence of different types of coherent limit cycles is in blue (bistable II). Trajectories on the complex plane and time series of the two oscillators when (k, m) = (0.5, 0.5) (panels c and g: ILC), (5.0, 0.5) (panels d and h: CLCi), (1.5, 0.5) (panels e and i: CLCo), and (2.0, 0.1) (panels f and j: AD). Panels (c-f) are the noiseless limit with black and red colors representing the first and second oscillator, respectively. The real parts Re(z) are the dashed lines in the temporal plots and the amplitudes |z| are represented by the solid lines. Panels (g-j) are the noisy oscillators with grey color being 10 selected trajectories of Re(z 1 ) of the first oscillator. The dashed black (red) ( Re(z 1 ) ) and solid black (red) (| z 1 |) represent trajectories of the first (second) oscillator averaged over 1000 realizations. Throughout we have set γ = 1 and the points marked c, d, e, and f in blue in panel b correspond to parameter values used in panels c, d, e, and f, respectively. The orange and red boundaries in panels a and b are determined analytically from linear stability analysis as shown in Appendix A. The black boundaries are numerically obtained using the bifurcation diagrams similar to Fig. 2. Temperature used in panels g, h, i is K B T L = k B T R = 4.0 and in j is 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 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 be-tween 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 π/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 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 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 Fig. 2). In the ILC regime (k < k f c ), 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 f c 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 b c 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 b c the system transits from CLCi to an ILC which is maintained for all k < k b c . 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 b c and k f c [29,30] that crucially depends on the choice of initial conditions indicating bistability.
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 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 = γ∆ω/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 γ∆ω = 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 0 c 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 f c . 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 b c at which the system transits to being out-of-phase (CLCo). At k 0 c 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). 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(Ωt+φ 1 ) and z 2 = r 2 e i(Ωt+φ 2 ) in Eq. (5) with ξ 1 (t) = ξ 2 (t) = 0 and solving for r i , Ω, and ∆φ = φ 1 − φ 2 we obtain, Equations (7)- (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Ω 2 + k[cos(∆φ) − 1] with Ω = (ω 1 + ω 2 )/2 and cos(∆φ) = 1 − γ 2 ∆ω 2 /4k 2 .
In the CLCi regime, the presence of an inertia-like term modifies only the radius (r > 0 since 4k 2 > γ 2 ∆ω 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. (7) and (8) which need to be satisfied simultaneously reduce to r 2 r 1 cos(∆φ) − 1 = 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 gj (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 gj) 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.

Stochastic Thermodynamics
We are now interested in studying the dynamics of two coupled inertial-like Stuart-Landau oscillators in presence of stochastic forcing. Similar study was done in case of single Duffing oscillator using the Volterra series approach in [55], where the authors were mainly interested in developing an analytical formulation for the time series analysis of a simplest bistable system. We, however, want to understand the interplay of noise and synchronisation in our model. A non-linear system is often influenced due to its contact with an environment e.g. a thermal bath. Hence, studying the thermodynamic quantities like heat and work and their effect on the behaviour of the system becomes important. We are interested in understanding this interplay within the stochastic thermodynamics framework [7,8,9].
In order to build a thermodynamic picture of this model, we need to be able to identify meaningful physical quantities such as heat and work, as well as the energy of the system as a function of its coordinates. For an equilibrium system coupled with a single bath and an external forcing through a conservative interaction, those quantities would be straightforward to define: the energy could be identified in the logarithm of the equilibrium distribution of the system with constant external force, and subsequently the heat and work could be defined as the changes of energy of the system related, respectively, to stochastic transitions, and to external parameter variations. Our first result in this section will be to show that, perhaps surprisingly, the ISLD does in fact fit in that framework, with each separate dimer being equivalent to a charged particle in a heat bath.

Inertial-like SLO as a charged particle
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Ẋ andẎ , 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 where ξ X (t) and ξ Y (t) are Gaussian white noises with variance ξ X (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 and ii) the third component of the angular momentum 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 with the natural normalisation N = dX dY dẊ dẎ exp[−β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 ξ(t) = ξ R (t) + iξ I (t), the equations of motion (9) and (10) become withZ being the complex conjugate of Z. Let us define new coordinates z = x + iy in the rotating frame through Placing ourselves in the rotating frame of angular velocity ω simplifies the velocitydependent 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 forz: 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ω 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: We thus recover the uncoupled inertial-like SLO equation (5) for k = 0. Under these new coordinates, the energy becomes giving us the corresponding stationary distribution, Eq. (13), and the third component of the angular momentum simplifies to 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.

Heat and work
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 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 {ẋ, x,ẏ, 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 γ → 0, as expected. In terms of the rotating coordinates z i , this becomes As a side note, the angular momentum exchanged by oscillator i with the reservoir is given by 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 Eq. (14) with the appropriate frequency ω i . The driving consists in adding the terms proportional to k in Eq. (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 Case II: Alternatively, we can consider the interaction term in Eq. (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 In the rotating frames, this extra term becomes more complicated because of the different frequencies: 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 Eq. (5) amounts to having a workperforming external driving. In this case, when computing work in the same way as case I above, the first term of Eq. (23) disappears, and we are left with which is surprisingly simple. Note that this case has the added advantage of producing no work in more situations (δW (k) = 0 for ω 1 = ω 2 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.

Numerical Analysis of Thermodynamic Work
We first analyze the zero temperature (T = 0) noiseless limit. In this regime, both definitions of work [Eqs. (23) and (26)] for CLC and AD regions coincide with Above we have used z j = r j e i(Ωt+φ j ) with the parameters r j obtained using Eqs. (7)- (8).
With the damping coefficient γ > 0 in the CLC regime the system always acts like a machine [18] with δW < 0 whereas no work is done in the AD regime since the system relaxes to a fixed point solution at the origin. Specifically, for the CLCi regime the average work can be computed analytically using r 1 = r 2 = r provided below Eq. (8).
For the general ILC regime wherein the amplitudes can be time dependent it is less clear whether the system works as a machine or a dissipator. In the noisy regime, we integrate Eqs. (1) -(4) numerically to study the thermodynamic response of the system. Simulations are performed using a modified velocity Verlet algorithm by Ermak [56]. The algorithm was historically introduced to include hydrodynamic interactions in Brownian dynamics simulations, but for our purpose even the standard fourth order Runge-Kutta method would suffice. Thus, we now look at the average work W in different bifurcation regimes described in Fig. 2. The system is first evolved for a transient time t = 20 (see Figs. 1 gj to gauge an estimate of the transient time). The average work at a specific k is calculated after averaging over a duration of 8 × 10 6 with the incremental time step dt ∼ 10 −3 . The system is initiated in the same forward and reverse direction initial conditions we used for the noiseless case in Fig. 2. The time averaged work W obtained from Eq. (23) at different temperatures is shown in Fig. 3.
At extremely low temperatures (black lines) in Fig. 3 the simulation time is not long enough for the system to relax to its asymptotic state. Hence, in this case the system remains in a metastable state for extremely long times that mimics the noiseless system perfectly. In the metastable regime, the average work shows bistability with the forward (solid black lines) and backward (dashed black lines) initial condition results being distinctly different between the critical backward k b c and critical forward k f c couplings. In the CLC regime, the work output ( W < 0) of the machine matches perfectly with the analytic results obtained in the noiseless limit, Eq. (27).
At moderate and high temperatures, we reach the asymptotic state in our simulation time and in this case the bistability disappears, i.e., the solid and dashed (blue and green) lines in Fig. 3 perfectly overlap. Thus, even though the system is in a unique steady state, the thermodynamic variable, work, retains critical information about the underlying phases obtained in the noiseless limit. Specifically, when the oscillators are incoherently synchronized the average work output is the minimum. Whereas, the magnitude of work output for the in-phase synchronized regime (k > k f c ) is always larger than the incoherently or the out-of-phase coherently synchronized oscillators. Moreover, even the reliability [57,58], that measures how much the average dominates over the fluctuations, is the highest when the oscillators of the working substance synchronize in phase. In this phase, both work output and reliability decrease as a function of temperature indicating that the desirable machine is only obtained in the low temperature limit, concurrent with the notion of inphase coherent synchronization. The reliability of the machine is always larger when the oscillators are coherently synchronized. At moderate temperatures k B T 1 = k B T 2 = 1, a sudden dip in the reliability is observed in the bistable II (CLCi + CLCo) regime around k ≈ 2.5 ( Fig. 3 d). At these temperature values, the work distribution becomes bimodal with the position of the peaks being determined by the deterministic CLCi and CLCo work given by Eq. (27) [see description below on the distributions and Fig. 4 c].
For the given choice of parameters, since the deterministic work values for CLCi and CLCo are well separated, the variance for the moderate temperature is large giving rise to a low reliability. Despite these occurrences that depend on the parameter choice, overall an in-phase synchronizing working substance always leads to a powerful and reliable machine. As a side remark, we would like to note that since we evaluate the time-averaged work output, our different definitions of work, Eqs. (23) and (26), match each other exactly, as expected. Another question of interest, when looking at the inertial-like Stuart-Landau dimer as an engine, is that of efficiency. At first glance, given that the system works with a single temperature, using the first law of thermodynamics over a cycle (i.e. any trajectory of the system in its stationary or periodic regimes) implies that the heat absorbed is exactly converted to work, leading to unit efficiency. That being said, the concept of efficiency relies on a meaningful definition of absorbed heat and of useful work, which are unambiguous in the context of traditional heat engines, but less so for microscopic machines. We will therefore leave this question for more concrete situations, Panel a is within the ILC region with k = 1.0, b is for the CLCo region with k = 1.6, and c represents the bistable (CLCi + CLCo) region with k = 2.6. Panel d is with m = 2.0, ∆ω = 0.5, and temperatures k B T 1 = k B T 2 = 10 −4 (black), 0.5 (green), and 1.0 (blue) in the bistable (ILC + CLCi) region with k = 0.3. Solid and dashed lines represent probability distributions of work W obtained from the forward and backward initial conditions, respectively. Our low temperature results (black lines) are in the metastable regime in which the bistability is clearly reflected in the distinct distributions obtained for different initial conditions (panels c and d). For the green and blue lines, the solid and dashed lines perfectly overlap. and simply note that, according to the noiseless expressions (27) that we obtained for the average work, this system operates on average as a useful machine in at least one regime, namely CLC. The stochastic thermodynamics framework described in Sec. 3.2 allows us to go beyond the averages and calculate the full probability density function (PDF) of work as evaluated numerically in Fig. 4. The system is evolved for a transient time (t = 20) as described earlier and then the accumulated stochastic work is evaluated over a duration of period τ = 40, which spans over multiple periods in the ILC regime for very low temperatures. This is repeated for N = 2 × 10 5 times, for a total duration of N τ , giving us the N stochastic work values to form the distribution for particular values of the trap strength k and temperatures. At small temperatures (black lines in Fig. 4), our distributions are obtained in the metastable regime and can be fully explained using the noiseless solutions. In this regime, the noise only perturbs the work distribution slightly and gives it a finite width proportional to the strength of the noise, as expected. In cases wherein the CLC solutions hold, since the work (noiseless asymptotic limit) only depends on the time-independent limit cycle radii r 1 and r 2 , Eq. (27), the distributions are peaked at W = δW CLC = −γ∆ω 2 r 2 1 r 2 2 /(r 2 1 + r 2 2 ). For the ILC regimes, the broader distributions observed in panels a and d result from the fact that instantaneous work is not constant along the ILCs, so that the average work over a time-span τ will depend on the initial condition on the cycle. The range of values taken by W will be of order τ −1 and will become negligible for τ large enough, though still noticeable if W is small, as is the case in panel d.
At moderate temperatures (green curves in Fig. 4), the asymptotic distributions are bimodal, which is a signature of the underlying noiseless bistability. In the noiseless stable regimes, ILC and CLCo regime (Figs. 4 a and b) the work distributions broaden with temperature as expected, but in the bistable regimes, Figs. 4 c and d, the distributions become bimodal. In the bistable II regime (panel c), wherein the two different limit cycles coexist (CLCi + CLCo), the system settles in the CLCi often indicated by the slightly higher maximum of the distribution centered around the noiseless CLCi distribution (dashed black line). The CLCo solution, which was stable in the noiseless limit, is no longer equally probable to the CLCi solution, but influences the distribution strongly making it bimodal. For the bistable I regime (panel d), the system again tends to often be trapped in the CLCi solution. The ILC solution, which was stable in the noiseless limit, influences the overall distribution and leads to the overall distribution being bimodal. As temperature increases further (blue curves) all signatures of bistability are washed away and the variance of most observables including work scale with the temperature as can be observed from Fig. 4.

Conclusions
To conclude we studied the stochastic inertial-like Stuart-Landau dimer, whose coordinates x i and y i are independent of each other and represent the position of the oscillators. We extensively analysed our model in the zero temperature limit to uncover the rich phase space structure, which is not observed in the non-inertial system. We characterized various phases depending on the limit cycle behavior of the oscillators, namely, incoherent (ILC), in-phase coherent (CLCi), out-of-phase coherent (CLCo), and amplitude death (AD). Interestingly, we found the existence of bistable phases in which different limit cycles can coexist. Using a combination of analytic (linear stability analysis) and numeric (bifurcation diagrams) tools we obtained the phase boundaries that demarcate the different co-existing phases. At finite temperatures, the system always relaxed to a unique fixed point solution, erasing all information about the underlying phases.
Furthermore, we developed a physically meaningful stochastic thermodynamics framework that required a frame transformation due to which the Stuart-Landau oscillator transformed to a magnetic charged particle in a quartic potential. Depending on the physical interpretation of the interaction between the two magnetic charged particles, stochastic heat and work were identified. More specifically, we considered two meaningful equilibrium systems of which the inertial-like Stuart-Landau dimer is then a driven version: one where the two charged particles are independent, and another where they interact through a spring potential.
In the zero temperature regime, wherein the phase diagram displayed distinct phases analytic expressions for work were obtained in the CLC and AD regime. Surprisingly, the system output work like a machine in this regime and the behavior persisted for higher temperatures. In the ILC regime, although we were not able to prove this analytically, our extensive numerical simulations showed that the dimer continues to act as a machine albeit with a lower work output as compared to the in-phase coherently synchronized oscillators (CLCi). Other performance metrics, such as the reliability, corroborated that when the working substance is in-phase coherently synchronized we have the most desirable machine with a high reliability.
The signatures of the underlying bistability persisted for the averages in the temporally metastable regime at low temperatures. This led to a hysteresis curve for the average work output as the interaction between the oscillators was varied. Although such signatures were completely wiped out once the system reached the asymptotic state, but the distribution of work continued to display remnants of the underlying bistability. This was observed in the bimodal structure of the work distributions at moderate temperatures, but in the stable asymptotic limit. At extremely high temperatures, the thermal fluctuations dominated all processes and all signatures of bistability either in the averages or in the distributions disappeared.
Our work highlights the importance of a synchronizing working substance in the performance of a machine. It not only shows the effect of phase transitions on thermodynamic observable, such as work, but also stresses that the stochastic thermodynamic framework need not be necessarily unique given the dynamical equations. In other words, stochastic thermodynamics should be a physically inspired and not a mathematically constructed framework using the equations of motion that govern the dynamics. Although, we considered a simple dimer in this study, extending our system to a network of oscillators to investigate the effect of long-range interactions [59], the influence of network topology on thermodynamic observables, or to understand how biological systems make use of synchronization [20], would be fascinating future directions.