Current-mediated synchronization of a pair of beating non-identical flagella

The basic phenomenology of experimentally observed synchronization (i.e., a stochastic phase locking) of identical, beating flagella of a biflagellate alga is known to be captured well by a minimal model describing the dynamics of coupled, limit-cycle, noisy oscillators (known as the noisy Kuramoto model). As demonstrated experimentally, the amplitudes of the noise terms therein, which stem from fluctuations of the rotary motors, depend on the flagella length. Here we address the conceptually important question which kind of synchrony occurs if the two flagella have different lengths such that the noises acting on each of them have different amplitudes. On the basis of a minimal model, too, we show that a different kind of synchrony emerges, and here it is mediated by a current carrying, steady-state; it manifests itself via correlated"drifts"of phases. We quantify such a synchronization mechanism in terms of appropriate order parameters $Q$ and $Q_{\cal S}$ - for an ensemble of trajectories and for a single realization of noises of duration ${\cal S}$, respectively. Via numerical simulations we show that both approaches become identical for long observation times ${\cal S}$. This reveals an ergodic behavior and implies that a single-realization order parameter $Q_{\cal S}$ is suitable for experimental analysis for which ensemble averaging is not always possible.

(1), the dimensionless functions θ 1,2 (t) ∈ (−1/2, 1/2) are the two phases mentioned above, the dot denotes the time derivative, ν 1 and ν 2 are the natural frequencies (Hz) of the flagella 1 and 2, respec-tively, and ε (Hz) is the amplitude of the coupling between the flagella. This phenomenological parameter accounts for the fact that fluid flows driven by beating flagella provide a hydrodynamic coupling between the latter. Lastly, ζ 1,2 (t) are delta-correlated, Gaussian white noises with zero mean and identical covariances [15,16]: ζ i (t)ζ j (t ) = 2T eff δ i,j δ(t − t ), where the bar denotes the average over realizations of the noises; i, j ∈ 1, 2; δ i,j is the Kronecker symbol while T eff (Hz) can be considered as an effective "temperature", because it defines the amplitude of the noise terms. Importantly, the major contribution to T eff stems from the fluctuations of the rotary motors of flagella [14,17,18]. Indeed, these active fluctuations are several orders of magnitude larger than the thermal noise [14,18], so that the flagella are not in thermal equilibrium with its bath.
We note that the sine terms in Eq. (1), which link the time evolutions of θ 1 (t) and θ 2 (t), describe the actual coupling due to hydrodynamic interactions between the two flagella [6][7][8][9][10][11][12][13] only in an effective way. However, up to now no other good and justified alternative to such a phenomenological description has emerged. In particular, explicit results presented in, e.g., Ref. [13], are based on the assumption that the distance between the two flagella is much larger than their length, which is not the case for the system studied in Ref. [14]. We also note that in more complex situations of unicellular algal species bearing multiple flagella, both the effective coupling introduced in Ref. [14] and the results in Refs. [6][7][8][9][10][11][12][13] may turn out to be insufficient to describe properly all facets of the synchronized behavior: a different synchronization scenario may be realized which is provided, e.g., by contractile fibers of the basal apparatus [19,20].
Solving Eq. (1) numerically for given realizations of the noise terms, Goldstein et al. [14] have found consistency between ∆ t evolving according to Eq. (1) and the experimentally observed behavior of the flagella of the biflagellate alga [1], i.e., the calculated trajectories of the phase difference exhibit essentially the same noisy synchronization interrupted by occasional phase slips.
The comparison with experimental data has facilitated to identify the physically relevant values of the parameters entering the effective Langevin equations. In particular, for flagella of length l 12 µm, observations based on the dynamics of 21 individuals and the comparison with the time series for ∆ t , spanning over an interval of 10 2 seconds (i.e., containing several thousands of beats), have shown that ε lies within the range 0.14 − 0.7 Hz and that the effective temperature T eff is within the range 0.05 − 0.28 Hz, while ν = (ν 1 + ν 2 )/2 47 Hz and δν/ν = |ν 1 − ν 2 |/ν 0.004. The experimentally obtained values of ε appear to be in line with the theoretical prediction in Ref. [13]. Importantly, the results of Ref. [14] have emphasized for the first time the essential role played by the biochemical noise in the dynamics of eukaryotic flagella as manifested by its realization-torealization fluctuations.
A more sophisticated experimental analysis has been performed in Ref. [21], which is focused on the dependence of the coupling parameter ε and of the effective temperature T eff on the length l of the flagella. This enhanced experimental analysis took advantage of the ability of the Chlamydomonas reinhardtii alga to shed its flagella and to regrow them after a deflagellation has occurred [22]. The flagella of length l 10.82 µm have been first clipped by a micropipette [21] and then left to regrow. Within 90 minutes the flagella reached the length l 11.48 µm, which, surprisingly, exceeded the original one. The dynamics of the slowly growing flagella has been recorded every ten minutes within time intervals two minutes long (within which the length of the flagella did not appreciably change).
From 19 such experiments it was inferred [21] that, for progressively longer flagella, the periods of the synchronous beating of both flagella become more pronounced. Analyzing the data, it was shown that, as the flagella grow, the beating frequency ν decreases ∝ 1/l, implying that the beating mechanism operates at constant power output per length (see Ref. [21]). The coupling parameter ε turned out to be linearly proportional to l, which explains the trend for a progressive increase of the synchronization periods. Therefore, the proportionality ε ∝ l is in agreement with the elastohydrodynamic scaling ε ∝ ν 2 l 3 as predicted in Ref. [13]. Lastly, a variation of T eff with l has been observed. In particular, for l 6 µm, T eff was found to lie within the range 0.06−0.09 Hz, while for l 8 µm it lies within the range 0.04 − 0.06 Hz, which is not overlapping with the previous interval.
For larger values of l, T eff was shown to saturate at a constant value of the order of 0.04 Hz. Therefore, T eff evidently depends on l, at least for sufficiently short flagella, such that it is larger for shorter flagella. In view of the active nature of the fluctuations of the rotary motors, this is in line with the intuitively expected behavior [23].
Once noise appears to be a physically relevant parameter, it is natural to explore a wider range of possible effects. In this sense the conceptually important question arises what kind of synchronization, if any, may take place in situations in which the length of the two flagella differ. Such a situation may apparently be realized experimentally by amputating just one flagellum of the biflagellate alga, and by leaving the second one intact, as described in Ref. [22]. We note that in this case the regeneration scenario is more complicated, as compared to the case when both flagella are removed (see Ref. [21]). Here, the intact flagellum first shortens linearly in time while the amputated one regenerates. This way the two flagella attain an equal, intermediate length. Then both grow and eventually approach their initial length at the same rate. However, the time required to reach an equal intermediate length can be quite long, i.e., 20 to 40 minutes [22]. It can become even longer if certain chemicals (e.g., colchicine) are added after a deflagellation, which inhibit the regeneration process [22,24]. Therefore, there is a time window in which both flagella have distinctly different lengths. Within this time window, the coupled phases θ 1 (t) and θ 2 (t) will undergo a stochastic evolution -each at its own temperature T (1) eff or T (2) eff , respectively. Such a system is no longer characterized by a unique effective temperature T eff . One expects that the Fokker-Planck equation [16] associated with the Langevin equations (1) will have a non-trivial, current-carrying steady-state solution. In the following we shall refer to the case of unequal temperatures as an out-of-equilibrium case, keeping in mind, of course, that the original physical system is not in equilibrium with its bath, even for T Viewed from a different perspective, which is perhaps equally important due to certain other applications (see Ref. [26] for a discussion) we note that the minimal model in Eq. (1) with a unique effective temperature T eff represents the so-called Sakaguchi model [27] which is a noisy version of the celebrated Kuramoto model of coupled oscillators (see, e.g., Refs. [2][3][4][5]). Its generalization to the case of two different temperatures [28] emerges naturally within the present context of the synchronization of beating, non-identical flagella. We note parenthetically that recently the stochastic evolution of systems with several temperatures was intensively studied and a wealth of interesting out-of-equilibrium phenomena has been predicted (see, e.g., Refs. [29][30][31][32][33][34][35][36] and references therein). To the best of our knowledge, the issue of synchronization under out-of-equilibrium conditions in general, and in a system with two degrees of freedom exposed to two different effective temperatures in particular, has not yet been addressed. Inter alia, this motivates our quest for synchrony in a minimal model with two different effective temperatures.
Here, we focus on the stochastic evolution of the phases θ 1 (t) and θ 2 (t) of two coupled oscillators, which obeys the minimal model in Eq. (1) with the covariance functions of the noise terms of the form where T (1) eff and T (2) eff are, in general, not equal. For simplicity, we assume that the natural frequencies of both oscillators are the same, ν 1 = ν 2 = ν. On one hand, this assumption appears to be justified because the experimentally observed difference of the natural frequencies is indeed rather small (see above) [14], so that in a first approximation it can be neglected. On the other hand, this assumption allows us to disentangle the effects of an out-of-equilibrium active noise from the effects caused by a possible, albeit small, difference of the natural frequencies ν 1 and ν 2 .
We demonstrate, both analytically and numerically, that in such a system an emerging steady-state is characterized by a nonzero current j(θ 1 , θ 2 ) in the frame of ref- eff ). This current, which is the same for both phases, with an amplitude depending on the instantaneous values of θ 1 (t) and θ 2 (t), sustains a synchronized time evolution of the ratesθ 1 (t) andθ 2 (t) at which the phases change, i.e., it produces correlated drifts of phases. At the same time, we realize that the stochastic phase locking seen in Ref. [14] degrades for unequal effective temperatures (see Fig. 1) and is weakest in the case that one of the effective temperatures equals zero, i.e., that the corresponding rotary motor does not fluctuate. In order to quantify the degree of synchronization, based on the correlated drifts of phases, we define a characteristic order parameter Q, which measures the relative amount of the novel, out-of-equilibrium synchronization mechanism and vanishes if the effective temperatures of the noise terms become equal. This definition of Q is based on the explicit expression derived here for the steady-state current j(θ 1 , θ 2 ) and hence, it represents a property averaged over the statistical ensemble of the trajectories of θ 1 (t) and θ 2 (t). In experiments, however, it is often not possible to garner a sufficiently large statistical sample in order to carry out this kind of averaging. For this reason, we propose an analogous order parameter Q S defined on the level of a single-realization of the trajectories θ 1 (t) and θ 2 (t) tracked within a time interval (0, S). We show, via numerical simulations, that both definitions lead to consistent results, i.e., Q S → Q in the limit of an unlimited long observation time S. This result, inter alia, shows that the system under study is ergodic, which cannot be expected a priori, especially if T eff . The outline of the paper is as follows. In Sec. II we present our main results obtained for the model defined by Eqs. (1) and (2). We present here explicit expressions for the probability density function in the steady-state, the steady-state current and an ensemble-averaged order parameter. Further on, we introduce analogous quantities for individual realizations of θ 1 (t) and θ 2 (t). In Sec. III we discuss the behavior of the ensemble-averaged quantities and of their counterparts defined for a single realization of noises, and outline some perspectives for future research. Details of calculations are relegated to the Appendix A. Here, we provide the Fokker-Planck equation associated with the minimal Langevin model in Eqs.
(1) and (2), and present its solution in the limit t → ∞. We also describe our numerical approach, which is based on the discretization of the Langevin equations in Eq. (1).

II. RESULTS
A. Ensemble-averaged properties.
Probability density function in the steady state. Our main analytical result is an exact expression for the joint probability density function (pdf) P (θ 1 , θ 2 ), which is the steady-state solution of the Fokker-Planck equation (see Eq. (A1) in the Appendix A) associated with a system of two coupled Langevin equations (Eq. (1)), and with the noise terms defined by Eq. (2): where T eff is the "mean" effective temperature [37] T eff = This normalization constant can be calculated exactly and is given by Z = I 0 ε/(2T eff ) , where I 0 (x) is a modified Bessel function of the first kind. We note that Eq. 3 here represents a particular case of a more general result derived in Ref. [31].
Since both natural frequencies ν 1 and ν 2 are taken to be equal, the steady-state solution P (θ 1 , θ 2 ) depends on the phases only via the phase difference, and thus becomes independent of ν = ν 1 = ν 2 (see Eq. (3)). At first glance, the latter property appears to be somewhat astonishing, but it can be readily understood once one notices that the time evolution of the phase difference in the Langevin equations (1) becomes independent of ν if both natural frequencies are equal to each other. This means that P (θ 1 , θ 2 ) is defined in the frame of reference rotating with the unique frequency ν. Naturally, the maximum of P (θ 1 , θ 2 ) occurs for θ 1 = θ 2 , regardless of the relation between the temperatures T Individual realizations of (short) trajectories θ1(t) and θ2(t) (recall that θ1(t), θ2(t) ∈ (−1/2, 1/2) are periodic quantities), defined in the reference frame rotating with the frequency ν, for T Hz. The individual trajectories θ1(t) and θ2(t) exhibit a noisy dynamics, but nonetheless evolve alongside each other for rather extended periods of time. This is precisely the stochastic phase locking phenomenon described in Ref. [14] for the case of equal temperatures. It is inferred by following the time evolution of the phase difference ∆t = θ1(t) − θ2(t) in experiments and in numerical simulations. Here, we observe that this kind of stochastic synchronization [14] degrades if T (1) eff = T (2) eff . Indeed, the stochastic phase locking is seemingly strongest in the case of equal temperatures (panels (c) and (f)). It is less pronounced for the combination T of the minimal model with two different temperatures is, that in the non-equilibrium steady-state a nonzero current J occurs. This is a well-known aspect for stochastic dynamics of coupled components, each evolving at its own temperature (see, e.g., Refs. [29][30][31][32][33][34][35][36]). However, in the case at hand this nonzero current has a peculiar form due to the fact that the coupling term in Eq. (1) is a periodic function of the phase difference. The components J 1 and J 2 of this current can be inferred directly from the Fokker-Planck equation (A1) (see Appendix A, Eq. (A2)). They obey The expression on the right-hand-side (rhs) of Eq. (4) contains the trivial term νP (θ 1 , θ 2 ), which is the same for both components as could be expected on general grounds. It appears due to the constant drift term ν = ν 1 = ν 2 on the rhs of the Langevin equations (1). In addition, there is a non-trivial contribution j(θ 1 , θ 2 ), which is a steady-state current in the frame rotating with the unique frequency ν; it reads eff . Rather unexpectedly, j(θ 1 , θ 2 ) appears also to be the same for both components J 1 and J 2 of the current J, due to the form of the pdf in Eq. (3).
One can straighforwardly check that  [14,21]. The curves correspond to the analytical prediction made in Eq. (3). The symbols represent the results of the numerical simulations (see Sec. III), based on a single-realization, time-averaged P S (θ1, θ2) (Eq. (12)) with S = 10 4 sec (such that for a typical value of the frequency ν 47 Hz [14], each flagella makes, on average, 4.7 × 10 5 full beats). Panel (b): Sign of the out-of-equilibrium current j(θ1, θ2) (Eq. (5)), which causes synchronized "drifts" of phases θ1(t) and θ2(t), on the periodic (θ1, θ2)-plane for T eff . The curves correspond to the analytical prediction made in Eq. (7), which defines the ensemble-averaged order parameter Q. The symbols represent the results of the numerical simulations for the time-averaged order parameter Q S (Eq. (14)), based on a single realization of the noises with S = 10 4 sec (i.e., approximately 2.8 hours). Note that for T (2) eff = 0.5 Hz, i.e., in the case of equal effective temperatures, both the ensemble-averaged order parameter Q and its time-averaged counterpart Q S are equal to zero. θ 2 ). On the other hand, j(θ 1 , θ 2 ) is not equal to zero locally (except for θ 1 = θ 2 and |θ 1 −θ 2 | = 1/2, where the current changes sign), and its sign and amplitude depend on the precise values of the phases θ 1 and θ 2 . In Fig. 2(b), for a particular example with T (1) eff > T (2) eff , we present a "phase chart" for the sign (i.e., the direction) of the out-of-equilibrium current j(θ 1 , θ 2 ) in the periodic (θ 1 , θ 2 )-plane. Further on, in Fig. 2(c) we show the current j(θ 1 , θ 2 ) (Eq. (5)) as a function of θ 1 for several fixed values of θ 2 , which also provides insight into its amplitude.
Out-of-equilibrium synchronization. Equations (4) and (5) demonstrate that in a minimal model with two different effective temperatures, in addition to a stochastic phase locking of the coupled phases θ 1 and θ 2 (as observed in Ref. [14]), there is a different synchronization mechanism (based on the out-of-equilibrium current j(θ 1 , θ 2 )) which manifests itself via drifts of the phases. These drifts are correlated in that they have the same sign (i.e., direction) and the same amplitude for both phases. The actual direction of such drifts depends on the sign of ∆T eff , as well as on the relative positions of θ 1 and θ 2 with respect to each other. In order to illustrate this behavior, we suppose ∆T eff > 0 and 0 < θ 2 − θ 1 < 1/2. In this case, according to Eq. (5), both θ 1 and θ 2 experience a drift in the negative direction up to the time at which, due to the thermal noise, the phase difference exceeds the value 1/2 so that θ 2 − θ 1 > 1/2. Then, both θ 1 and θ 2 revert the direction of their drift. Once θ 1 and θ 2 interchange their positions, such that −1/2 < θ 2 −θ 1 < 0, the current j(θ 1 , θ 2 ) changes sign and turns positive, so that both θ 1 and θ 2 drift in the positive direction. Once θ 2 − θ 1 < −1/2, the drift direction again changes sign and becomes negative.
An order parameter in the steady-state. In order to quantify this novel synchronization mechanism, and also in order to render it observable either experimentally or in numerical simulations, one has to introduce a meaningful order parameter. In view of the above discussion, for a minimal model with two effective temperatures the latter should be associated with the steady-state current j(θ 1 , θ 2 ). As in general, there is, however, some liberty in choosing this parameter. Here we define an order parameter by integrating the out-of-equilibrium current j(θ 1 , θ 2 ) over θ 2 across half of the domain in which this variable is defined, and dividing the result by the mean effective temperature. This gives the following dimensionless order parameter (see Eqs. (3) and (5)) This order parameter Q is a function of the phase θ 1 and depends on the coupling parameter ε as well as on the values of the effective temperatures. It vanishes for ε → 0, for T eff → ∞, and for T (1) eff → T (2) eff , the latter limit being characteristic of the transition to the equilibrium setup. The order parameter also vanishes for θ 1 = ±π/4.
It is rewarding to determine the asymptotic behavior of Q in several particular limits. For instance, in the high-temperature limit one has which reduces to if one of the effective temperatures is much higher than the other one.
In the opposite limit T eff (ε cos(2πθ 1 ))/2, which can be reached either via a sufficiently strong coupling ε (and for such values of θ 1 for which cos(2πθ 1 ) is nonzero) or if both temperatures are sufficiently small. In these cases the system is close to the realm of the standard noiseless Kuramoto model and one finds directly from Eq. (7) that, in leading order in the parameter ε/(2T eff ) → ∞, the order parameter Q varies as In these limits T (1,2) eff → 0 and for any θ 1 = 0, Q(θ 1 ) is exponentially small. For θ 1 = 0, the exponential factor in the latter expression equals 1 and hence the order parameter Q varies algebraically as function of the effective temperatures and the coupling parameter ε: i.e., it is a non-analytic function of the coupling parameter and the effective temperatures.

B. Time-averaged properties for individual realizations of noises.
It is not always possible to generate a statistical sample of large enough size, either in experiments or in numerical simulations, which allows one to average over an ensemble of trajectories. To this end, we present alternative definitions for the pdf P (θ 1 , θ 2 ), the current j(θ 1 , θ 1 ), and the order parameter Q, based on their time-averaged counterparts.
The pdf for a single realization of noises. The pdf P (θ 1 , θ 2 ) (Eq. (3)) obeys P (θ 1 , θ 2 ) = lim S→∞ P S (θ 1 , θ 2 ) with where θ 1 (t) and θ 2 (t) are two individual realizations of the trajectories of the phases, corresponding to the solutions of the Langevin equations (Eq. (1)) for a given realization of the noises ζ 1 (t) and ζ 2 (t) in Eq. (2). In Eq. (12), P S (θ 1 , θ 2 ) is the total number of simultaneous occurrences, within the time interval (0, S), of two given realizations of the trajectories θ 1 (t) and θ 2 (t) at the positions θ 1 and θ 2 , respectively, divided by the observation time S. If the system under study is ergodic, as it is the case (see Sec. III), the ensemble average and the time average yield identical results, such that, in the limit S → ∞, P S (θ 1 , θ 2 ) should attain P (θ 1 , θ 2 ). Time-averaged current. We introduce the current j S (θ 1 , θ 2 ) as an average over the observation time S: where θ 1 (t) and θ 2 (t) obey Eq. (1) with ν set to zero. Note that the expressions in the first and the second line in Eq. (13) correspond to the components of the current J and differ with respect to the time derivative of the phases, i.e.,θ 1 (t) orθ 2 (t). As we have shown above, the components of the ensemble-averaged current are exactly equal to each other. We thus expect (and verify via numerical simulations) that the same holds for their introduced time-averaged counterparts.
Time-averaged order parameter Q S . We integrate the expression in the first line on the rhs of Eq. (13) over θ 2 across one half of the domain in which this variable is defined. Dividing the result by the mean effective temperature (see the definition of Q in Eq. (7)), we obtain where θ H (z) is the Heaviside theta-function, which is zero for z < 0 and 1 for z > 0. We expect that, similarly to the time-averaged quantity P S (θ 1 , θ 2 ) (Eq. (12)) and to the time-averaged current in Eq. (13), for large observation times S → ∞, Q S converges to Q given in Eq. (7).

III. DISCUSSION
In Fig. 1 we present appropriately discretized, individual realizations of the trajectories θ 1 (t) and θ 2 (t) which consist of n = 10 6 steps with a discrete time step δt = 10 −6 sec (see the Appendix A for more details) for two distinct values of the coupling parameter (ε = 0.5 Hz for the upper row and ε = 1 Hz for the lower row), for the fixed effective temperature T eff corresponds to the original one considered in Ref. [14]. In contrast, the trajectories in pan-els (a) and (d) of Fig. 1 correspond, within the present choices, to the extreme case of maximal disparity between the temperatures. In (a), θ 2 (t) corresponds to zero temperature (i.e., a perfect rotary motor operating with no noise) and is entrained in random motion by θ 1 (t), which is subject to random noise. Interestingly, the stochastic phase locking described in Ref. [14] is seemingly strongest in the case of equal temperatures (panels (c) and (f)), is less pronounced for the combination T (2) eff = 0 (panels (a) and (d)) for which the periods of synchrony are hardly visible. Therefore, the synchronization observed in Ref. [14] degrades if the effective temperatures become unequal.
The pdf P (θ 1 , θ 2 ) as a function of θ 1 for several values of θ 2 is illustrated in Fig. 2(a) together with the results of numerical simulations for the time-averaged, singletrajectory quantity P S (θ 1 , θ 2 ) in Eq. (12). The very nice agreement between P (θ 1 , θ 2 ) and P S (θ 1 , θ 2 ) shows indirectly that the system is indeed ergodic. Such an agreement is, however, achieved for trajectories which are substantially longer than the ones shown in Fig. 1.
Here we have used the same δt = 10 −6 sec but a larger value n = 10 10 , so that the observation time is S = 10 4 sec.
We use next the trajectories provided in Fig. 1 in order to obtain the introduced time-averaged current j S (θ 1 , θ 2 ) (Eq. (13)) for individual realizations of θ 1 (t) and θ 2 (t). The results (see the Appendix A for more details) are presented in Fig. 2(c) together with the ensemble-average of j(θ 1 , θ 2 ) (Eq. (5)) as obtained from the solution of the Fokker-Planck equation. The agreement between the two results is very satisfactory for n = 10 12 and δt = 10 −6 sec, such that the observation time S = 10 6 sec. For smaller n the data appear more noisy and no conclusive statement on the convergence of j S (θ 1 , θ 2 ) to j(θ 1 , θ 2 ) can be made.
In Fig. 2(d) we show Q obtained from Eq. (7) together with Q S following from Eq. (14). The latter is obtained from the trajectories depicted in Fig. 1 (with δt = 10 −6 and n = 10 10 , such that S = 10 4 sec), as function of θ 1 for ε = 0.5 Hz, T (1) eff = 0.5 Hz, and three values of T (2) eff . We observe full agreement between our theoretical prediction in Eq. (7), which is defined for an ensemble of trajectories, and Q S as introduced in Eq. (14), which is defined for a single realization of noises. This implies that the latter can be conveniently used for a single-trajectory analysis of corresponding experimental and numerical data. Finally, we note that a rather long observation time S = 10 4 sec has been used in Fig.  2(d) in order to demonstrate convergence of the timeaveraged order parameter to the ensemble-averaged one. The observation that the order parameter Q S deviates from zero in out-of-equilibrium conditions can be made already for more moderate values of n, although the data will look more noisy.
In summary, we have presented a generalization of a minimal model introduced in Ref. [14] to the case in which the phases in Eq. (1) are subject to noises with different amplitudes. This can be thought of as a noisy Kuramoto (or Sakaguchi) model of two coupled oscillators with distinct effective temperatures. From a physical point of view, the original model in Ref. [14] has been introduced in order to describe the noisy synchronization of two identical flagella of a biflagellate alga. Our generalized model is expected to be appropriate for the description of a noisy synchronization of two flagella having different lengths. Indeed, the analysis in Ref. [21] has revealed that the noise amplitudes depend on the length of the flagella. Viewed from a different perspective, our study provides an, apparently first, solvable example for the synchronization of coupled oscillators under out-ofequilibrium conditions. Hence, it opens new perspectives for a similar analysis of more complicated models, such as a FitzHugh-Nagumo model (see, e.g., Ref. [26]). Note that in the example studied here the difference between the effective temperatures is not artificially imposed but emerges naturally.
We have shown, both analytically and numerically, that in such a system a very peculiar form of a synchronization of two coupled oscillators takes place. It is mediated by an emerging, current-carrying steady-state. More specifically, we have shown that, on top of the synchronization of the phases as observed in Ref. [14], i.e., a stochastic phase locking, an additional synchronized rotation (drifts) of the phases takes place. In order to quantify the degree of such a synchronization, we have introduced a characteristic order parameter, which vanishes if the effective temperatures become equal to each other. This order parameter has been determined as the average over an ensemble of realizations of the stochastic evolution of phases, as well as on the level of an individual realization. The latter makes the order parameter suitable for experimental and numerical analyses, for which a sufficiently large statistical sample cannot be formed. Via numerical simulations we have shown that both definitions become equivalent in the limit of sufficiently long observation times, which also demonstrates the ergodicity of the system under study. Finally, we remark that we expect a much richer behavior for the relevant situation in which, in addition to unequal effective temperatures, the natural frequencies are also different. This is a challenging subject for future research. Analytical approach. We provide the Fokker-Planck equation for the joint pdf P (θ 1 , θ 2 ), associated with the system of two coupled Langevin equations (Eq. (1)) with the noise terms defined by Eq. (2). The associated Fokker-Planck equation is derived by standard means [16] and readṡ where J = (J 1 , J 2 ) is the probability current. The steadystate solution of Eq. (A1) can be determined analytically and is given by Eq. (3) in Sec. II.
From Eq. (A1) we infer the following expressions for the components J 1 and J 2 of the out-of-equilibrium current J: Inserting the explicit expression of the pdf in Eq. (3) into Eq. (A2), and performing differentiations, we obtain Eq. (4).
Numerical approach. Here, we provide a brief description of our numerical algorithm. To this end we rewrite the Langevin equations (Eq. (1)) in the frame of reference rotating with the frequency ν, i.e., we change variables according to θ 1,2 =θ 1,2 (t) + νt, and then we discretize the time variable t = nδt, where n is an integer; δt is the time-interval between the consecutive steps. Without loosing generality, in what follows we set ν = 0 and, in order to avoid a clumsy notation, we drop the tilde mark. Lastly, we recall the standard scaling properties of Gaussian delta-correlated noises ζ 1,2 (t) in order to cast the noise terms into a different (but equivalent) form, in which the effective temperatures appear explic-itly as amplitudes of the noises [16]. This turns Eq. (1) into recurrence relations of the form where ∆ t = θ 1 (t) − θ 2 (t) is an instantaneous phase difference. The above recursions (with dimensionless prefactors of the sine and of the noise terms) allow us to define the values of θ 1 (t + δt) and θ 2 (t + δt) through the values of ∆ t and the values of the noise terms η 1,2 (t) at the previous moment of time. This permits us to sequentially generate individual realizations of the phases θ 1 (t) and θ 2 (t) of arbitrary duration t. The noises η 1 (t) and η 2 (t) in Eq. (A4) are dimensionless Gaussian random variables, uncorrelated for distinct values of n, with zero mean and variances σ 2 1,2 ≡ 1, such that the probability density function is given explicitly by P (η 1,2 ) = exp(−η 2 1,2 /2)/ √ 2π. This choice ensures that for ε = 0 the phases θ 1 (t) and θ 2 (t) undergo standard diffusive motion on the unit circle with diffusion coefficients T (1) eff and T (2) eff , respectively. Adopting δt = 10 −6 sec (which is sufficiently short such that for a typical beating frequency of ν 47 Hz (see Ref. [14]) each flagella makes a full beat within roughly 2.13× 10 4 intervals δt), we generate two individual realizations of noises, thereby building up individual realizations of the trajectories θ 1 (t) and θ 2 (t) which consist of n = 10 6 steps. As a consequence, within the full observation time (here S = 1 sec) each flagella makes, on average, 47 full beats. These trajectories are depicted in Fig. 1.
In order to determine numerically the introduced timeaveraged current j S (θ 1 , θ 2 ) (Eq. (13)) for individual realizations of θ 1 (t) and θ 2 (t), we first appropriately discretize the expression on the rhs of Eq. (13) and replace the time-derivativeθ 1 (t) (orθ 2 (t)) by the finite difference given by Eq. (A4). Then, we use the trajectories provided in Fig. 1. The results of such a procedure are presented in Fig. 2(c) together with the ensemble-average of j(θ 1 , θ 2 ) (Eq. (5)) as obtained from the solution of the Fokker-Planck equation. The order parameter Q S , as introduced in Eq. (14), is determined numerically in a similar way.