Jérôme Dorignac

Fred Geniet

Estelle Pitard

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

TASEP Exit Times

Jérôme Dorignac, Fred Geniet, and Estelle Pitard Laboratoire Charles Coulomb (L2C), Université de Montpellier, CNRS, 34095 Montpellier, France (Dated: October 13, 2023)

We address the question of the time needed by N particles, initially located on the first sites of a finite 1D lattice of size L, to exit that lattice when they move according to a TASEP transport model. Using analytical calculations and numerical simulations, we show that when N ≪ L, the mean exit time of the particles is asymptotically given by TN (L) ∼ L + βN √ L for large lattices. Building upon exact results obtained for 2 particles, we devise an approximate continuous space and time description of the random motion of the particles that provides an analytical recursive relation for the coefficients βN . The results are shown to be in very good agreement with numerical results. This approach sheds some light on the exit dynamics of N particles in the regime where N is finite while the lattice size L → ∞. This complements previous asymptotic results obtained by Johansson in [1] in the limit where both N and L tend to infinity while keeping the particle density N/L finite.

I. INTRODUCTION

The TASEP model (Totally Asymmetric Simple Exclusion Process) is a unidirectional model of transport of particles with exclusion on a one dimensional lattice [2]. It has various interesting applications in traffic on lanes, waiting times lists, directed transport of particles through channels and more [START_REF] Schadschneider | Stochastic transport in complex systems: from molecules to vehicles[END_REF][START_REF] Arita | [END_REF][START_REF] Arita | the Conference on Pedestrian and Evacuation Dynamics 2014 (PED 2014)[END_REF]. It can also be mapped on models of interface growth [1,[START_REF] Sasamoto | [END_REF], providing alternate interpretations of its results. Originally introduced in the context of the kinetics of biopolymerization, it has also been a paradigmatic model in the field of biological transport since [7,8].

Most theoretical investigations of the TASEP model have been dedicated to obtaining results at stationarity when the flux of particles entering and exiting the lattice has reached a stationary value. In that respect, particle density and current properties have been thoroughly studied [9][10][11][12][13]. But some results have also been obtained in non-stationary regimes, especially in infinite lattices. For instance, the exact Green functions of the continuous time TASEP model on Z have been obtained by Schütz in [14]. Related quantities have subsequently been used to determine some asymptotic features of the time evolution of the particle density when starting from a stepinitial condition where particles initially populate the left half of the lattice only [15]. Much in the same vein, the statistical features of the motion of certain (tagged) particles along the lattice have been elucidated as well [START_REF] Sasamoto | [END_REF]16]. The question we address here pertains to that class of non-stationary problems: how a set of particles, transported according to the TASEP rules, evacuate a finite lattice, especially when they start from a "steplike" configuration where all of them are located on the leftmost sites of that lattice? To answer that question, we shall study the distribution of their exit time and, more specifically, their mean exit time.

Studies on exit times (also called evacuation times or escape times) in single-file systems, that is in 1D sys-tems where particles cannot pass each other, generally involve bidirectional motion like in single-file diffusion (SFD) problems (see for instance [17,18] and references therein). In this context, exit time distributions may be analyzed via the first passage time density of a "tracer" (or tagged) particle moving within a crowd of like particles (see e.g. [19]). Analysis of these SFD problems shows that the tracer position x(t) has a subdiffusive behaviour leading to a mean-squared displacement that scales as ⟨(x(t)-x 0 ) 2 ⟩ ∝ t 2H at long times where H is the Hurst exponent [20,21]. This behaviour, due to crowding effects generated by 1D confinement at a given density of particles, contrasts with the ⟨(x(t)-x 0 ) 2 ⟩ ∝ t scaling typical of a diffusive behaviour for which H = 1/2. In the particular case of the symmetric exclusion process (SEP) for instance -a 1D hard-core lattice gas problem equivalent to the TASEP model but where particles may equally jump to the right or to the left provided the corresponding site is empty (see e.g. [22]) -the Hurst exponent is H = 1/4 and the SEP problem has been shown to be equivalent to a fractional Brownian motion (fBm) that depends on the particle density [23].

The exit of particles following TASEP transport rules from a finite size lattice share some similarities with SFD systems. In particular, the motion of a given particle is hindered by others (exclusion) and therefore cannot perform a simple independent random walk. But there are two main differences between SFD problems and the question we address in this paper. First, the particle density does not remain constant over time because particles progressively leave the system they start from and free the motion of those that remain within the lattice. In that respect, the escape of colloidal particles from microfluidic channels studied in ref. [24] is a problem closer to ours. Second, the motion of particles is unidirectional (particles only move to the right) that is, transport is totally biased. This situation is similar to emergency evacuation in trains or aircrafts where individuals have to quickly walk down a narrow seat aisle [25,26].

The evacuation of particles according to TASEP rules might therefore provide some insight on emergency evacuation although pedestrian dynamics has a quite complex stochastic structure [27,28].

In this paper, we focus on a special setting of the TASEP model: particles start from a step-like initial state and no particle is injected at site 1. Moreover, all particles exit the lattice as they reach site L + 1 (absorbing condition after the Lth site). At time t = 0, the N particles are located on sites 1, ..N with N ≤ L, as displayed in Figure 1. We are interested in the emptying time of this model which is equal to the exit time of the leftmost particle of the lattice. In particular, we shall use analytical calculations and numerical simulations to determine the mean exit time (MET).

After introducing the model and quantities of interest in Section II, we present some exact results for 1 and 2 particles in Section III and Section IV. We then use a continuous space and time description of the relative motion of the particles with respect to the leading one to calculate the exit time in the large L limit in Section V. This approach provides a simple physical approximate solution of the problem, yielding a recursive expression for T N (L) in the large L limit. These results are then compared to Gillespie simulations in Section VI. In section VII, we discuss our asymptotic results and compare them to those of Johansson [1] obtained in the finite density regime.

II. TRANSPORT MODEL AND ITS EXIT TIME DISTRIBUTION

A. The TASEP Model

The Totally Asymmetric Simple Exclusion Process (TASEP) is a paradigmatic dynamical model for the unidirectional transport of particles on a lattice that takes into account exclusion. It is generally defined by the following rules: a particle may be loaded on the lattice with probability rate α provided the first site is empty. It then proceeds forward with a hopping rate p provided the neighbouring site (on the right) is empty and leaves the last lattice site with a probability exit rate β. In what follows, we shall study the exit time distribution of N particles initially located on the first (leftmost) N sites of a lattice containing L ≥ N sites, see figure 1 for a pictorial view. We shall moreover assume that the hopping rates on the lattice are homogeneous, p = β and using 1/p as unit of time, we shall simply set p = β = 1. Finally, the incoming rate α is set to zero in such a way that no particle enters the lattice from t = 0 onward.

The TASEP model is a Markov process governed by the master equation,

d|P (t)⟩ dt = M |P (t)⟩ (1)
where the probability vector may be written as

|P (t)⟩ = σ P σ (t)|σ⟩ . (2)
Here, the configuration vector is

|σ⟩ = |σ 1 ⟩ ⊗ • • • ⊗ |σ L ⟩ with column vectors |σ i ⟩ = (1 -σ i , σ i ) T where σ i = 1
when site i is occupied by a particle and σ i = 0 otherwise. The sum runs over all 2 L particle configurations. Within our settings, the Markov matrix M reads [29]

M = L-1 i=1 1I i-1 ⊗ m ⊗ 1I L-1-i + 1I L-1 ⊗ b . ( 3
)
where 1I is the 2 × 2 identity matrix and where

m =    0 0 0 0 0 0 1 0 0 0 -1 0 0 0 0 0    ; b = 0 1 0 -1 . ( 4
)
It is worth noting that, as the incoming rate α has been set to zero, M is a 2 L × 2 L upper triangular matrix.

B. Exit time distribution

As the TASEP model is a random process, the time t taken by N particles to empty an L-site lattice is a random variable. We shall denote by p N,L (t) its probability density function (PDF). In the terminology of the previous section, the probability that the lattice is empty at time t is given by P 0 (t) where 0 = (0, . . . , 0) is the configuration where the L sites are empty. Now, the lattice is empty at time t if the N particles have evacuated it by a time τ ≤ t. Then, the probability Pr(τ ≤ t) that the exit time of the N particles is less than t is exactly equal to P 0 (t). The exit time PDF, p N,L (t) = dPr(τ ≤ t)/dt, is therefore given by p N,L (t) = Ṗ0 (t) = P (0,...,0,1) (t) ,

where the dot denotes the time derivative and where the last equality is obtained from the master equation (1).

It is thus sufficient to evaluate the probability that only site L is occupied at time t to obtain the exit time PDF.

Taking the Laplace transform of the master equation ( 1), we obtain the following algebraic system

(s -M )| P (s)⟩ = |P (0)⟩ , ( 6
)
where s is the parameter of the Laplace transform defined by

f (s) = ∞ 0 dt f (t)e -st . (7)
In Eq. ( 6), |P (0)⟩ is the initial probability vector with a single nonzero component: P (1,...,1,0,...,0) = 1 with N 1's and (L -N ) 0's. Solving the triangular algebraic system (6), one obtains the Laplace transform of P (0,...,0,1) (t) and, from there, the PDF p N,L (t) itself. Of particular interest is the Mean Exit Time (MET) of that distribution. In the next sections, we shall be mainly interested in the asymptotic behavior of this quantity as L becomes large while N remains finite. For that reason, we shall denote the mean exit time of N particles from a lattice {1, L} with L sites as T (f )

N (L). The superscript (f ) emphasizes the fact that this MET is obtained for a finite lattice with L sites and not for a section with L sites [1, L] embedded in an infinite lattice. We shall come back to that point in section IV B. Let us just note for now that T (f ) N (L) may be directly derived from the Laplace transform of p N,L (t) as

T (f ) N (L) = - dp N,L (s) ds s=0 . ( 8
)
A word is in order here. The solution of ( 6) is technically immediate, both because the system is triangular and because the results are rational fractions in s whose inverse Laplace transforms are straightforward. Nonetheless, it allows for an analytical determination of the exit time distribution p N,L (t) and its MET for small lattices only. The size of the Markov matrix grows indeed exponentially fast with L and we have not found any compact way to express the analytical solution in the general case (N ≤ L). Of course, starting from N particles, only configurations with at most N particles contribute to the dynamics of the system. The dimension of the Markov matrix reduced to these configurations is much smaller than 2 L : for instance for N = 2, the total number of configurations with at most two particles is L(L + 1)/2 + 1 which grows algebraically as L 2 /2 for large L. For L = 20, the reduced Markov matrix is then roughly 200×200 vs 10 6 × 10 6 for the full one. However, in spite of this drastic reduction, analytical expressions become very lengthy whenever L > 20 and, although exact, they are not particularly helpful in determining asymptotic behaviors for large L. They provide results that can be used as benchmarks for simulations though. Examples of such results for N = 2, 3 are provided in appendix A. In the next section, we shall therefore use a different method to tackle the determination of the MET for arbitrary large lattices.

III. ONE PARTICLE: BALLISTIC REGIME

We briefly treat here the exit time distribution of a single particle initially located on site 1 of an L-site lattice {1, L}. As it is more convenient, we switch from the "Eulerian" description based on particle configurations, that we have used so far to express probabilities, to a "Lagrangian" approach where particles are traced. Let us then call P (n; t) the probability that the particle lies on site n ∈ {1, . . . , L + 1} at time t. The addition of a virtual (L + 1)th site allows the particle to exit the lattice. This site is "absorbing" and P (L + 1; t) is then the probability that the lattice {1, L} is empty. According to Eq. ( 5), we then have p 1,L (t) = P (L; t) = Ṗ (L + 1; t). In the Lagrangian terminology, the master equation (1) translates into Ṗ (1; t) = -P (1; t) (9) Ṗ (n; t) = P (n -1; t) -P (n; t) , n ∈ {2, . . . , L} .

Taking the Laplace transform of equation ( 9) with an initial condition given by P (1; 0) = 1 (all other component being zero) yields P (n;

s) = (1 + s) -n . Hence, p1,L (s) = P (L; s) = (1 + s) -L , (11)
which upon inversion yields the exit time distribution of a single particle out of the {1, L} lattice,

p 1,L (t) = t L-1 (L -1)! e -t . ( 12
)
This distribution is of the Poisson type, as expected: a single particle indeed never experiences exclusion and spends on each site a time that follows the same exponential distribution (e -t ). Therefore, the total amount of time it spends on the lattice {1, L} is nothing but the sum of L exponentially distributed variables which leads to the Poisson distribution (12). Moreover, according to equations ( 8) and (11), the mean exit time of the particle is

T (f ) 1 (L) = L . (13)
The particle spends on average a unit of time on each site and thus travels at constant velocity. In that respect, its motion is ballistic. The purpose of the next section is to detail how this motion is hindered when another particle seats initially next to its right side.

IV. TWO PARTICLES : EXACT AND ASYMPTOTIC EXPRESSIONS FOR THE MET.

A. Finite lattice

Let us label the particles in their exiting order, namely from right to left, and consider first the same problem as in the previous section but with 2 particles initially located on site 2 (first particle) and on site 1 (second particle) of the lattice {1, L}. We can show (see appendix B) that the mean exit time of these two particles is exactly given by

T (f ) 2 (L) = L + L -1 4 L-2 × 2L -3 L -1 . ( 14
)
where m n = m!/(n!(m -n)!) is the binomial coefficient. Asymptotically, for large L, we then find

T (f ) 2 (L) = L + 2 √ π √ L + O(L -1/2 ) . (15)
Comparing this expression to the 1 particle MET (13) shows that the main effect of adding a particle next to the first one at the start of the process is to delay its exit by an amount of time that is proportional to the square root of the distance it has to travel to exit. In the next section, we shall interpret that result as a consequence of the random motion of the second particle confined on its right side by the random motion of the first one that it cannot overtake.

B. Infinite lattice

We now consider a problem related to the previous one although slightly different: what is the time T 2 (L) necessary for 2 particles to exit the section [1, L] of an infinite lattice, with the same initial positions as for the finite lattice {1, L}? This problem has much simpler boundary conditions than in IV A as particles keep moving on the infinite lattice instead of being absorbed at site (L + 1). This will enable us to develop a connection with a diffusion equation. The physical difference between the two situations is that in this present case, a particle having exited the [1, L] section still hinders the previous ones, whereas in the finite domain problem, the dynamics of a particle changes to a ballistic one each time its predecessor exits the lattice {1, L}. However in the large L limit, we expect the particles mean relative distances to become large, and the additional constraint provided by the following particles to be weak. Our following results will sustain this claim.

We first use results developed in [14], which provides an exact formula for the probability of 2 particles to be at positions X 1 and X 2 at times t knowing their initial positions Y 1 and Y 2 at time 0. From that we are able to deduce (see appendix C):

T 2 (L) = L + 2 √ π × Γ(L + 1/2) Γ(L) (16)
and this again yields the same asymptotic behaviour than Equation (15) :

T 2 (L) = L + 2 √ π √ L + O(L -1/2 ) ( 17
)
therefrom showing the equivalence of the finite and infinite formulation of the exit problem in the large L limit.

We now make a connection between this problem and a diffusion equation. The master equation for this two particles case writes [14]:

∂ t P (k 2 , k 1 ; t) = +P (k 2 -1, k 1 ; t) + P (k 2 , k 1 -1; t) -2P (k 2 , k 1 ; t) P (k, k + 1; t) = P (k, k; t) P (1, 2; t = 0) = 1,
or else 0 (18) valid for any k 2 < k 1 , the positions of the trailing and leading particles, respectively. The boundary condition elegantly accounts for the special case k 1 = k 2 + 1. The probability P(χ; t) that the distance between the 2 particles be χ at time t then follows from

P(χ; t) = k2≥1 P (k 2 , k 2 + χ; t) (19)
and satisfy

∂ t P(χ; t) = P(χ + 1; t) + P(χ -1; t) -2P(χ; t) (20)
with the boundary condition

J (χ = 0; t) ≡ P(0; t) -P(1; t) = 0 (21)
We notice that this is a discretized version of the diffusion equation with a no flux condition (J = 0) originating from the exclusion constraint, and this will enable us to develop an approach based on this equation in the next section. Solving equation (20) we obtain for the mean distance (see appendix D)

⟨χ(t)⟩ = e -2t 2 [(4t + 1)I 0 (2t) + 4tI 1 (2t)] + 1 2 (22)
where I k is the modified Bessel function of order k. Consequently, for large times t, the distance between the 2 particles varies like ⟨χ(t)⟩ ≃ 2 √ π √ t. In the context of the exit of 2 particles initially at k 2 = 1 and k 1 = 2, particle 1 reaches the end of the lattice after a time L, time at which particle 2 is on average at a distance ⟨χ(L)⟩ ∝ L 1/2 behind particle 1. Then particle 2 reaches the end of the lattice with a delay ⟨χ(L)⟩. Finally, for L large, the exit time of the 2 particles is T 2 (L) ≃ L + 2 L/π, which is consistent with our previous exact results.

V. THE DIFFUSION APPROXIMATION

The former calculations suggest the following simple physical approach : since the leading particle has on average a ballistic motion with a constant velocity, it is convenient to study the motion of the rear particles in the reference frame of the leading one. As we saw in equation ( 20) this leads to a diffusion equation for the motion of the second particle with a no flux boundary condition accounting for the exclusion. This can be generalized to any of the (N -1) trailing particles, the preceding particle acting as an impenetrable wall due to exclusion, to recursively find the average position of the n-th particle with respect to the leading one.

FIG. 2. Mean exit times of the first two particles.

We denote by x the relative position of a particle with respect to the leading one (particle 1), and X its absolute position in the lab frame. Let us first consider the occupancy probability P 2 (x; t) of the second particle. As we saw, this relative motion is simply described by a continuous diffusion equation. We are then left with the following set of equations in the domain x < 0 :

∂ t P 2 (x; t) = ∂ xx P 2 (x; t) , x < 0 P 2 (x; t) -----→ x→-∞ 0 P 2 (x; t = 0) = δ(x -0 -) (23) J (0; t) ≡ -∂ x P 2 (x = 0 -; t) = 0
Some comments are in order : here the space domain extends from x = 0 corresponding to the position of the leading particle up to x = -∞ when the trailing particle stays at rest in the lab frame. At t = 0 particle 2 is situated next to the leading particle, which in the continuous limit gives the stated initial condition. Finally the exclusion caused by the leading particle is described by a no flux condition J (x = 0) = 0. The solution of this set of equations is elementary and is twice the fundamental solution of the 1D diffusion equation. This immediately leads to the average position for the second particle with respect to the first one: ⟨x 2 (t)⟩ = -2 t/π. When the leading particle exits at a mean time T 1 = L the second one therefore sits at a position ⟨X 2 (L)⟩ = L -2 L/π in the lab frame. Since we demonstrated the equivalence of the finite and infinite lattice frames for the exit problem, we can assume particle 2 to be then unconstrained. Hence, it needs an additional time T 2 -T 1 = 2 L/π to exit, (see Figure2 for a pictorial view). This reproduces the result obtained previously by our exact algebraic computations (15), and this validates our continuous approach.

Encouraged by this first result, we seek a recursive scheme to obtain the mean position of the (n + 1)th particle with time, assuming an average position ⟨x n (t)⟩ = -β n √ t of the previous one, always relative to the first particle. The set of evolution equations for the (n + 1)-th particle can then be written as:

∂ t P n+1 (x; t) = ∂ xx P n+1 (x; t) , x < ⟨x n (t)⟩ P n+1 (x; t) -----→ x→-∞ 0 P n+1 (x; t = 0) = δ(x -0 -) (24) [∂ x P n+1 (x; t) + ⟨ ẋn (t)⟩P n+1 (x; t)] x=⟨xn(t)⟩ = 0
The last condition can be established for example by imposing d dt ⟨xn(t)⟩ -∞ P n+1 (x; t)dx = 0 which results from the normalization condition. It enforces a no flux condition at the moving boundary x n (t). The solution of Eqs ( 24) is simply

P n+1 (x; t) = 1 √ πt erfc(β n /2) e -x 2 /4t ( 25
)
One can check that both the normalization and the noflux condition, which are related, are satisfied by the above solution when the boundary is moving ∝ √ t. The average velocity d dt ⟨x n+1 (t)⟩ = -β n+1 /(2 √ t) can finally be calculated by taking the mean of the diffusion equation as d dt ⟨x n+1 (t)⟩ = -P n+1 (0; t), and this allows us to write the following recursion relation :

β n+1 = 2 √ π exp (-β 2 n /4) erfc(β n /2) ( 26
)
initiating at β 1 = 0. This is the central analytical result which allows to calculate the actual average position of the nth particle as a function of time in the lab frame as

⟨X n (t)⟩ = t -β n √ t.
In Figure 3, we have plotted the resulting mean trajectories of the first 5 particles in the lab frame as a function of time, (see dashed lines). For comparison we also plotted in solid lines the mean trajectories obtained by a direct Gillespie-type simulation of the TASEP over 1000 replicas and L = 2000 sites, see section VI. As is clear in the insert, the continuous diffusion approximation breaks down at short times, since TASEP particles can not have negative velocities in the lab frame. It however works remarkably well at larger times/positions, giving the exact large L asymptotic for N = 1 and 2 and an error of the order of a fraction of percent for small values of N (see below).

Finally, for N particles and large L the asymptotic behavior for the mean exit time is :

T N (L) ≃ L + β N √ L ( 27
)
where β n is given by Eq. ( 26).

FIG. 3. Mean trajectories ⟨Xn(t)⟩ = t -βn √ t of successive n = 1, 2, 3, 4, 5 TASEP particles from the continuous approach (dotted lines). Inset is a zoom near the origin. In comparison, the Gillespie-simulated trajectories in solid lines for 1000 different histories show that that the diffusion approach is quite good, see section VI (note that the artefact at the end of the simulation near X = 2000 is due to taking a mean position over a non constant ensemble of particles, since exiting particles are disappearing the simulation at L = 2000.

VI. GILLESPIE SIMULATIONS OF EXIT TIMES

Numerical simulations were also done to compute directly the exit times. In order to compare our results with real data, we simulated the emptying of a TASEP from an initial step condition using a continuous time Gillespie algorithm, see appendix E for details. The simulations were done with L = 300, 600 and 1000 sites, and with up to N = 50 particles. The rather modest number of copies was generally enough to ensure a reasonable error on the mean exit time, since this quantity is in itself a mean of different Gillespie times along one single history. The values of the coefficients β G N (L) were then computed using the definition

β G N (L) ≡ T G N (L) -L √ L (28)
and compared to the values obtained by our recursion equation ( 26) in figure 4, see the red-green-blue stars and black solid curve respectively.

FIG. 4. Values of β G

N obtained directly by Gillespie simulations (stars : red, green, blue L = 300, 600, 1000 respectively) compared to the values of βN obtained by the recursive scheme based on the diffusion equation approach, equation (26) (solid curve). The dashed line correspond to the limit of vanishing densities of reference [1], see section VII below.

The agreement between these two independent methods (diffusion approximation calculation and Gillespie simulations) is excellent for L = 1000 and for small values of N , with a vanishing relative error at N = 1 and 2 since the diffusion method gives the exact result there, and a relative error ranging from +1% for N = 3 to -0.7% for N = 10. This ultimately validates our diffusion continuous approach. We also note that for large N values, the Gillespie simulations are getting closer to our estimation (27) when increasing L, with an error of only -7% for N = 50 and L = 1000. Actually we conjecture that our diffusion scheme produces an exact result in the limit L → ∞ and N → ∞ with N/L → 0 as we discuss in the next section.

VII. DISCUSSION AND CONCLUSIONS

In this work, we have considered the question of the mean time taken by N particles to empty a lattice with L sites while being transported according to the rules of the TASEP model and starting from the leftmost sites of that lattice (step initial condition). We have investigated two slightly different versions of that problem: A) particles definitively exit the lattice as they leave the site L and B) particles keep moving along an infinite lattice after they have crossed the Lth site. For N = 2 particles, we have found the exact mean exit time for both problems and we have shown that they have a common asymptotic behaviour at large

L equal to T 2 (L) = L + 2 √ π √ L + O(L -1/2
). Still for N = 2 particles and within the framework of problem B, we have revisited that result by showing that the probability distribution of the distance between the particles obey a master equation that is a discrete version of a diffusion equation. From there, we have calculated the mean distance as a function of time and rederived the asymptotic behaviour of the mean exit time. Then, generalizing this approach to N ≥ 3 particles, we have devised an approximate diffusion model for the relative positions of consecutive particles that leads to a mean exit time for N particles that behaves for large L as T N (L) ∼ L + β N √ L where β N can be calculated recursively. Finally, we have confirmed the validity of this approximation by Gillespie simulations for values of N ≪ L.

Our diffusion model seems to work particularly well for a small finite number of particles. In the limit where the lattice size becomes infinite, L → ∞, the average particle density N/L tends therefore to zero. It is nonetheless tempting to try to extrapolate our results to a number of particles proportional to the lattice size, N = µL (as L → ∞) with a proportionality coefficient µ ≪ 1 in order to keep N ≪ L. Assuming equation (26) to be still valid for large values of N ≪ L, the asymptotic behaviour β N ∼ 2N 1/2 for large N obtained from Eq. ( 26) would yield

T N (L) ≃ L + 2 √ N L. (29)
Letting N = µL then provides the following asymptotic behaviour

T µL (L) ∼ (1 + 2 √ µ)L, (L → ∞, µ ≪ 1) . ( 30
)
This is to be compared to the exact known result T µL (L) = (1+ √ µ) 2 L obtained by Johansson for N = µL in the limit L → ∞ and µ finite (see ref. [1], theorem 1.6, Eq. (1.19) in which γ = 1/µ [30]). The corresponding value of β J N defined as in equation ( 28) reads β J N = (2 + √ µ) √ N . Using this result in the limit µ = 0 corresponding to our vanishing density regime, we have also plotted in Figure 4 the corresponding β J N = 2 √ N (black dashed line). We can see that our approximation gives a much better estimate of T L (N ) in the finite N regime and behaves decently at N large, with the same asymptotic value of β N . This leads us to conjecture that equation ( 29) is exact in the limit L → ∞ and N → ∞ with the density µ = N/L → 0, a region outside of the scope of ref. [1].

Another problem of interest is the exit time of N particles transported without exclusion. In that case, particles are all independent. They wait for a time t distributed according to the exponential distribution e -t between two consecutive jumps to the right, may overtake each other and occupy the same lattice site as others. The particle to last exit the lattice among the N , irrespective of its initial location, sets the exit time. Evaluating the distribution obeyed by the latter thus simply amounts to finding the distribution of the maximum of the individual exit times of each of the N particles (that depend on their initial location). For N = 2 particles starting respectively from sites 1 and 2 of an L-site lattice, it can be shown that the exit time asymptotically reads T

2 (L) ∼ L + 1 √ π √ L + O(1), for large L [31]. Strik-
ingly, we see that the √ L correction to the mean exit time of a single particle is not solely attributable to the exclusion effect of the TASEP model. It also occurs in independent particles as a by-product of the distribution of the maximum of their individual exit times, although with a different prefactor (half of the TASEP one for 2 particles). Preliminary analytical and numerical results seem to show that for a large number N of particles, all starting from site 1 of the lattice, the prefactor of the √ L correction of the exit time is proportional to ln N as N → ∞. This behaviour is to be contrasted with the √ N correction obtained in presence of exclusion for the TASEP model.

Finally, this study can be seen as a step towards the calculation of exit times of some more refined transport models. For example, one could try to test the diffusion approximation used in this paper to compute probabilities of interest studied in the clearance problem of [32]. Queuing problems [START_REF] Arita | [END_REF] or experimental microfluidic setups [33] could also benefit from our approach (e.g by relaxing the exclusion constraint for the queuing problem, or allowing for bidirectional transport like in the ASEP or SEP models, see [16,22]).

Appendix A: Some exact results for small lattices In table I, we list some exact results for the exit time distributions and their MET that can be obtained from the method exposed in section II B. As may readily be checked from the third column of this table, the Mean Exit time of N = 2 particles on a finite lattice {1, L} (L ≤ 10) agrees with the exact formula provided in (14). Laplace transforms of the time distributions have been given up to L = 5 only for they then become somewhat lengthy. From L ≥ 3 onwards, the denominator of p2,L (s) is (s + 1) L+1 (s + 2) 2L-5 . The constant coefficient of the numerator polynomial is 2 2L-5 and its highest degree coefficient (s L-3 for L ≥ 3) is the Catalan number C(L) = 2L L /(L+1) (valid for L ≥ 2). As for p3,L (s), its denominator is given by (s+1) L+2 (s+2) 2L-3 (s+3) 3L-14 for L ≥ 5.

Exact results for small lattices (up to L = 20) with N = 2, 3 particles are typically obtained by Maple on a basic laptop within less than a minute computation time. These results may serve as benchmarks for simulations. As expected, the time needed by the two particles to exit the section [1, L] of an infinite lattice is slightly longer than the time needed to exit the finite lattice {1, L} given that when the rightmost particle has gone out of {1, L}, the last one is free to move ahead while it can still be hindered by the front particle on the infinite lattice.

Appendix D: Exact mean relative distance for 2 particles

To obtain the mean relative distance between two particles, initially side by side, on an infinite lattice, we first take the Laplace transform of Eqs. (20) and ( 21): s P(1; s) -1 = -P(1; s) + P(2; s) s P(χ; s) = -2 P(χ; s) + P(χ + 1; s) + P(χ -1; s), where χ ≥ 2. Solving for P(χ; s) and taking into account the fact that χ≥1 P(χ; t) = 1 (i.e. χ≥1 P(χ; s) = 1/s), we obtain

P(χ; s) = 1 -λ s λ χ-1 , (D1)
where • for N ≤ 10 the relative error of equation ( 26) with respect to our best estimate of the exact β N , obtained with the highest copy number and the high-est L, is less than 1%, (0 for N = 1 and 2 since our expression is then exact), and reaches -7 % for N = 50.

λ = 1 + s 2 - 1 + s 2

FIG. 1 .

 1 FIG. 1. Initial, intermediate and final configurations.

2 - 1 .

 21 inverting that expression, we finally obtain⟨χ(t)⟩ = e -2t 2 [(4t + 1)I 0 (2t) + 4tI 1 (2t)] + 1 2 (D4)where I k is the modified Bessel function of order k:I 0 (x) = n≥0 x 2n /[4 n (n!) 2] and I 1 (x) = xI ′ 0 (x) where the prime denotes derivatives wrt x.

TABLE I .

 I Exact Laplace transform and Mean Exit Time of the exit time distribution of a finite lattice {1, L} for L ≤ 10 and N = 2, 3.

	L 2	p2,L(s) 1 (s + 1) 3	T 2 (L) (f ) 3	p3,L(s) n.a	T 3 (L) (f ) n.a
	3	2 (s + 1) 4 (s + 2)	9 2	2 (s + 1) 5 (s + 2)	11 2
	4	5s + 8 (s + 1) 5 (s + 2) 3	47 8	12s 2 + 39s + 32 (s + 1) 6 (s + 2) 5	233 32
	5	2(7s 2 + 21s + 16) (s + 1) 6 (s + 2) 5	115 16	110s 3 + 495s 2 + 751s + 384 (s + 1) 7 (s + 2) 7 (s + 3)	3409 384
	6 7 8 9 10	. . . . . . . . . . . . . . .	1083 128 2485 256 11195 1024 24867 2048 437075 32768	. . . . . . . . . . . . . . .	107617 10368 13237775 1119744 2132010983 161243136 254084494957 17414258688 7491745364599 470184984576

Appendix B: Exact MET for 2 particles, finite lattice

To find the Mean Exit time of two particles initially located on site 1 and 2 of the finite lattice {1, L}, it is sufficient, according to Eqs. ( 5) and (8), to find the Laplace Transform (LT) of the probability that particle 2 is on site L of that lattice while particle 1 has left it. We shall denote that quantity by Po (L; s) where the subscript o indicates that particle 1 has left the lattice. We shall denote by P (k 2 , k 1 ; s) the LT of the probability that particle 1 is at site k 1 and particle 2 at site

We shall now take advantage of the fact that P (n, L; t) is known exactly for it is the probability that two particles located on sites 1 and 2 at t = 0 be located at site n and L, respectively, at time t. This transition probability is provided by Schütz in [14] who has solved this problem on an infinite lattice. Yet, as none of the particles have left the section [1, L], this probability is exactly the same as for the finite lattice {1, L}. This makes the necessary connection between the finite and infinite lattice problems. According to [14], we have

where

and where F -1 (k; t) = F 0 (k; t) -F 0 (k + 1; t). Expanding P (n, L; t) in sums of products of exponentials and powers in t makes it easy to obtain its Laplace transform P (n, L; s). Reinstating the latter in Eq. (B2) and using

eventually yields, after a somewhat lengthy calculation,

Appendix C: Exact MET for 2 particles, infinite lattice

The easiest way to obtain the exact Mean Exit Time (MET) for 2 particles leaving a section [1, L] of an infinite lattice while being initially located on sites 1 and 2 of that section is probably to use the integral formula given by Rakos and Schütz [15] for the probability that the second leftmost particle of two initially side by side particles has carried out at least L steps to the right at time t. This probability, that is exactly the probability that the two particles have left the section [1, L] by time t, is given by

From Eq. (C1), the corresponding exit time distribution is given by p 2,L (t) = Ṗ (L, 2, t), whence the MET

The probability P (L, 2, t) may be evaluated in terms of incomplete gamma functions γ(n, t) = t 0 x n-1 e -x dx as

Using this expression and (C3), T 2 (L) can eventually be cast into the simple form

where Γ(L) = γ(L, ∞) is the complete gamma function.

If we compare the asymptotic expressions of T 2 (L) and T (f ) 2 (L) (the MET of two particles leaving a finite lattice {1, L} -see Eq. ( 14)), we find that they differ at order L -1/2 . More precisely,

Appendix E: Gillespie simulations details Numerical simulation were performed using Octave on a DELL XPS13. The continuous time Gillespie method was used in order to produce an in-silico realization of equation 1. In this method, each history simulate a stochastic trajectory associated with the TASEP master equation. Most of the simulations were done using 10 3 histories, in order to keep the computation time manageable on a laptop, especially for large values of N ≃ 50 and L ≃ 1000. To estimate the error for mean values such as T G N (L), we performed 20 independent simulations of 10 3 histories and obtained a dispersion of values of the order of ∆T ∼ 0.5 for T G N (L) ∼ 300 -1000. The same procedure was then used with 10 4 histories and, as expected, lowered this figure to ∆T ∼ 0.15.

The precision obtained with 10 3 copies was usually enough to compare the simulations results with our theoretical value β N . For small N < 5 however, the values of β N (L) estimated by the different methods and for different length L = 300 -1000 are very close, and it was necessary to use 10 4 copies to order the different values of β N (L) properly. It was found that :

• for fixed N the values of β G N (L) are systematically decreasing when increasing L, as seen in Figure 4.

• for values of N > 10 our result (26) underestimates the value of the coefficient, β N (L), while for small values it is overestimating.


Supravat Dey
email: supravat.dey@gmail.com

Gladys Massiera
email: gladys.massiera@umontpellier.fr

Estelle Pitard
email: estelle.pitard@umontpellier.fr

Role of cilia activity and surrounding viscous fluid on properties of metachronal waves

Large groups of active cilia collectively beat in a fluid medium as metachronal waves, essential for some microorganisms motility and for flow generation in mucociliary clearance. Several models can predict the emergence of metachronal waves, but what controls the properties of metachronal waves is still unclear. Here, we investigate numerically a simple model for cilia in the presence of noise on regular lattices in one-and two-dimensions. We characterize the wave using spatial correlation and the frequency of collective beating. Our results clearly show that the viscosity of the fluid medium does not affect the wavelength; the activity of the cilia does. These numerical results are supported by a dimensional analysis, which is expected to be robust against the model for active force generation, unless surrounding fluid influences the cilia activity. Interestingly, enhancement of cilia activity increases the wavelength and decreases the beating frequency, keeping the wave velocity almost unchanged. These results might have significance in understanding paramecium locomotion and mucociliary clearance diseases.

The emergence of phase-travelling waves in dense arrays of active beating cilia, known as metachronal waves, is a complex multiscale physics problem [START_REF] Gray | Ciliary movement[END_REF][START_REF] Sleigh | The Biology of Cilia and Flagella[END_REF][START_REF] Blake | [END_REF][4][5][6][7][8] and is nonequilibrium because of the internal activity-driven movements of cilia [4]. The active beating of each cilium arises from the sliding of microtubules by thousands of molecular motors, and the subsequent interaction with the surrounding fluid medium. The coupling of a large number of these oscillators lead to synchronized dynamics over larger length scales. Illustrations are abundant in nature with ciliary living systems differing by cilia assembly geometry, cilia activity, or properties of the surrounding fluid. In respiratory tissues, the continuous cleaning of our lungs is provided by cilia beating waves that generate mucus flow [START_REF] Wanner | American journal of respiratory and critical care medicine[END_REF][START_REF] Smith | Respiratory physiology and Nurobiology[END_REF]. For certain microorganisms such as paramecium, synchronized beating of cilia help in their efficient locomotion [START_REF] Lauga | The fluid dynamics of cell motility[END_REF]. The complexity of cilia active beating pattern and their interaction with each other through a complex environment makes it difficult to predict the emergent wave properties, despite recent theoretical and experimental advancements. Models of cilia arrays [START_REF] Lagomarsino | [END_REF][START_REF] Elgeti | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Golestanian | [END_REF][15][16], aim to identify the conditions required for such a coordinated state and to comprehend the physical parameters that govern the properties of the metachronal wave and the subsequent mucus transport. Several models have been proposed [START_REF] Lagomarsino | [END_REF][START_REF] Elgeti | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Golestanian | [END_REF]17], wherein the coupling is primarily described as a viscous hydrodynamic coupling. In these models, different types of active forces -from simple to complex, successfully generate continuous beating of a cilium. Numerical simulations enable to investigate the intricate structure of cilia by considering their beating as a filament bending wave [START_REF] Elgeti | Proc. Natl. Acad. Sci. USA[END_REF]18,19]. Another approach is to model cilia by actuated micron-sized beads called rotors [20][21][22][START_REF] Kanale | Proceedings of the National Academy of Sciences[END_REF] or rowers [START_REF] Lagomarsino | [END_REF][START_REF] Kotar | [END_REF]25]. For a large group of cilia arranged in arrays, it has been shown that hydrodynamic coupling can lead to metachronal waves for various models of cilia [START_REF] Elgeti | Proc. Natl. Acad. Sci. USA[END_REF]15,16].

Recently, the influence, on these collective behaviors, of several physical parameters such as noise [25,26] and disorder in the arrangement and orientation of cilia has been investigated both numerically [25,27,28] and experimentally [28,29], showing that spatial heterogeneity favors transport. Other important physical quantities that may play a role on the coordination are the activity and the dissipation, that will have opposite impacts on the metachronal waves emerging from cilia beating. Experimentally, a decrease in beating frequency with viscosity was found [30,31], whereas the beating amplitude and the metachronal wavelength were found constant up to ≈ 50 times the viscosity of water [32,33]. Theoretically, the mutual influence of activity and dissipation were almost not explored [34]. Here, our fundamental inquiry pertains to the interplay between cilia activity and fluid medium and its impact on the overall properties of metachronal waves.

To investigate this, we study the metachronal waves in the rower model of cilia in viscous fluid for one and two-dimensional regular lattices in the presence of thermal noise. In the rower model [START_REF] Lagomarsino | [END_REF][START_REF] Kotar | [END_REF]25], the complex active beating of a cilium is simplified into the back and forth motion, along an axis, of a micron-sized bead immersed in a viscous fluid, thus ensuring a low Reynold's number regime. Such an oscillating motion is driven by two harmonic potential branches, corresponding to the stroke and anti-stroke of the cilia beating, with a geo-

(A) y V(y) σ = 1 σ = -1 -𝒜 𝒜 μ × × × × × × × × × × × × × × × × rij = rj -ri r i r j (L × L) (C) (B)
FIG. 1. Rower model of cilia [START_REF] Lagomarsino | [END_REF]. (A) The motion of a micron-sized bead in a viscous medium under two harmonic potential branches, corresponding to σ ± 1, represents the stroke and anti-stroke beating of a cilium. (B) The bead switches branches when it reaches terminal position y = ±A. (C) Rowers beating along y-axis on two dimensional L × L regular lattices. Hydrodynamic interaction between rower i and j is modelled by the Oseen coupling that depends on the separation vector rij and viscosity of the medium. metric switching mechanism. The bead moves downhill of a potential until it reaches one of the two terminal positions for which switching to the second branch occurs (Fig. 1A,B). This switching is like pumping energy to lift the bead on the upper side of the other potential at the terminal. At a given time, the bead can be found in one of these two states, the stroke and anti-stroke of the cilia beating, represented by a discrete σ = ±1. The driving force for a bead displacement y for a given σ can be written as

f (y, σ) = -dV (y,σ) dy = -k(y -σµ/2), (1)
where k is the force constant associated with the harmonic potentials, µ is the distance between minima of two potentials, and A is the beating amplitude. The supply of energies during each downhill motion in a harmonic potential, kA 2 /2, and during each switch, the pumping energy kµA, keep the bead oscillating in the dissipating media. Therefore, for a given µ, the 'activity' of the bead depends on values of k and A. Because of its simplicity and ability to capture the two-stroke beating of cilia, the rower model has become a method of choice for theoretical and experimental studies of synchronization in ciliary systems [16,29].

We consider a system of N rowers beating in the y direction in a viscous medium. Rowers are placed regularly in one-or two-dimensions (square) lattices (see Fig. 1C) at fixed positions r i (for i = {1, 2, 3, ..., N }). The displacement, y i , of a rower i is hydrodynamically coupled with the others and is given by where γ = 6πηa is the viscous drag coefficient for a bead with radius a moving in a fluid medium with viscosity η and O(i, j) the coupling strength between rower i and j. In the far-field hydrodynamic coupling approximation, for which both the distance from the surface and the distance between two adjacent rowers (lattice spacing ℓ) are large compared to a, O(i, j) is set by the Oseen tensor:

dy i dt = - f i γ + j̸ =i O(i, j) f j + ξ i , (2)
O(i, j) = 1 8πηrij 1 + ( yij rij ) 2
, with i ̸ = j, and r ij = r j -r i , the separation vector between rowers i and j. The last term is due to the thermal noise, obeying the following delta-correlation: ⟨ξ(t)⟩ = 0, ⟨ξ(t 1 )ξ(t 2 )⟩ = 2 D δ(t 1 -t 2 ). For simplicity, we assume no correlation between the noise acting on each of the rowers as in [25]. The noise strength or diffusivity is equal to D = k B T /γ, k B and T being the Boltzmann constant and the temperature. The displacement of a single isolated bead shows sustained oscillations with the frequency ν 0 = 1/(2 τ d log [(µ + 2A)/(µ -2A)]), where τ d = γ/k is the relaxation time for the bead to reach equilibrium in a harmonic potential [START_REF] Kotar | [END_REF]. Such two coupled rowers beat collectively with antiphase synchronization [START_REF] Lagomarsino | [END_REF]. For many rowers, the interplay between the activity of the rowers and coupling through the medium generates metachronal waves [START_REF] Lagomarsino | [END_REF].

Simulation details -The Euler method with an integration step equal to 5×10 -3 s has been used to evolve the coupled dynamical equation (Eq. 2), starting from random initial values for {σ i , y i }. The open boundary condition is implemented. Parameters are chosen within the experimentally relevant range [START_REF] Kotar | [END_REF]35], as follows: a=1.5 µm, ℓ=8 µm, k=2.6 pN •µm -1 , µ=2 µm, A=0.56 -0.8 µm, η=2 -20 mP a•s, and T =300K. Results presented here for large system sizes; N =L=200 (for 1d) and N =L 2 =1600 (for 2d). Comparing results with smaller systems (not shown here), we confirm that the presented results have no system size dependence.

Results -Fig. 2 shows the metachronal waves in the onedimensional lattice. The beads' displacement against the rowers' position for a given time displays two spatial waves that are visualized by connecting displacements y i by lines for all rowers at the even and odd lattice sites separately (Fig. 2A) in agreement with [START_REF] Lagomarsino | [END_REF]36]. This is a unique feature of the rower model, and arises due to a degree of anti-phase synchronization between two adjacent rowers. The wave propagation is illustrated in Fig. 2C by the kymograph obtained for σ i (t). To characterize it, we compute the spatial correlation function between two rowers as a function of their separation vector r:

C(r) = ij ⟨σ i (r i , t)σ j (r j , t)⟩δ(r -r ij ) ij δ(r -r ij ) . (3)
As the rowers are placed on a regular lattice with lattice spacing ℓ, the coordinates of r are discrete and can be written as (mℓ, nℓ) with m, n ∈ {0, 1, 2, ..., L}. The measurement is done after a large equilibration time t 0 , where the system is assumed to reach a steady state. Brackets ⟨•⟩ represent average over times and ensembles. An ensemble is the collection of 5000 sets of {σ i (t)} recorded every 2 seconds after t 0 =2500 seconds. Fig. 2B shows the variation of C(x) = C(x, y = 0) in one dimension. For odd and even m values, two oscillating curves decay to zero as the distance between rowers x = mℓ increases. While the oscillations indicate the wave nature of the collective beating, the loss of correlations at larger x suggest a damping in the coordination on a characteristic length scale l d . C(x) can be fitted with the simple function ±e -x/l d cos(2πx/λ), the + (-) sign being for even (odd) m. This fit estimates the wavelength λ and decay length l d . For Fig. 2, λ ≃ 13.7ℓ and l d ≃ 9.0ℓ. In a recent work, wavelength, and decay length were measured experimentally for metachronal waves on the human bronchial epithelium, and these two lengthscale values are comparable [37]. Our results are consistent with the experiment. The ensemble and spatial average of the beating frequency was computed: ν ≃ 3.4 Hz and combined with λ to infer the metachronal wave velocity V = νλ ≃ 370 µm.s -1 . These values are consistent with estimates that can be inferred directly from the slopes in the kymograph Fig. 2C.

We then investigate the effect of viscosity of the fluid medium and activity of the cilia on the metachronal waves quantities: λ, l d , ν, and V =νλ, by computing C(x) for various η and A. The plot of C(x) for different A values shows that both the wavelength λ and the decay length l d increase with A (Fig. 3A andC), whereas the averaged beating frequency ν decreases with A, keeping the velocity of the wave V , almost constant (Fig. 3C andE). As ν 0 , the natural frequency of a single rower, decreases with A, the decrease of ν is expected. On the same line, increasing A, which is a characteristic length of the problem, may naturally increase the length scale of the emerging collective dynamics. Thus the respective variation of ν and λ can be generally expected. What is remarkable though is that they compensate to result in an almost constant metachronal wave velocity. Interestingly, C(x) does not depend on the values of η (see Fig. 3B), meaning λ and l d are independent of η and implying that the spatial behavior of emergent waves does not depend on the fluid viscosity and are only determined by cilia activity parameters, in agreement with experimental observations [30,33]. Below, we argued that such behavior is generally characteristic of a hydrodynamically coupled system. Finally, the frequency ν decreases as a function of η and so does V , as measured experimentally in [30,33].

On a square lattice, we find that the metachronal wave propagates along the beating direction y whereas no wave is obtained in the perpendicular direction (Fig. 5), suggesting longitudinal waves. In Fig. 4A andB, we plot the correlation function along y-direction C(0, y = nℓ) against n for various values of A and η. Similar to 1d, two spatial waves can be seen for even and odd values of n for a given parameter set. For a fixed η, λ increases and ν decreases with A, keeping the wave velocity V almost constant (Fig. 4C andE). On the contrary, λ remains constant and ν decreases with η, leading to a decrease in the velocity V with η (Fig. 4D andF). These results are consistent with the 1d results. We further note that although the qualitative behavior of metachronal waves in 1d and 2d are similar, the values of λ and V are relatively larger in 1d. This result raises interesting questions on the implications of the geometry of realistic ciliated tissues, which are mostly organized in 2d groups of cilia bundles.

In the direction perpendicular to beating, no oscillation is obtained (Fig. 5). C(x, 0) either monotonically decays to zero as for large A or shows a negative correlation for small y = mℓ with odd m that eventually approaches zero for large m. For a given A, C(x, 0) does not depend on η (Fig. 5B), although odd and even m can follow different curves, reminiscent of C(0, y). We compare the decay lengths of correlations along x and y directions l d,x and l d,y . The decay length for the damped oscillations along y, l d,y , can be estimated from the fitting method discussed above. The decay length l d,x is estimated from the exponential fit of the C(x = mℓ, 0) for even m values. The ratio l d,y /l d,x is plotted in the insets, one notes that l d,y /l d,x ≳ 2. For a fixed A, it remains unchanged with η. However, for a given η, the ratio increases with cilia activity A, which implies an enhancement of coherence along the beating direction compared the perpendicular one. This anisotropic response may be related to the anisotropy of the interaction strength. Indeed, considering the same r ij value, O(i, j) is two times larger along the y-axis than along the x-axis.

The fact that we obtain metachronal waves with spatial properties unaffected by viscosity has not been emphasized by previous studies, to our knowledge. Nevertheless, this remarkable numerical observation is robust on an order of magnitude of η obtained with both 1d and 2d simulations. To rationalize this result, one needs to look into the details of characteristic length, and timescales of the system set by the activity and the surrounding viscous medium. The relaxation time τ d = 6πηa/k for the bead motion in the viscous medium under a harmonic driving potential, which also determines the natural frequency ν 0 of a single rower, is a crucial timescale in our problem. In the rower model, there are two lengths scales, the amplitude A, and µ, the distance between the two branches of the potential (Fig. 1). Since we only vary A, we chose it as the typical length scale. We note that our conclusion below, however, does not depend on the choice of the length scale. Multiplying both sides of Eq. 2 by τ d /A leads to an adimensional equation for the collective beating dynamics:

dy ′ i dt ′ = -f ′ i - i̸ =j 3 a 4r ij 1 + (y ij /r ij ) 2 f ′ j + ζ i (t ′ ), (4)
where t ′ and y ′ are dimensionless time t ′ = t/τ d and displacement y ′ = y/A, f ′ j = -(y ′ j -µσ j /(2A)) is the dimensionless force acting on rower j, and ⟨ζ(t

′ 1 )ζ(t ′ 2 )⟩ = 2k B T /(kA 2 )δ(t ′ 1 -t ′ 2 )
is the adimensional noise correlation. As activity parameters A, k, and µ are constant, Eq. 4 is η independent. The latter suggests that the spatial properties are independent of η. However, as τ d is affected by η, it impacts the dynamical properties of the system. If any parameter A, k, and µ are influenced by the medium, then our observation will break down.

We argue that the independence on fluid viscosity of the spatial emergent properties is more generic to systems operating at low Reynold's numbers, irrespective of the model details. For systems at a low Reynold's number, both the viscous drag and hydrodynamic coupling between two objects are inversely proportional to η. The thermal noise strength D is also inversely proportional to η. As τ d is inversely proportional to η, this also means that the normalization of Eq. 2 by τ d will lead to the same conclusion. Therefore, one can get a similar adimensional equation as Eq. 4 for any model of active cilia coupled by a viscous fluid at low Reynolds numbers. Hence, the spatial properties of the emergent waves are expected to be independent of viscosity. If active beating is strongly dependent on the fluid rheology, then the active parameters of the rower model would depend on the viscosity, and this result would fail. Experimental results [31] seem to indicate though a very small dependence of cilia beating amplitude with the liquid medium. We note that additional sources of deviation from this result could be the viscoelastic nature of the fluid or a non-thermal noise, which we have not addressed in this paper. Finally, although this simple dimensional analysis cannot predict the occurrence or nature of emergent behavior, it is powerful in predicting that the wavelength or other spatial properties will be viscosity independent in general.

In conclusion, we have presented simple generic results about complex dynamics of hydrodynamically coupled model cilia. Using a very simple rower model of coupled oscillators, we have focused our study on the influence of activity and dissipation on the spatial and temporal synchronization properties of cilia assemblies. Enhancement of cilia activity increases the wavelength and beating period, keeping the wave velocity almost unchanged. On the other hand, viscosity does not affect spatial patterns characterizing metachronal waves such as wavelength or correlation lengths. On the contrary, the beating frequency and the wave velocity indeed decrease with viscosity. The deviation from such a behaviour may indicate the influence of additional properties of the medium not taken into account in the coupling description of the current model and could also be a signature of viscoelasticity or elasticity of the tissue itself. This could pave the way to the study of the emergence of specific functions of cilia in pathological contexts for example [START_REF] Kanale | Proceedings of the National Academy of Sciences[END_REF].

FIG. 2 .

 2 FIG. 2. Metachronal waves in 1d. A. Snapshot of displacements of the first 100 rowers. The displacement of even (odd) sites is plotted in light (dark) color. B. Correlation function C(x=mℓ) between two rowers is plotted against separation distance x. C. Kymograph of the beating state σ, with white color (black) representing σ=1 (σ= -1). Parameters: A = 0.8, η = 6 mPa•s, and N = 200.

FIG. 3 .

 3 FIG.3. Effect of viscosity and beating amplitude on metachronal waves in 1d (N =200). Correlation function C(x=mℓ) as a function of the distance between rowers, for three different A values (for η=6 mP a.s) (A) and for three different η values (for A=0.8 µm) (B). The wavelengths of metachronal waves λ are plotted against A (C), and η (D), together with the corresponding beating frequency ν. Propagation velocity V =λν is plotted against A (E), and η (F).

FIG. 4 .

 4 FIG.4. Effect of viscosity and beating amplitude on metachronal waves in 2d (N =L 2 =1600). Correlation function along beating direction C(x=0, y=nℓ) as a function of the distance between rowers, for three different values of A (for η=6 mP a.s) (A) and for three different η values (for A=0.8 µm) (B). The wavelengths of metachronal waves λ are plotted against A (C), and η (D), together with the corresponding beating frequency ν. Propagation velocity V =λν is plotted as a function of A (E), and η (F).

BFIG. 5 .

 5 FIG. 5. Effect of viscosity and beating amplitude on correlations along the direction perpendicular to the beating direction C(x=mℓ, y=0). A. C(x=mℓ, y=0) for three different A for a fixed η=6 mPa•s (A), and for three values of η for a given A=0.69µm (B).