Current-mediated synchronization of a pair of beating non-identical flagella
V. Dotsenko A. Maciolek G. Oshanin O. A. Vasilyev S. Dietrich















Sorbonne Université, CNRS, Laboratoire de Physique Théorique de la Matière Condensée, LPTMC (UMR CNRS 7600),
75252 Paris Cedex 05, France
L.D. Landau Institute for Theoretical Physics,
           119334 Moscow, Russia

Max-Planck-Institut für Intelligente Systeme,
  Heisenbergstraße 3, D-70569 Stuttgart, Germany
IV. Institut für Theoretische und Angewandte Physik,
  Universität Stuttgart, D-70569 Stuttgart, Germany
Institute of Physical Chemistry, Polish Academy of Sciences,
  Department III, Kasprzaka 44/52, PL-01-224 Warsaw, Poland
  
Sorbonne Université, CNRS, Laboratoire de Physique Théorique de la Matière Condensée, LPTMC (UMR CNRS 7600),
75252 Paris Cedex 05, France
oshanin@lptmc.jussieu.fr

Max-Planck-Institut für Intelligente Systeme,
  Heisenbergstraße 3, D-70569 Stuttgart, Germany
IV. Institut für Theoretische und Angewandte Physik,
  Universität Stuttgart, D-70569 Stuttgart, Germany

Max-Planck-Institut für Intelligente Systeme,
  Heisenbergstraße 3, D-70569 Stuttgart, Germany
IV. Institut für Theoretische und Angewandte Physik,
  Universität Stuttgart, D-70569 Stuttgart, Germany




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_ S - for an ensemble of trajectories and
for a single realization of noises  of duration S, respectively.
Via numerical simulations we show that both approaches become identical for long observation times S.  This 
reveals an ergodic behavior and
implies
that a single-realization order parameter Q_ S is suitable for experimental analysis for which ensemble averaging is not always possible.



      87.16.Qp,   
      05.45.Xt,  
      47.63.-b, 
      87.18.Tt 
      

Current-mediated synchronization  of a pair of beating non-identical flagella
    S. Dietrich
    December 30, 2023
=============================================================================








§ INTRODUCTION


There is  experimental evidence that two beating flagella, extending from one end of the biflagellate
alga Chlamydomonas reinhardtii,
synchronize their dynamics.
Analyzing the oscillatory intensity signals
x_1,2(t) = Γ_1,2(t) sin(2 πθ_1,2(t))
(where Γ_1,2(t)  and θ_1,2(t)
are the amplitudes and  the instantaneous phases of the periodic motion
of  flagella 1 and 2, respectively), which are 
obtained by local sampling of the video light intensity near the two flagella, Polin et al. <cit.> observed that the
phase difference _t = θ_1(t) - θ_2(t) contains periods of synchrony (i.e., the so-called phase locking behavior with Δ_t ≈ const <cit.>),
interrupted by sudden drifts of either sign. 
Referring to earlier ideas, that the hydrodynamic interactions between
eukaryotic flagella or cilia may underlie
their synchronization <cit.>, Goldstein et al. <cit.>
proposed a  phenomenological, minimal, stochastic model in which the motion of a
flagella pair is  described by two noisy phase oscillators
which move on circular trajectories and are
coupled via an antisymmetric function of  the phase difference _t.
In terms of the notations used in Ref. <cit.>, the equations of motion read

    θ̇_1(t)    = ν_1 + πεsin[2 π(θ_2(t)  - θ_1(t))] + ζ_1(t)  , 
    θ̇_2(t)    = ν_2 + πεsin[2 π(θ_1(t)  - θ_2(t))] + ζ_2(t)  .

In Eq. (<ref>), 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, respectively, 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 <cit.>:
ζ_i(t) ζ_j(t') = 2 T_ 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  <cit.>. 
Indeed, these active fluctuations are several orders of magnitude 
larger than the thermal noise <cit.>, so that the flagella 
are not in thermal equilibrium with its bath.


We note that the sine terms in Eq. (<ref>), which  link
the  time evolutions of  θ_1(t) and θ_2(t), describe the actual coupling
due to hydrodynamic interactions between the two flagella <cit.> 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. <cit.>, 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. <cit.>. We also note that in more complex situations of unicellular algal species bearing multiple flagella, both the effective coupling introduced in Ref. <cit.> and the results in Refs. <cit.>
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  <cit.>.



Solving Eq. (<ref>) numerically for given realizations of the noise terms, Goldstein et al. <cit.> have found consistency between  Δ_t evolving according to Eq. (<ref>) and
the experimentally observed behavior of the flagella  of the biflagellate
alga <cit.>, 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. <cit.>. Importantly, the results of Ref. <cit.>
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-to-realization fluctuations.

A more sophisticated experimental analysis has been performed in
Ref. <cit.>, 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 <cit.>. 
The flagella of length l ≃ 10.82 μm have been first clipped by a micropipette <cit.>
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 <cit.> 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.  <cit.>). 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. <cit.>. 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 <cit.>.


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. <cit.>. 
We note that in this case the regeneration scenario is more complicated, as compared to the case when both flagella are removed (see Ref. <cit.>). 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 <cit.>. It  can become even longer if certain chemicals (e.g., colchicine) are added after a deflagellation, which inhibit the regeneration process <cit.>. 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 <cit.> associated with the Langevin equations (<ref>)
 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^(1)_ eff = T^(2)_ eff.


 Viewed from a different perspective,  which is perhaps equally important
 due to certain other applications (see Ref. <cit.> for a discussion)
 we note that  the minimal  model in Eq. (<ref>)
 with a unique effective temperature T_ eff represents the so-called Sakaguchi model <cit.> which is
a noisy version of the celebrated Kuramoto model of coupled oscillators
(see, e.g.,  Refs. <cit.>).
Its generalization to the case of
two different temperatures <cit.> 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. <cit.> 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. (<ref>) with the covariance
 functions of the noise terms of the form

    ζ_i(t) ζ_j(t) = 2  δ_i,j  T^(i)_ eff δ(t-t')  ,     i,j = 1,2  ,

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)  <cit.>, 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 reference rotating with  frequency ν (note that j(θ_1,θ_2) ≡ 0 
  when  T^(1)_ eff =  T^(2)_ 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. <cit.>  degrades
 for unequal effective temperatures (see Fig. <ref>) 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^(1)_ eff≠ T^(2)_ eff.


The outline of the paper is as follows. In Sec. <ref> we present our main results obtained for the model defined by Eqs. (<ref>) and (<ref>). 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. <ref> 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 <ref>. Here, we
provide the Fokker-Planck equation associated with the minimal Langevin model in 
Eqs. (<ref>) and (<ref>), 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. (<ref>).





§ RESULTS





 §.§ 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. (<ref>) in the Appendix <ref>) associated with a system of two coupled Langevin equations (Eq. (<ref>)), and
with the noise terms defined by Eq.  (<ref>):

    P(θ_1,θ_2) = 1/Zexp(ε/2 T_ effcos(2 π(θ_1 - θ_2)))  ,

where T_ eff is the "mean" effective temperature <cit.>
T_ eff = (T^(1)_ eff + T^(2)_ eff)/2, while
Z insures that P(θ_1,θ_2) is properly normalized, i.e.,
∫_-1/2^1/2∫^1/2_-1/2 dθ_1 dθ_2 P(θ_1,θ_2) = 1.
This normalization constant can be calculated exactly and is given by
Z  =  I_0(ε/(2 T_ eff) ),
where I_0(x) is a modified Bessel function of the first kind.  We note that  Eq. <ref> here
represents a particular case of a more general result derived in Ref. <cit.>.


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. (<ref>)). 
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 (<ref>) 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^(1)_ eff and
 T^(2)_ eff. For ε/(2 T_ eff) →∞, the pdf turns into 
 a delta function of the difference θ_1 - θ_2.  Figure <ref>(a) provides 
 the pdf P(θ_1,θ_2) 
  (Eq. (<ref>)) as a function of θ_1 for several fixed values of θ_2.




Out-of-equilibrium current. A remarkable feature 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. <cit.>). However, in the case at hand this nonzero current has a peculiar form due to the fact that the coupling term in Eq. (<ref>) 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 (<ref>) (see Appendix <ref>, Eq. (<ref>)).
They obey

    J_1    =  J_2 = j(θ_1,θ_2) + ν P(θ_1,θ_2)   .


The expression on the right-hand-side (rhs) of Eq. (<ref>) 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 (<ref>). 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

    j(θ_1,θ_2) 
        = - πεΔ T_ eff/2 T_ effsin(2 π (θ_2 - θ_1)) P(θ_1,θ_2)   ,

with Δ T_ eff = T^(1)_ eff - T^(2)_ 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. (<ref>).

The mean out-of-equilibrium current,

    ⟨ j(θ_1,θ_2) ⟩≡∫^1/2_-1/2∫^1/2_-1/2 dθ_1 dθ_2   j(θ_1,θ_2)  = 0  ,

vanishes such that, due to ⟨ P(θ_1,θ_2) ⟩ = 1, 
⟨ J_1⟩ = ⟨ J_2 ⟩ = ν. 
One can straighforwardly check that  ∫^1/2_-1/2 dθ_1 j(θ_1,θ_2)   =   0   =  ∫^1/2_-1/2 dθ_2 j(θ_1,θ_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. <ref>(b), for a particular example with T_ eff^(1) >  T_ eff^(2), 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. <ref>(c) we show the current j(θ_1,θ_2) (Eq. (<ref>)) as a function of θ_1 for several fixed values of θ_2, which also provides insight into its amplitude.



 Out-of-equilibrium synchronization. Equations (<ref>) and (<ref>)  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. <cit.>), 
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. (<ref>),
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. (<ref>) and (<ref>))

    Q(θ_1)    = 1/T_ eff∫^1/2_0 dθ_2    j(θ_1,θ_2) 
       =  Δ T_ eff/T_ eff sinh((ε/(2 T_ eff) cos(2 πθ_1) )/ I_0(ε/(2 T_ eff) ) .

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_ eff^(1)→ T_ eff^(2),  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

    Q(θ_1)  ≃ Δ T_ eff/2 T^2_ effεcos(2 πθ_1)  ,    T_ eff≫ε/2,

which reduces to

    Q(θ_1) ≃2   sign(T_ eff^(1) - T_ eff^(2))/max(T_ eff^(1),T_ eff^(2)) εcos(2 πθ_1)  ,

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. (<ref>) 
that, in  leading order in the parameter ε/(2 T_ eff) →∞,  the order parameter Q  varies as

    Q(θ_1)  ≃ Δ T_ eff/2 T^3/2_ eff (πε)^1/2exp(- ε/T_ effsin^2(πθ_1)).

In these limits T_ eff^(1,2)→ 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 ε: 

    Q(θ_1=0) ≃Δ T_ eff (πε)^1/2/2 T^3/2_ eff ,

i.e.,  it is a non-analytic function of the coupling parameter and the effective temperatures.




 §.§ 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. (<ref>))  obeys P(θ_1,θ_2) = lim_ S→∞ P_ S(θ_1,θ_2) with

    P_ S(θ_1,θ_2) = 1/ S∫^ S_0 dt  δ(θ_1(t) - θ_1) δ(θ_2(t) - θ_2)  ,

 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. (<ref>)) for
a given realization of the noises ζ_1(t) and ζ_2(t) in Eq. (<ref>). 
In Eq. (<ref>), 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. <ref>), 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:

    j_ S(θ_1,θ_2)    = -   1/ S∫^ S_0 dt  θ̇_1(t) δ(θ_1(t) - θ_1) δ(θ_2(t) - θ_2) 
        = -  1/ S∫^ S_0 dt  θ̇_2(t) δ(θ_1(t) - θ_1) δ(θ_2(t) - θ_2)  ,

where θ_1(t) and θ_2(t) obey Eq. (<ref>) with ν set  to zero. Note that the expressions in the first and the second line in Eq. (<ref>) 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. (<ref>) 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. (<ref>)), we obtain

    Q_ S   = -   1/ ST_ eff∫^ S_0 dt  θ̇_1(t) δ(θ_1(t) - θ_1)  
       ×∫_0^1/2 dθ_2  δ(θ_2(t) - θ_2) 
       = -   1/ ST_ eff∫^ S_0 dt  θ̇_1(t) δ(θ_1(t) - θ_1) θ_H(θ_2(t))  ,

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. (<ref>))  and to the time-averaged current in Eq. (<ref>), 
for large observation times S→∞,
Q_ S converges to Q given 
in Eq. (<ref>).  




§ DISCUSSION



In Fig. <ref> 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 <ref> 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^(1) = 0.5   Hz of flagellum 1, and
for  three  temperatures of  flagellum
2 (T_ eff^(2) = 0, 0.1   Hz, and 0.5   Hz).
The case  T_ eff^(2) = 0.5   Hz, i.e., T_ eff^(1) = T_ eff^(2)
 corresponds to   the original one
considered in Ref. <cit.>. In contrast,
the trajectories in panels (a) and (d) of Fig. <ref> 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. <cit.> is seemingly
strongest in the case of equal temperatures (panels (c) and (f)),
is less pronounced for the combination  T_ eff^(1) = 0.5   Hz and T_ eff^(2) = 0.1   Hz
(panels (b) and (e)) and it is weakest for T_ eff^(1) = 0.5   Hz and T_ eff^(2) = 0 (panels (a) and (d))
for which the periods of synchrony are hardly visible. Therefore, the synchronization observed
in Ref. <cit.> 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. <ref>(a)
together with the results of numerical simulations for the time-averaged, single-trajectory quantity  P_ S(θ_1,θ_2) in Eq. (<ref>).
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. <ref>. 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. <ref> in order to obtain the introduced  time-averaged current j_ S(θ_1,θ_2) (Eq.  (<ref>)) for
  individual realizations of θ_1(t) and θ_2(t). The results (see the Appendix <ref> for more details) are presented 
in Fig. <ref>(c) together with
the ensemble-average of  j(θ_1,θ_2) (Eq. (<ref>)) 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. <ref>(d)
we show Q obtained from Eq. (<ref>) 
together with Q_ S following from Eq. (<ref>). 
The latter is obtained 
from the trajectories depicted in Fig. <ref>  (with δ t = 10^-6 and n=10^10, such that S=10^4 sec), as function of 
θ_1 for ε = 0.5   Hz, T_ eff^(1) = 0.5   Hz,   
and three  values of T_ eff^(2). We observe full agreement between 
our theoretical prediction in Eq. (<ref>), which is defined for an ensemble of trajectories, 
and Q_ S as introduced in Eq. (<ref>), 
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. <ref>(d) 
 in order
to demonstrate convergence of the time-averaged 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. <cit.> to  the case in which the phases
in Eq. (<ref>) 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. <cit.> 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. <cit.>
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-of-equilibrium conditions.  
Hence, it opens new perspectives
for a similar analysis of more complicated models, such as a FitzHugh-Nagumo model (see, e.g., Ref. <cit.>).
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. <cit.>,  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. 





§ DETAILS OF CALCULATIONS.




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. (<ref>)) 
with the noise terms defined by Eq.  (<ref>).
The associated Fokker-Planck equation is derived by standard means <cit.> and reads

    Ṗ(θ_1,θ_2) =  -  div J =
     - ∂/∂θ_1  J_1  -  ∂/∂θ_2  J_2
       = (∂/∂θ_1[T^(1)_ eff∂/∂θ_1 + (πεsin(2 π (θ_1 - θ_2))  - ν)] 
        + ∂/∂θ_2[T^(2)_ eff∂/∂θ_2 - (πεsin(2 π (θ_1 - θ_2))  + ν)] ) P(θ_1,θ_2)  ,

where J = (J_1, J_2) is the probability current.
The steady-state solution of Eq. (<ref>) can be determined analytically and is given by Eq. (<ref>) in Sec. <ref>.

From Eq. (<ref>) we infer the following expressions for the components J_1 and J_2 of the out-of-equilibrium current J:

    J_1    =  -(T^(1)_ eff∂/∂θ_1 + (πεsin(2 π(θ_1 - θ_2)) - ν) )P(θ_1,θ_2) 
    
    J_2    =  - (T^(2)_ eff∂/∂θ_2 - (πεsin(2 π(θ_1 - θ_2)) + ν) ) P(θ_1,θ_2) .

Inserting the explicit expression of the pdf in Eq. (<ref>) into Eq. (<ref>), and performing differentiations, we obtain Eq. (<ref>).

Rewriting the components J_1 and J_2 of the out-of-equilibrium current J in that frame of reference 
which rotates 
with the frequency ν, for j(θ_1,θ_2) (see the definition in Eq. (<ref>)) we find 

    j(θ_1,θ_2)    =  -(T^(1)_ eff∂/∂θ_1 + πεsin(2 π(θ_1 - θ_2)) ) P(θ_1,θ_2) 
       = - (T^(2)_ eff∂/∂θ_2 - πεsin(2 π(θ_1 - θ_2))) P(θ_1,θ_2)   .
       = - πεΔ T_ eff/2 T_ effsin(2 π (θ_2 - θ_1)) P(θ_1,θ_2)   .

 The expression in the last line in Eq. (<ref>) corresponds to our Eq. (<ref>).


Numerical approach. Here, we provide a brief description of our numerical algorithm.
To this end we rewrite the Langevin equations (Eq. (<ref>)) 
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 explicitly as amplitudes of the noises <cit.>.
This turns  Eq. (<ref>) into recurrence relations of the form

    θ_1(t+ δ t)    = θ_1(t) - δ t (πεsin[2 πΔ_t] +  √(2  T_ eff^(1)/δ t)η_1(t))  , 
    θ_2(t + δ t)    = θ_2(t) + δ t (πεsin[2 πΔ_t] +   √(2  T_ eff^(2)/δ t)η_2(t))  ,

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. (<ref>) 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(-η_1,2^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. <cit.>) 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. <ref>. 

In order to determine numerically the introduced  time-averaged current j_ S(θ_1,θ_2) (Eq.  (<ref>)) for
  individual realizations of θ_1(t) and θ_2(t), we  first appropriately discretize the expression on the rhs of Eq. (<ref>) 
and replace the time-derivative θ̇_1(t) (or θ̇_2(t)) by the finite difference given by  Eq. (<ref>). 
Then,  we use
the trajectories provided in Fig. <ref>. The results of such a procedure are presented 
in Fig. <ref>(c) together with
the ensemble-average of  j(θ_1,θ_2) (Eq. (<ref>)) as obtained from the solution of the Fokker-Planck equation. The order parameter
 Q_ S, as introduced in Eq. (<ref>), is determined numerically in a similar way.





§ ACKNOWLEDGMENTS


The authors acknowledge valuable discussions with M. Dykman, F. Ritort, M. Rosenblum, and S. Ruffo.
The work by A.M. has been  partially  supported by the Polish National Science Center
(Harmonia Grant No. 2015/18/M/ST3/00403).




99

1 M. Polin, I. Tuval, K. Drescher, J. P. Gollub, and R. E. Goldstein, Chlamydomonas
 swims with two "gears" in a eukaryotic version of run-and-tumble locomotion,
 Science 325, 487 (2009).
 
 11a S. H. Strogatz, From Kuramoto to Crawford: exploring the onset of
synchronization in populations of coupled oscillators, Physica D 143, 1 (2000).

m2 A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization:
A Universal Concept in Nonlinear Sciences
(Cambridge University Press, Cambridge, 2003).

11   J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler, The Kuramoto model: A simple paradigm for synchronization phenomena,
Rev. Mod. Phys. 77, 137 (2005).

11b S. Gupta, A. Campa, and S. Ruffo, Statistical Physics of Synchronization (Springer, Berlin, 2018).

2 G. I. Taylor,  Analysis of the swimming of microscopic organisms, Proc. R. Soc. A 209, 447 (1951).

3 K. E. Machin, The control and synchronization of flagellar movement, Proc. R. Soc. B 158, 88 (1963).

4 M. Cosentino Lagomarsino, P. Jona, and B. Bassetti, Metachronal waves for deterministic switching 
two-state oscillators with hydrodynamic interaction,
Phys. Rev. E 68, 021908 (2003).


5 P. Lenz and A. Ryskin, Collective effects in ciliar arrays,
Phys. Biol.  3, 285 (2006).


6 S. Gueron and K. Levit-Gurevich,  Energetic considerations of ciliary beating and the advantage of metachronal coordination,
Proc. Natl. Acad. Sci. USA,  96, 12240 (1999).


7 A. Vilfan and F. Jülicher,
Hydrodynamic low patterns and synchronization of beating cilia,
Phys. Rev. Lett. 96, 058102 (2006).

8 B. Guirao and J. F. Joanny, Spontaneous creation of macroscopic flow and metachronal waves in an array of cilia,
Biophys. J.  92, 1900 (2007).


9 T. Niedermayer, B. Eckhardt, and P. Lenz, Synchronization, phase locking, and metachronal wave formation in ciliary chains, Chaos 18, 037128 (2008).


10 R. E. Goldstein, M. Polin, and I. Tuval,
Noise and  synchronization in pairs of beating eukaryotic flagella,
Phys. Rev. Lett. 103, 168103 (2009).

remark Note that the definition of the covariance of the noise terms used
in Ref. <cit.> does not contain the usual prefactor 2 of the temperature  (see, e.g., Ref. <cit.>). Thus,
in our settings T_ eff is  half the analogous  quantity in Ref. <cit.>.
The values of T_ eff, which we present in the main text, take this property into account.

kalm C. W. Gardiner, Handbook of Stochastic Methods (Springer, Berlin, 1983).

10a H. C. Berg, The rotary motor of bacterial flagella, Annu. Rev. Biochem. 72, 19 (2003).

10b R. Ma, G. S. Klindt, I. Riedel-Kruse, F. Jülicher and B. M. Friedrich, 
Active phase and amplitude fluctuations of flaggelar beating, Phys. Rev. Lett. 113, 048101 (2014).

9a K. Y. Wan and R. E. Goldstein, Coordinated beating of algal flagella is mediated by basal coupling, 
Proc. Natl. Acad. Sci. USA 113 (20), E2784 (2016).

21 Y. Liu, R. Claydon, M. Polin, and D. R. Brumley,
 Transitions in synchronization states of model cilia through basal-connection coupling, J. R. Soc. Interface 15, 20180450 (2018).


c R. E. Goldstein, M. Polin, and I. Tuval, Emergence of synchronized beating during the regrowth of eukaryotic flagella, Phys. Rev. Lett.
107, 148203 (2011).

non-equal1 J. L. Rosenbaum, J. E. Moulder and D. L. Ringo, Flagellar elongation and shortening in Chlamydomonas, J. Cell. Biol. 41, 600 (1969).

22 We note that the amplitude
of the noise of a single flagellum as a function of its length l might be different from that
inferred from the experimental data in Ref.  <cit.> for a pair of flagella.

non-equal2 N. F. Wilson, 
J. K. Iyer, J. A. Buchheim, and W. Meek, Regulation of flagellar length in Chlamydomonas,
Seminars in Cell & Developmental Biology 19, 494  (2008). 

m1  A. Malakhov, Influence of the intrinsic noise of oscillators on their mutual synchronization, Investiya VUZ Radiofizika 8,  1160 (1965).

m3 R. M. Amro, B. Lindner, and A. B. Neiman, Phase diffusion in unequally noisy coupled oscillators, Phys. Rev. Lett. 115, 034101 (2015).

12a H. Sakaguchi,  Cooperative phenomena in coupled oscillator systems under external fields, Prog. Theor. Phys. 79, 39  (1988).

25 The dynamics of two coupled Kuramoto oscillators
exposed to two unequal noises
has been already studied (see Refs. <cit.>).
However, for this situation the only (albeit quite complicated)
issue, which was considered, concerns the form of the
diffusion coefficient of the phase difference Δ_t.

12 R. Filliger and P. Reimann, Brownian  gyrator: a minimal heat engine on the nanoscale,
Phys. Rev. Lett. 99, 230602 (2007).

13 V. Dotsenko, A. Maciolek, O. Vasilyev, and G. Oshanin, Two-temperature Langevin dynamics in a parabolic potential,
Phys. Rev. E 87, 062130 (2013).

14 A. Y. Grosberg and J.-F. Joanny,  Nonequilibrium statistical mechanics of mixtures of particles in contact with different thermostats,
Phys. Rev. E 92, 032118 (2015).

15 Y. Murashita and M. Esposito, Overdamped stochastic thermodynamics with multiple reservoirs, Phys. Rev. E 94, 062148 (2016).

16 S. Ciliberto,
Experiments in stochastic thermodynamics: short history and perspectives,
Phys. Rev. X 7, 021051 (2017).

17 H. C. Fogedby and A. Imparato, A minimal model of an autonomous thermal motor, EPL 119,  50007 (2017).

18 V. Mancois, B. Marcos, P. Viot, and D. Wilkowski, Two-temperature Brownian dynamics of a particle in a confining potential, Phys. Rev. E 97, 052121 (2018).

19 S. Cerasoli, V. Dotsenko, G. Oshanin, and L. Rondoni, Asymmetry relations and effective temperatures for biased Brownian
 gyrators, Phys. Rev. E 98, 042149 (2018).

35 
The effective temperatures enter the steady-state pdf in Eq. (<ref>) only 
in an additive way, such that there is  no singular behavior in the case that
 they are unequal;  here the only pertinent quantity
 is the mean effective temperature T_ eff. We note that this seems to be a general rule (see Eq. (3) in Ref. <cit.>)
 for any form of the coupling term between the phases θ_1(t) and θ_2(t) in Eq. (<ref>),  
 which depends on the phase difference only.
 As a consequence,  the  standard complex
 order parameter 
 r exp(i Ψ) ≡ (exp(2π i θ_1) + exp(2π  i θ_2))/2, where r defines the phase coherence and Ψ the average phase
  <cit.>, averaged over the pdf given in Eq. (<ref>), is exactly equal to zero. 




