S Chateau

J Favier

U D'ortona

S Poncet

TRANSPORT EFFICIENCY OF METACHRONAL WAVES IN 3D CILIA ARRAYS IMMERSED IN A TWO-PHASE FLOW

The present work reports the formation and the characterization of antipleptic and symplectic metachronal waves in 3D cilia arrays immersed in a two-fluid environment, with a viscosity ratio of 20. A coupled lattice-Boltzmann -Immersed-Boundary solver is used. The periciliary layer is confined between the epithelial surface and the mucus. Its thickness is chosen such that the tips of the cilia can penetrate the mucus. A purely hydrodynamical feedback of the fluid is taken into account and a coupling parameter α is introduced allowing the tuning of both the direction of the wave propagation, and the strength of the fluid feedback. A comparative study of both antipleptic and symplectic waves, mapping a cilia inter-spacing ranging from 1.67 up to 5 cilia length, is performed by imposing the metachrony. Antipleptic waves are found to systematically outperform sympletic waves. They are shown to be more efficient for transporting and mixing the fluids, while spending less energy than symplectic, random, or synchronized motions.

Introduction

Cilia and flagella are contractile hair-like structures, put into motion by biochemical energy, and protruding on the free surface of eukaryotic or prokaryotic cells. While prokaryotes are single-celled organisms, eukaryote cells have membrane-bound organelles and are found in every mammals. Many living organisms, going from the prokaryotic bacteria to mammals, use ciliary and/or flagellar propulsion as a swimming mechanism. Usually, flagella are external appendices used by micro-swimmers such as the algae Chlamydomonas reinhardtii for locomotion purposes, while cilia are generally internal appendices, shorter and more numerous, used for moving materials such as nutrients, dusts, or proteins into living organisms. Ciliary propulsion is a universal phenomenon, and many examples could be cited. For the particular case of the human body, cilia are responsible for the left-right asymmetry of the heart in the early embryonic development, for the transport of nutrients in the brain, and for the transport of mucus in the mucociliary clearance process which is the background of the present work (see [START_REF] Satir | Overview of structure and function of mammalian cilia[END_REF] for a review about the structure and function of mammalian cilia).

During the breathing process, a large number of foreign particles (bacteria, dusts, pollutants or allergens) can penetrate the organism. The human body has then developed three mechanisms to protect itself from these particles: coughing, alveolar clearance and mucociliary clearance which occurs on the epithelial surface of the respiratory system. In order to trap the particles, a layer of fluid called the Airways Surface Liquid (ASL) is covering the epithelial surface. Due to differences in the concentration of mucins inside the ASL, it is generally assumed to be the superposition of two different layers: the periciliary layer (PCL) and the mucus. The 7 µm depth PCL is located between the epithelial surface from where the cilia protrude and the mucus phase just above it. Mainly composed of water and of a few mucins of low molecular weight, this phase, often considered as being a Newtonian fluid similar to water, is a kind of lubricant allowing the mucus to slip on it and the cilia to beat without too much viscous resistance.

The mucus is a highly non-Newtonian fluid, with a strong visco-elastic behaviour. Additionally, the inner structure of the macromolecules composing it confers to the mucus a clear thixotropic behaviour, which is currently being studied and characterised to understand its inner rheological properties which exhibit a large variability [START_REF] Lafforgue | Rheological characterization of macromolecular colloidal gels as simulant of bronchial mucus[END_REF]. Mucus is indeed composed of 95% of water, but also contains macromolecules called the mucins [START_REF] Lai | Micro-and macrorheology of mucus[END_REF]. It serves as a physical barrier against infectious agents and dusts, but also to humidify the air flowing into the respiratory system and to catch the particles. Its height varies between 5 to 100 µm depending on many factors including the position in the respiratory system [START_REF] Widdicombe | Regulation of human airway surface liquid[END_REF], the pathology for a particular person, and the quality of the air inhaled. Its viscosity can also vary by several orders of magnitude within the same day [START_REF] Kirkham | Heterogeneity of airways mucus: variations in the amounts and glycoforms of the major oligomeric mucins MUC5AC and MUC5B[END_REF].

Cilia are the other protagonist of the mucociliary process. They are organized as tufts (around 200 to 300 cilia per tuft) at the epithelium surface. Formed by 9 pairs of microtubules placed regularly in a circle and 1 pair of microtubule at the center (which form the so-called axoneme), their purpose is to propel the surrounding fluid layers. Their diameter varies from 0.2 to 0.3 µm and their length from 6 to 7 µm [START_REF] Sleigh | The propulsion of mucus by cilia[END_REF].

The cilia motion can be decomposed into two steps: a stroke phase and a recovery phase. The stroke phase is characterized by almost straight cilia orthogonal to the flow in order to maximize the pushing effect; while the cilia are bending during the recovery phase in a more inclined plane to get closer to the epithelial surface in order to minimize the viscous resistance, and therefore to reduce their impact on the flow. During the stroke phase, which takes around 1/3 of the total beating period, the tips of the cilia enter the mucus phase [START_REF] Widdicombe | Regulation of human airway surface liquid[END_REF]. Their beating frequency is estimated to vary between 10 and 20 Hz. Note that the spatial asymmetry is essential for the cilia to generate propulsion in creeping flows, while the temporal asymmetry (recovery phase longer than stroke phase) is not necessary to induce a mucus motion [START_REF] Khaderi | Breaking of symmetry in microfluidic propulsion driven by artificial cilia[END_REF].

It has been experimentally observed that cilia usually do not beat randomly [START_REF] Sleigh | The biology of Cilia and Flagella[END_REF], but adapt instead their beatings accordingly to their neighbours, giving birth to the so-called metachronal waves (MCW) observed at the tips of the cilia. MCW occur when adjacent cilia beat with a constant phase lag ∆Φ between each other. For 0 < ∆Φ < π, the MCW move in the opposite direction as the fluid propelled and are called antipleptic MCW. On the contrary, for -π < ∆Φ < 0, the MCW move in the same direction to the flow, and are called symplectic MCW. When the phase lag ∆Φ is null, all cilia beat in a synchronized way. Finally, when the phase lag between neighbouring cilia is ∆Φ = ±π, a standing wave appears.

The universality of ciliary propulsion has intrigued scientists for decades, and several studies were conducted in order to understand it. One actual objective is to be able to mimic this process in order to, for instance, create cilia-based actuators for mixing, use them as flow-regulator in microscopic biosensors, or as micropumps for drug-delivery systems [START_REF] Li | Dynamics modeling and analysis of a swimming microrobot for controlled drug delivery[END_REF][START_REF] Chen | Magnetically actuated artificial cilia for optimum mixing performance in microfluidics[END_REF]. Moreover, diseases such as asthma or Chronic Obstructive Pulmonary Disease (COPD) are related to the mucociliary clearance process [START_REF] Gardiner | The importance of being cilia[END_REF]. The objectives are to understand the underlying mechanism that allows hundreds of cilia to act as a whole for the transport of mucus and how it affects the flow generated. The results could bring a deeper understanding in such pulmonary diseases.

Numerically, Gueron and co-authors [START_REF] Gueron | Cilia internal mechanism and metachronal coordination as the result of hydrodynamical coupling[END_REF][START_REF] Gueron | Energetic considerations of ciliary beating and the advantage of metachronal coordination[END_REF] showed in the 90's that two neighbouring cilia beating randomly will quickly synchronize within a few beating cycle due to hydrodynamic interactions and form antipleptic MCW. This has been recently confirmed by [START_REF] Elgeti | Emergence of metachronal waves in cilia arrays[END_REF], who observed the formation of MCW in a 3D single-phase environment. A certain degree of freedom in the beating pattern is required as shown in the theoretical work of [START_REF] Niedermayer | Synchronization, phase locking, and metachronal wave formation in ciliary chains[END_REF]. They modelled the beating pattern of cilia by circular trajectories. By allowing some flexibility in the radii, they managed to introduce the coupling leading to MCW formation. Among the different models used for the study of ciliary propulsion, the envelope model [START_REF] Taylor | Analysis of the swimming of microscopic organisms[END_REF][START_REF] Reynolds | The swimming of minute organisms[END_REF][START_REF] Tuck | A note on a swimming problem[END_REF][START_REF] Brennen | Fluid mechanics of propulsion by cilia and flagella[END_REF]Blake 1971a,b) assumed that cilia are so densely packed that it is possible to consider their tips as an oscillating surface. Nevertheless, such configuration has only been observed in nature for symplectic metachrony. Moreover, this technique is limited to small amplitude oscillations and impose no slip and impermeability conditions at the oscillating surface. In the sublayer (or stokeslets) model [START_REF] Keller | A traction-layer model for ciliary propulsion[END_REF][START_REF] Lighthill | Flagellar hydrodynamics[END_REF][START_REF] Phan-Thien | A boundary-element analysis of flagellar propulsion[END_REF][START_REF] Blake | A model for the micro-structure in ciliated organisms[END_REF][START_REF] Gueron | Cilia internal mechanism and metachronal coordination as the result of hydrodynamical coupling[END_REF][START_REF] Gueron | Energetic considerations of ciliary beating and the advantage of metachronal coordination[END_REF][START_REF] Niedermayer | Synchronization, phase locking, and metachronal wave formation in ciliary chains[END_REF][START_REF] Gauger | Fluid transport at low Reynolds number with magnetically actuated artificial cilia[END_REF][START_REF] Smith | Discrete cilia modelling with singularity distributions: Application to the embryonic node and the airway surface liquid[END_REF][START_REF] Ding | Mixing and transport by ciliary carpets: a numerical study[END_REF]) the cilia are modelled by a distribution of stokeslets, which impose a force on the surrounding fluid. A proper mirror image of the stokeslets is required to impose the no-slip condition on the surface from where the cilia protrude. However, the presence of a wall is known to alter the nature of the far-field of the stokeslets, and it is thought that this can have important consequences on the hydromechanics of the cilia near the wall [START_REF] Blake | Fundamental singularities of viscous flow[END_REF]. Moreover, stokeslets can only be used for fluids with constant viscosities. Results using this method tend to show that symplectic metachrony would be more efficient for the mucus transport than synchronously beating cilia, and that antipleptic metachrony would induce the lower flow rate. The opposite was nevertheless obtained by [START_REF] Gauger | Fluid transport at low Reynolds number with magnetically actuated artificial cilia[END_REF] who found that antipleptic metachrony was more efficient than symplectic metachrony for a particular cilia spacing. For a distance between two cilia of 1.5L, L being the length of the cilia, they obtained an increase in the pumping performance of 40% relative to a single cilium. However, while the beating pattern used in [START_REF] Gauger | Fluid transport at low Reynolds number with magnetically actuated artificial cilia[END_REF]) was realistic, they considered a slow stroke phase and a quick recovery phase which is the opposite of what is observed in nature. Some authors [START_REF] Gueron | Cilia internal mechanism and metachronal coordination as the result of hydrodynamical coupling[END_REF][START_REF] Gueron | Energetic considerations of ciliary beating and the advantage of metachronal coordination[END_REF][START_REF] Kim | Pumping fluids with periodically beating grafted elastic filaments[END_REF] showed that the energy spent by a cilium would decrease in the presence of metachronal motion. Others tried to model the internal axoneme of a cilium. Among them, [START_REF] Mitran | Metachronal wave formation in a model of pulmonary cilia[END_REF] used an overlapping fixed-moving grid formulation, coupling the finite volume and the finite element methods, to study the emergence of MCW. His model was very detailed and has lot of assets (two-layer flow, viscoelastic mucus) but required many experimental parameters and the impact of the MCW on the flow has not been considered so far. [START_REF] Niedermayer | Synchronization, phase locking, and metachronal wave formation in ciliary chains[END_REF] observed MCW formation in 1D cilia arrays. Those waves were stable only if the wavelength was 4 times higher than the cilia spacing. For the interested reader, a clear review of those computational modelings of the internal axoneme can be found in [START_REF] Fauci | Biofluidmechanics of reproduction[END_REF]. A different approach is to let the cilia adapt their motion in order to find the most energetically efficient ciliary beating pattern. In that context, [START_REF] Eloy | Kinematics of the most efficient cilium[END_REF] and [START_REF] Lauga | Shape of optimal active flagella[END_REF] computed the shape and energyoptimal kinematics of cilia from an energetic point of view. They found that the optimal kinematics strongly depend on the cilium bending rigidity, and closely resemble the twostroke ciliary beating pattern observed in natural cilia. Similarly, [START_REF] Osterman | Finding the ciliary beating pattern with optimal efficiency[END_REF] computed the ciliary beating pattern with optimal pumping efficiency of isolated cilia and arrays. Recently, new methods have been introduced such as the immersed boundary method (IB) used by [START_REF] Dauptain | Hydrodynamics of ciliary propulsion[END_REF] to model the swimming of the pleurobrachia. [START_REF] Lukens | Using lagrangian coherent structures to analyze fluid mixing by cilia[END_REF] used an IB method to study the mixing produced by a carpet of cilia in the context of the mucociliary clearance process. A coupled immersedboundary lattice Boltzmann method (LBM) was also used by [START_REF] Sedaghat | Numerical simulation of muco-ciliary clearance: immersed boundary lattice Boltzmann method[END_REF] to study several parameters in a 2D configuration using an Oldroyd-B model for the mucus rheology. They found that the transport of mucus was maximized when considering the mucus as a Newtonian fluid. While many studies regarding the emergence of MCW were conducted, only few addressed 3D configurations of the mucociliary clearance process. Among them, [START_REF] Elgeti | Emergence of metachronal waves in cilia arrays[END_REF] managed to observe symplectic and laeoplectic (perpendicular to the power-stroke direction) MCW formations. [START_REF] Ding | Mixing and transport by ciliary carpets: a numerical study[END_REF] did not study the emergence of MCW but performed a comparative study of antipleptic metachrony versus symplectic metachrony in terms of transport efficiency and mixing. Their results showed that both antipleptic and symplectic MCW enhanced the fluid transport and the mixing, with the antipleptic waves being the most efficient ones. However, the two aforementioned studies only considered a single fluid layer.

In this work, by using the solver developed by [START_REF] Li | An improved explicit immersed boundary method to couple with lattice Boltzmann model for single-and multi-component fluid flows[END_REF] and already validated in similar configurations, the focus is placed on the emergence of MCW in a 3D two-phase flow configuration with a viscosity ratio of 20 and a large number of cilia. In particular, it is shown that a simple hydrodynamical feedback based on mechanical concepts can trigger the emergence of both symplectic and antipleptic MCW, while usually only a single layer of fluid is considered and only antipleptic MCW are seen to emerge. From a numerical point of view, the main advantage of the present method is its ease of implementation. Additionally, the local character of the collisions in the LBM allows an easy and straightforward parallelization of the code, making the simulation of a large number of cilia possible. The numerical method also possesses the following advantages: (i) viscosity ratios up to O (10 2 ) can be achieved [START_REF] Porter | Multicomponent interparticle-potential lattice Boltzmann model for fluids with large viscosity ratios[END_REF], and (ii) the mucus-PCL interface emerges intrinsically from the model. To the knowledge of the authors, this solver is the only one that combines all these capabilities. The main contribution of this work is the thorough analysis of the advantages of the antipleptic and symplectic MCW over synchronized beatings by computing appropriate transport and efficiency ratios, mapping an inter-cilia spacing going from 1.67 to 5 cilia length. Finally, and for the first time, both the PCL and the mucus layer have been taken into account in the present study. The inherent advantages of the MCW for flow transport are studied by (i) considering the efficiency of the waves to transport the flow, (ii) comparing the flow generated and the volume transported by the different kinds of metachrony, (iii) comparing the energetic cost of the MCW, and (iv) analyzing the capacity of the waves to transport particles from an energetic point of view. The following analysis shows that antipleptic MCW are always the most efficient to transport mucus.

The remainder of this manuscript is organized as follows. In Section 2, details about the numerical method are given, and the different quantities used are introduced. In Section 3, the results are presented starting from the emergence of MCW by considering the fluid feedback onto the cilia; and then a parametric study is done to quantify the impact of both the MCW and the cilia spacing on the transport of mucus. A summary of the results and the future perspectives conclude this manuscript in Section 4.

Numerical method

From the numerical point of view, modeling the complete problem under realistic conditions remains a huge challenge and several assumptions need to be made to simplify the problem, while keeping the essential ingredients of the physical mechanisms in the model. In particular, it is assumed that there is no mass transfer at the wall, that both fluids are Newtonian and that the cilia are equally spaced filament-like structures. More assumptions are detailed in this section, together with the numerical setup and geometry.

Geometrical modeling and beating pattern

The computational domain is a box with a regular mesh composed of N x × N y × N z points. The cilia are equally spaced along the bottom (x-y) wall, such that their base points are located at z = 0 (see figure 1(a) for a schematic view of the geometry). Their motion is imposed to be in the x-direction only. The spacing between two adjacent cilia is denoted a in the x and y-directions. The wavelength λ of the MCW is set such as λ = N x when the metachrony is forced. For the case of synchronously beating cilia, N x is chosen to be 8a. The length of the cilia L is set to 15 lattice units (lu) except when stated differently; and 200 snapshots in time per beating cycle are uniformly distributed to model their motion. The ratio h/H between the PCL thickness and the height of the domain is fixed to 0.27 for all simulations.

The equations of motion for the cilia are inspired from [START_REF] Chatelin | Méthodes numériques pour l'écoulement de Stokes 3D : fluides à viscosité variable en géométrie complexe mobile[END_REF] and reproduce the beating pattern by resolving a 1D transport equation along a parametric curve. Let P (ζ, t) be the position of the curve at time t and at a normalized distance ζ from the base point of a cilium. With appropriate boundary conditions, a realistic beating pattern can be obtained using the following transport equation:

∂P ′ ∂t + ν(t) ∂P ′ ∂ζ = 0 BC: P (0, t) = (0, 0, 0) P ′ (0, t) = (2 cos(2πt/T ), 0, cos(2πt/T )) (2.1)
with ν(t) = [1 + 8 cos 2 (π(t + 0.25T )/T )]/T being the viscosity of the surrounding fluid, T the beating period, and P ′ = ∂ ζ P . The resulting angular amplitude between the beginning and the end of a stroke phase is θ = 2π 3 , which agrees well with experimental data [START_REF] Sleigh | The propulsion of mucus by cilia[END_REF]. Note that a 3D beating pattern would allow to achieve more realistic simulations while being more computationally expensive. Hence, the choice has been made to use this 2D beating pattern. It captures the essential ingredients of the beating, and as cilia have a diameter smaller than a lattice unit, the difference of induced flow between two cilia overlapping in 2D or slipping onto each other in 3D is very small. Figure 1(b) gives a view of the beating pattern obtained by resolving Eq. (2.1). Note that with the present model both phases take the same amount of time while moving in the same (x,y) plane. This choice has been made in order to only study spatially asymmetric motions, which are the only mechanisms effective at low Reynolds numbers. The temporal asymmetry is indeed a mechanism that can enhance the flow only when inertial forces are no longer negligible [START_REF] Khaderi | Breaking of symmetry in microfluidic propulsion driven by artificial cilia[END_REF]. Nevertheless, when the feedback of the fluid is taken into account (as in §3.1), a non symmetrical motion will develop, with a stroke phase slower than the recovery phase. More details regarding this temporal asymmetry are given in §2.3.

The PCL is set such that it fills the region going from the bottom of the domain (z = 0) up to an altitude of h. In all the simulations, the value of h = 0.9L has been used in order to allow the tips of the cilia to emerge into the mucus layer during their stroke phase as observed in real epithelium configurations.

Both PCL and mucus are considered to be Newtonian fluids. The kinematic viscosity of the mucus is ν m = 10 -3 m 2 /s, and the viscosity ratio r ν = ν m /ν P CL between the mucus and PCL is set to 20. It has indeed been recently shown [START_REF] Chatelin | A parametric study of mucociliary transport by numerical simulations of 3D non-homogeneous mucus[END_REF] that mucus transport was maximized for viscosity ratios ranging from 10 to 20 with a stiff transition between the two fluid layers. The beating period of the cilia is equal to N it × dt (with dt = 1 using the classical LBM normalization), N it being the number of iterations for a cilium to perform a complete beating cycle. An oscillatory Reynolds number Re osc , based on the velocity of the cilia's tips U cil = 2θL

Tosc can now be defined:

Re osc = U cil L ν mucus = 4π 3 Tosc L 2 ν mucus = ωL 2 ν mucus (2.2)
with ω being the angular frequency of the cilia beating. With realistic physical quantities corresponding to the ciliated epithelium surface, the value of Re osc is of the order of 10 -5 (using L ≈ 10 -5 m, ν mucus ≈ 10 -3 m 2 /s and U cil ≈ 10 -3 m/s).

To avoid running simulations at such a low Reynolds number which would require a very high number of iterations using a lattice Boltzmann scheme, a higher Reynolds number was chosen: Re = 20. Indeed, achieving simulations at Re < 1 using a LB scheme, while still describing well the cilia (L = 15 lu at least), requires a lot of CPU resources. Thus, inertial effects are introduced in the model but it has been carefully checked that the results remain the same for creeping flows (Re = O(10 -2 )) by comparing them to the results obtained with LBM formulation designed for Stokes flows [START_REF] Guo | Lattice Boltzmann Method and Its Applications in Engineering[END_REF][START_REF] Zou | A improved incompressible lattice Boltzmann model for time-independent flows[END_REF]. The differences are found to be less than 7% on the transported velocity for all phase lags ∆Φ = 0. The corresponding beating patterns of the cilia shown in the following are similar also in terms of vorticity generation and the same qualitative coordination concerning the antipleptic/symplectic behaviour is observed. However, for the particular case of synchronously beating cilia (i.e. ∆Φ = 0), the inertial effects will play a non-negligible role in the flow dynamics. Despite being weak, they will indeed cancel the reversal of the flow that should occur when the cilia are in the recovery phase. More details regarding the inertial effects will be given in §3.2.7.

Finally, it is worth noticing that the lattice has been chosen such that the numerical diffuse diameter of the cilia due to IB method corresponds to real cilia diameter (≈ 0.3µm). Hence, a realistic drag is taken into account in the present study.

Algorithm

The numerical model is described in [START_REF] Li | An improved explicit immersed boundary method to couple with lattice Boltzmann model for single-and multi-component fluid flows[END_REF], and validated on several configurations involving flexible and moving boundaries in multiphase flows, with a 2 nd order accuracy. Briefly, the idea is to add a forcing term

F σ i = F SC σ + F IB σ to the discrete LB equation for each σ fluid component. F SC σ
is an interparticle potential force that takes into account the fluid-fluid cohesion forces [START_REF] Shan | Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation[END_REF], and F IB σ is an IB-related force for ensuring the no-slip condition at the fluid-solid interface.

The fluid part is first solved on a Cartesian grid with LBM using the Bhatnagar-Gross-Krook (BGK) operator and a D3Q19 scheme. The collision and streaming steps proper to the LBM method are first performed. The model of [START_REF] Porter | Multicomponent interparticle-potential lattice Boltzmann model for fluids with large viscosity ratios[END_REF] is used to model the two-phase flow as it allows to minimize the magnitude of spurious currents near the fluid-fluid interface. More importantly, it also allows to consider higher density or viscosity ratios. Then, values for the fluid velocity are interpolated at the Lagrangian points. It allows to compute an IB force to be spread onto the neighbouring Eulerian fluid nodes in order to ensure the no-slip condition along the cilia. The macroscopic fluid velocity is then updated. The motion of each cilium is decomposed into a finite number of steps (snapshots) during a period. If necessary, an interpolation can be done in order to have the velocity values along the cilia in between two steps. Note that the geometric shape of the beating is fixed in all simulations and is not impacted by the feedback law introduced in §2.3, which only affects the duration of the recovery and stroke phases.

Since the model of [START_REF] Porter | Multicomponent interparticle-potential lattice Boltzmann model for fluids with large viscosity ratios[END_REF] uses a Shan-Chen (SC) repulsive force [START_REF] Shan | Lattice boltzmann model for simulating flows with multiple phases and components[END_REF][START_REF] Shan | Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation[END_REF], surface tension effects emerge intrinsically at the PCL-mucus interface. Hence, a sharp interface between the mucus and the PCL can be maintained at any time. However, small diffusion effects might occur on the lattices bordering the interface. Additionally, the cilia which enter the PCL-mucus interface may punctually induce a small mixing. This is however corrected by the SC force which "unmixes" the two fluids.

Periodic boundary conditions are used in the x and y directions, while no-slip and free-slip boundary conditions are used at the bottom and top walls respectively. The size of the computational domain ranges from 50 lu to 400 lu depending on the configuration considered, except for the size along the z direction which is always set to 50 lu.

Taking advantage of the local character of the LBM algorithm, the code is parallelized using MPI libraries (Message Passing Interface), by splitting the full computational domain into 9 sub-domains of size (N x /3, N y /3, N z ). More details on the numerical model can be found in [START_REF] Li | An improved explicit immersed boundary method to couple with lattice Boltzmann model for single-and multi-component fluid flows[END_REF].

Feedback of the fluids onto the cilia

The basic idea is to modulate the beating motion of the cilia as a function of the fluids motion. To do so, it is assumed that all cilia follow the same beating pattern, meanwhile a feedback of the fluids, which consists in accelerating or slowing down the motion of the cilia, is introduced.

Each cilium is discretized with N s = 20 Lagrangian points. Let s be the subscript corresponding to the s th Lagrangian point, starting from the base tip at s = 1, and V s i the velocity on the s th Lagrangian point of the i th cilium. For each cilium, we define the average velocity over all Lagrangian points V i , which is linked to the number of steps (snapshots) this cilium will skip during one iteration of the fluid solver. The fluid feedback onto the cilia thus consists in modifying the norm of the velocity vector V i , while its direction remains unchanged.

The feedback is computed in three steps. First, the IB forces corresponding to mucus and PCL are projected onto the corresponding velocity vectors for each Lagrangian point. Then, an estimate of the feedback is computed based on the torques of the forces for each Lagrangian point. Finally, the beating pattern of the cilia is adjusted at the beginning of the next time step. The different forces and geometrical variables used are illustrated in figure 2, and the three steps explained below.

(i) The interpolated IB forces applied by the i th cilium onto the fluids -respectively F i m for the force imposed on the mucus phase and F i P CL for the force imposed on the PCL phase-, are projected on V s i :

F i,proj * = F i * • V s i V s i 2 V s i , (2.3)
where * stands for "m" or "PCL" depending on the position of the Lagrangian node. In order to take into account the difference of viscosity between the two layers, the forces F i,proj m and F i,proj P CL are weighted by a term under the form:

τ * τm+τP CL . The total projected force of the fluids onto the s th segments of the i th cilium F i f luids→cilia writes:

F i f luids→cilia = - (τ m F i m + τ P CL F i P CL ) • V s i (τ m + τ P CL ) V s i 2 V s i (2.4) (ii)
For each Lagrangian point s, the norm of the torque of F i f luids→cilia with respect to the base point O of the i th cilium is computed by:

M i O (F i f luids→cilia ) = F i f luids→cilia L p (2.5) with L p = Xp⊗V s i V s i
being the lever arm. The total fluid feedback onto the i th cilium considered is then computed by summing the computed quantities over all Lagrangian points:

T i = Nsegments s=2 F i f luids→cilia L p (2.6)
The velocity of each cilium is finally modified as follows:

V i = V 0 + αT i (2.7)
where V 0 is the initial speed of the cilia, and α a coupling parameter allowing the tuning of the feedback strength.

(iii) Then, at the beginning of the next iteration, the beating pattern of the i th cilium is adjusted:

N t i = mod (N t-1 i + V i , N total i ) (2.8)
where N total i is the total number of snapshots defining the beating pattern of the cilia, N t-1 i the previous snapshot, and N t i the new position of the i th cilium (in terms of number of snapshots) at the current iteration.

Injection of passive tracers

In the context of real epithelial systems in human body, mucus acts as a barrier against particles and pollutants. To gain insight on how particles are dispersed and advected in the mucus and PCL layers, passive tracers are injected into the domain on a (y,z) plane, when the flow has reached an established regime (see figure 3 for a view of a domain filled with tracers). Their displacements are then computed and averaged over several cilia beating cycles. Their equations of motion are solved by a second-order Runge-Kutta (RK2) scheme, using the interpolated fluid velocity at each time step and the same procedure as in the IB method.

By taking a closer look at figure 3, one can observe that particles initially seeded into the PCL are greatly mixed while staying in the PCL. On the contrary, particles initially seeded into the mucus layer stay in the mucus layer and are not mixed. This is mainly due to the surface tension effects present at the mucus-PCL interface which prevent a too strong mixing at the interface. This shows how wetting particles that are deposited to the air-mucus interface might enter the mucus layer but never reach the PCL one.

Results & Discussion

Emergence of metachronal waves

Using the feedback force introduced in §2.3, randomly beating cilia are observed to synchronize with their immediate neighbours giving birth to symplectic or antipleptic metachrony. The time for the synchronization to occur depends on the set of parameters used. Both local and global synchronizations can be observed. The empirical parameter α plays a role in the emergence of the waves. When α is set to a too low value (|α| < 0.5), the cilia will beat randomly, while if set to a too high value (|α| > 7 for the case a/L = 1.67 for example), the cilia will fully synchronize with each other without any phase lag. Additionally, higher absolute values of α usually decrease the times needed for metachronal waves to emerge.

Figure 4 shows a symplectic MCW emerging from an initially random state (every cilium were initially beating at a random step of their beat cycle). In the present case, 3 rows of 8 cilia on a computational domain of size (N x = 361, N y = 136, N z = 50) are considered. The PCL is set such as h = 0.9L, and the ratio of viscosity between the PCL and the mucus phase is 15. The feedback coefficient α used is α = -3.5. The spacing between two neighbouring cilia is a/L = 3 in the x and y directions with L = 15 lu. Hence cilia never collide and overlapping of kernels from neighbouring cilia never occurs. One can clearly see the symplectic MCW that emerged from the initially random state of the cilia, with a wavelength λ = 180 lu, and a phase lag ∆Φ ≈ -π/2, confirming that hydrodynamic interactions only suffice to account for the emergence of MCW.

Figure 5 shows a similar configuration (h = 0.60L, r ν = 15, α = -3.5) where 1024 cilia arranged in a 32 × 32 square are considered on a computational domain of size (N x = 161, N y = 161, N z = 32). One can clearly see the formation of an antipleptic MCW with a wavelength λ = 80 lu, and a phase lag ∆Φ ≈ π/8. Notice that in the simulation presented in figure 5 the cilia spacing has been reduced (a/L = 0.23 and L = 22 lu), to take into account a larger number of cilia. Thus, neighbouring cilia may overlap, which can cause spurious numerical effects, due to the 2D beating pattern of the cilia, when two Lagrangian points with opposite velocities occupy the same mesh point. However, it has been verified on some simulations, that these effects do not play any role in the emergence of the MCW by setting to zero the IB forces contribution of these Lagrangian points, therefore canceling the spurious effects that may occur at the overlapping points of neighbouring cilia. In other words, the fluid velocity at the Eulerian node surrounding these particular Lagrangian points was not modified by the IB method, proving that it is not the source of the cilia synchronization. As shown in figures 4 and 5, weakly 3D effects may be observed but do not play a key role in the physics of emerging waves. In all the simulations presented, the synchronization along the x-direction is quite strong. In some simulations, particular cilia may lose synchronization over time, but quickly readjust their beating accordingly to the others. It also appears that once an equilibrium is reached, the synchronization along the y-direction is stable too. A remarkable fact is that for particular configurations (as the one displayed in figure 5), the interface between the mucus and the PCL also forms a wave traveling in the same direction as the MCW. If the mucus viscosity is increased, the torques felt by the cilia become stronger. Hence, the cilia are either more accelerated (α > 0) or decelerated (α < 0), and the mucus will adapt its displacement accordingly. However, when the stationary regime is reached, the clearance velocity should not be strongly modified as shown by [START_REF] Chatelin | A parametric study of mucociliary transport by numerical simulations of 3D non-homogeneous mucus[END_REF]. According to these authors, the mucus velocity is decreased by only 50% when the viscosity ratio is increased by a factor of 100000. Preliminary results (not shown) seem to indicate that, assuming a least effort behaviour of the cilia (meaning α < 0), antipleptic waves are obtained for small cilia spacings (a/L 1.5) while symplectic waves are seen to emerge for higher cilia spacings (a/L 1.5). Since in nature, antipleptic waves are often observed for densely packed cilia, it implies that natural cilia adopt a least effort behaviour. It also suggests the existence of a critical value for the cilia density, ρ c = ρ(α c ) where α c is the corresponding critical value of the coupling parameter, from which the kind of waves emerging can be controlled. On the contrary, with the present model, assuming that the cilia beat faster when encountering a resistance (meaning α > 0), symplectic waves are seen to emerge for small cilia spacings (a/L 1.5) while antipleptic waves emerge for higher cilia spacings (a/L 1.5), which is not observed in nature. Movies of the MCW are provided as online supplementary materials.

As already explained in §2.1, when the feedback is taken into account, a non symmetrical motion develops, with a stroke phase slower than the recovery phase for both antipleptic and symplectic MCW. For the antipleptic case, V i is smaller during the stroke phase, and since α is negative and the torques T i are always positive, it means that the values of T i computed during the stroke phase are larger than the values computed during the recovery phase. It results in a weaker velocity for the cilia which cover less snapshots during the stroke phase. It is the opposite for the symplectic case where α is positive: T i takes larger values during the recovery phase as V i . In other terms, the feedback depends on the clustering character of the motion: when cilia are clustered, the torque exerted by the fluids onto the cilia is weaker, and when they are far from each other (during the stroke phase for antipleptic motion and during the recovery phase for symplectic motion), the torques are stronger. This makes sense since the cilia encounter more viscous resistance when they are not clustered, as it will be discussed further down in §3.2.4. The resulting motion of the cilia then differs from what is observed in nature, but it indicates that cilia experience stronger stresses during the stroke phase of antipleptic motion and during the recovery phase of symplectic motion. It supposes that the beating kinematics of the real cilia are not dictated only by hydrodynamical interactions and suggests that other biological parameters or functions (such as sensing) may play a role. This is in agreement with [START_REF] Guo | Cilia beating patterns are not hydrodynamically optimal[END_REF] who compared the performance of pumping-specialized cilia and swimming-specialized cilia as a function of the metachronal coordination, and found that the later almost always outperforms the pumping-specialized cilia. As it will be further detailed in §3.2.5, their results are also in accordance with the outcome of the present paper: antipleptic waves are the most efficient ones to transport fluids. Finally, it is worth noticing that the degree of asymmetry is much higher for antipleptic MCW compared to symplectic MCW.

Quantitative study of the metachronal waves

This section presents the results obtained, once the flow is well-established, for the three configurations studied: (a) synchronized case (all cilia beat together with no phase lag); (b) symplectic MCW where two neighbouring cilia beat with a negative phase lag, i.e -π < ∆Φ < 0; and (c) antipleptic MCW where two neighbouring cilia beat with a positive phase lag, i.e 0 < ∆Φ < π. Note that in all following results, contrary to §3.1, the metachrony is imposed in order to study specific phase lags ∆Φ. Hence, the size of the domains must be changed accordingly, since different phase lags ∆Φ imply different wavelengths, and hence more or less cilia. However, the quantities presented here have been averaged over the domain, and there are no effects in changing the size of the box. In order to study only spatially asymmetric motions, the recovery and stroke phases now take the same amount of time in all following results. Also note that the standard deviations are not displayed in any of the figures that will be presented as they are extremely small (less than 0.01%) for all considered quantities. Hence they do not give additional information. Schematic views of synchronized beating and metachronal motions are displayed in figure 6. Another configuration, corresponding to randomly beating cilia is also represented in figure 6. One can observe that the synchronized motion of cilia creates vorticity only in the periciliary zone, whereas metachronal motions also induce vorticity into the mucus layer. It is worth noticing that for symplectic motion (case (b)), vortex trails emerge from the cilia tips performing their stroke phase. The same phenomenon is observed in the antipleptic case during the recovery motion of cilia.

The presence of coherent vortices is clearly visible in both cases.

Transport and mixing zone

The displacement field is calculated such as d(x)= T 0 u(x(t), t)dt, where x is the position vector, and u the fluid velocity. Its component over the x-direction as a function of the axial position is averaged over several beating cycles and displayed in figure 7. The transport and mixing zones defined in [START_REF] Ding | Mixing and transport by ciliary carpets: a numerical study[END_REF] are clearly visible. It shows that the antipleptic MCW perform better than the other kind of coordination to transport the particles. Note that in figure 7b, the two curves with ∆Φ = ±π/4 have a similar shape, whereas in figure 7a symplectic MCW induce a negative transport (i.e. a counter-flow) in the mixing zone for the same phase lag. For a given value of phase lag, the transport depends on the density of cilia. One can also notice the importance of the phase lag by looking at figure 7a: for a phase lag ∆Φ = -π/9, the symplectic MCW induce no 7b), the present results show that the displacement induced by the symplectic wave has the same shape as the synchronized and antipleptic cases, i.e it induces transport in the mucus phase even if it is smaller. On the contrary, [START_REF] Ding | Mixing and transport by ciliary carpets: a numerical study[END_REF] obtained no transport at all for the symplectic wave in the mucus phase, but instead a peak of transport under the cilia tips in the PCL at z = 0.7L due to the presence of a vortex-like structure both below and above the cilia tips for this particular value of the phase lag ∆Φ. Since the main difference between both studies is the beating pattern of the cilia (same phase lag and cilia spacing), it highlights the sensitivity of the system to the beating pattern.

The displacement in the z and y-directions have also been analyzed. The results show that the y-component of the displacement could be neglected and that the z-component of the displacement was almost null above the cilia tips, but not zero under. As it will be discussed later, the z-component of the displacement reaches its maximal value in regions where the shear rate is maximal too. It agrees with the existence of a mixing zone under the cilia tips.

By doing an analogy with the strain-rate tensor classically used in solid mechanics, the gradient of the displacement field ∇d can be computed by considering each Eulerian node of the Cartesian grid as passive tracers. Since the displacement in the transverse direction is very small, the y-component of the displacement field can be neglected. Henceforth, the gradient is computed over every (x,z) plane, and an average is then done over the y-direction. Following the methodology described in [START_REF] Ding | Mixing and transport by ciliary carpets: a numerical study[END_REF], ∇d is decomposed into an antisymmetric component R = (∇d -(∇d) T )/2 corresponding to rotation, and a symmetric component S = (∇d + (∇d) T )/2 corresponding to shear deformation. The two eigenvalues of S are of the form ±γ and indicate the rates of stretching (+γ) and compression (-γ). The unit eigenvector e γ corresponding to the positive eigenvalue indicates the direction of stretching, and the other eigenvalue the direction of compression. Because of the incompressibility condition, both eigenvectors are orthogonal, and so plotting only one of them is sufficient to have the complete set of information. In figure 8, the stretching rate and its direction are presented in a (x,z) plane for the three cases studied. The stretching rate is maximal near the upper part, and at little distance above the cilia tips in all three cases. For the synchronized motion (case (a)) there is almost no stretching away from the cilia during the recovery phase whereas for cases (b) and (c), corresponding respectively to symplectic and antipleptic motion, a weak stretching can always be observed (see the right side of cases (b) and (c) in figure 8). A complex shape of the stretching rate is observed at the cilia tips during the stroke phase in antipleptic motion (see cilium 3 for instance); during the stroke phase in symplectic motion (see cilia 6 and 7) and during the stroke phase in synchronized motion (results not shown). From their orientations, one can expect an enhancement in the mixing in this region. As in [START_REF] Ding | Mixing and transport by ciliary carpets: a numerical study[END_REF], the stretching direction is a nonlinear function of space. Except for zones where the shear rate is maximal, the γ e γ field is then oriented at 45 • compared to the x-direction, and almost uniform. This is reminiscent of a linear shear profile d x = cz where c is a constant. The nonuniform aspect of the stretching orientation indicates the presence of "folding" in the displacement field d which also plays a role in the mixing as explained by [START_REF] Kelley | Separating stretching from folding in fluid mixing[END_REF]. [START_REF] Ding | Mixing and transport by ciliary carpets: a numerical study[END_REF] obtained that the stretching rate was maximal for the antipleptic case for ∆Φ = π/4 which is exactly what is obtained here, even if the beating patterns used are different. One can then expect that this value of phase lag ∆Φ is the most efficient in mixing fluids by enhancing the stretching near the upper part of the cilia.

In order to assess the reliability and robustness of the solver, a comparison with experimental data is also performed. The average clearance velocity was computed on the plane (x,y,3.2L) for all simulations. The highest clearance velocity is reached for the case a/L = 1.67 with ∆Φ = π/4 and equals 33.47 µm/s. It is in good agreement with the experimental results of [START_REF] Matsui | Coordinated clearance of periciliary liquid and mucus from airway surfaces[END_REF] who observed a clearance velocity of 39.8 ± 4.2 µm/s for the mucus. [START_REF] Matsui | Coordinated clearance of periciliary liquid and mucus from airway surfaces[END_REF] also observed that the PCL flows in the same direction as the mucus, which is the case in the present study. Finally, the PCL-mucus interface remains approximately flat in most of the simulations presented, as it is the case in the micrographs performed by [START_REF] Sanderson | Ciliary activity of cultured rabbit tracheal epithelium: Beat pattern and metachrony[END_REF] on rabbit tracheal epithelium.

Directional pushing efficiency

Inspired by the works of [START_REF] Khaderi | Microfluidic propulsion by the metachronal beating of magnetic artificial cilia: a numerical analysis[END_REF][START_REF] Gauger | Fluid transport at low Reynolds number with magnetically actuated artificial cilia[END_REF], and [START_REF] Kim | Pumping fluids with periodically beating grafted elastic filaments[END_REF], a positive flux Q p and a negative flux Q n are defined as follows: the x-component of the velocity field is considered over the (N x ,y,z) plane and, at each time step of a full beating cycle, the negative and positive velocity values are separated. The instantaneous negative and positive fluxes are then computed over a sufficiently large number of cycles. The difference (Q p -Q n ) gives the net instantaneous flux; and the directional pushing efficiency ǫ P N is defined as :

ǫ P N = (Q p -Q n )/(Q p + Q n ).
Note that here, contrary to the previous works of the aforementioned authors, the domain is not restricted to the transport area but covers the whole flow region instead. Such a choice is justified by the fact that it has been experimentally observed, using confocal microscopy and fluorescent markers, that PCL and mucus are in reality transported at approximately the same rate [START_REF] Matsui | Coordinated clearance of periciliary liquid and mucus from airway surfaces[END_REF]. As a matter of fact, the PCL transport seems to depend on the presence of a mucus layer above it, and ciliary mixing is thought to be responsible for the diffusion of momentum from mucus to PCL. Taking into account the transport zone and the mixing zone in the computation of ǫ P N then allows one to collect information about how transport and mixing work together during a beating cycle.

The variation of ǫ P N over one beating cycle for five different phase lags ∆Φ, all for antipleptic MCW with a cilia spacing such as a/L = 2, is illustrated in figure 9. Let us remind that, although the quantities considered here are the results of the motion of all cilia, the calculations have been made through the beating cycle of a single cilium of reference. Therefore, the fact that the corresponding curves have a huge decrease around t = π is related to the choice of this particular cilium. Indeed, the plane used to compute the positive and negative fluxes is located close to a particular cilium. This cilium is in the recovery phase at t = π, a small value of ǫ P N at this particular time. Choosing another cilium or another plane would shift horizontally the drop observed around t = π.

To interpret the evolution of ǫ P N in figure 9, it is necessary to examine the topology of the flow during this beating cycle, shown in figure 10. As cilium 1 begins its recovery phase (case (a) of figure 10), there is no reversal flow at the periodic boundary frontier. Hence, ǫ P N = 1, as it is visible in figure 9 at the point (a). Then, when cilium 1 carries on its recovery phase (case (b) of figure 10), no reversal flow can be observed since cilium 1 is still at the beginning of its recovery phase and cilium 8 is finishing its stroke phase. At t = π/2 (case (c) of figure 10), both cilia 1 and 8 are in their early recovery phase and the directional pushing efficiency begins to decrease (see points (b) and (c) on figure 9). At t = 3π/4 (case (d) of figure 10), cilium 1 is reaching the maximal speed of its recovery phase. The negative flow generated is important, and the efficiency quickly drops: see point (d) on figure 9. Between t = 3π/4 and t = π, the directional efficiency weakly increases. This is due to the positive velocity near the base of cilium 1 which is doing a "whip-like" motion: while its free end is still going to the left, its base starts moving on the right. Case (e) of figure 10 corresponds to the end of the recovery phase of the first cilium and at the same moment, cilium 8 is at its maximal negative velocity: henceforth the directional pushing efficiency is still decreasing. At t = 5π/4 (case (f) of figure 10), cilium 8 has finished its recovery phase and cilium 1 is in the middle of its stroke phase. It counteracts the reverse flow created earlier and ǫ P N finally increases. At t = 3π/2 and t = 7π/4 (cases (g) and (h) of figure 10), both cilia 1 and 8 are in their stroke phases, and the flow generated is purely in the positive x-direction: see points (g) and (h) on figure 9.

It is now interesting to compare the average directional pushing efficiency < ǫ P N > of the different kinds of synchronization over a full beating cycle (see figures 11 (a) and (b)). Before going further into details, let us recall the fact that the current efficiency < ǫ P N > is a criterium qualifying the non-isotropy of the transport. Thus, it only gives an insight of the capacity of the cilia to transport the flow in an unidirectional way, but does not give any quantitative information on the volume of fluid flow, which is effectively displaced. Figures 11 (a) and (b) show an interesting phenomenon: systems with larger cilia spacing have a better capacity to transport the flow in the same direction. Another intriguing fact is that, for the smallest cilia spacing (a/L = 1.67, figure 11), antipleptic MCW seem to have a better ability to create an unidirectional flow compared to synchronized or symplectic motion, with a peak for the average efficiency around π/2. It agrees particularly well with the results of [START_REF] Gauger | Fluid transport at low Reynolds number with magnetically actuated artificial cilia[END_REF] who reported that antipleptic MCW were more efficient than the synchronized motion of cilia, itself being more efficient that the symplectic case. It also partially agrees with the results of [START_REF] Ding | Mixing and transport by ciliary carpets: a numerical study[END_REF] who obtained peaks of efficiency for phase lags around ±π/2, with a stronger one for the antipleptic case. In the present study, the minimal efficiency for metachronal motion is reached for a symplectic MCW with a phase lag of approximately -π/2.

When the cilia spacing is increased to a/L = 2 (figure 11), a peak of efficiency appears for the antipleptic MCW for ∆Φ = 2π/3; and the lower values are reached in the symplectic case for ∆Φ ≈ -π/2 and ∆Φ ≈ -π/3. These trends will be found again in the results shown in §3.2.3. Then, if the cilia spacing is increased above a limit value (a/L = 3.33, figure 11), the directional pushing efficiency becomes equal to 1 for all cases (antipleptic, symplectic, and synchronized). It agrees well with the results obtained by [START_REF] Khaderi | Microfluidic propulsion by the metachronal beating of magnetic artificial cilia: a numerical analysis[END_REF] who concluded that "the amount of flow enhancement depends on the inter-cilia spacing [but] the efficiency is not significantly influenced". Indeed, as the cilia are set away from each other, the influence of their neighbours become negligible. Nevertheless, the authors reported a very low efficiency for synchronized cilia, whereas the present results always show a positive flow in the mucus phase for synchronized motions. As mentioned in §2.1, this is a direct consequence of the inertial effects (Re > 1) in the present simulations. It is recalled that they only affect the synchronized case. It has carefully been checked that all other conclusions drawn for antipleptic and symplectic metachrony remain the same for Re < 1 and Re = 20. For a detailed study of the inertial effects at Re = 20, see §3.2.7. The present results are different from the results of [START_REF] Gueron | Cilia internal mechanism and metachronal coordination as the result of hydrodynamical coupling[END_REF], as it appears that even when cilia spacing is higher than 2 cilium lengths, the influence of neighbouring cilia cannot be neglected (see the case a/L = 2.53 for example).

It is important to remember that this efficiency does not take into account the actual net flow volume transported.

Compared hydrodynamic efficiency

After having investigated the capacity to transport the flow in the desired direction, it is now interesting to look at the actual flow volume displaced to see if a better efficiency to transport the flow directionally corresponds to a better flow transport. To do so, the global volumetric flow rate Q v over a unit volume of size (1 × 1 × N z ) is introduced:

Q v = N z U * dx 2 L 2 (3.1)
where dx = 1 using the classical LBM normalization, and U * = U av /U ref , with U ref =(λ/N cil )/T being the reference velocity and U av = 1 ninj n k i,j,k U ijk the average fluid velocity inside the whole domain.

On figure 12, the dimensionless volumetric flow rate Q v is plotted over one beating cycle for arrays of cilia with different phase lags (∆Φ = 0, ∆Φ = -π/4, and ∆Φ = π/4) and different cilia spacings. For cilia spacings equal to a/L = 1.67 and a/L = 2, the antipleptic MCW is the most efficient for transporting fluids. However, note that as the cilia spacing is increased, the ability of the antipleptic wave to transport fluids exhibits a huge decrease compared to the symplectic wave. For a larger cilia spacing (a/L = 3.33), the fluid transport is not impacted anymore by metachrony. Finally, no flow reversal occurs for the synchronized cases: this is the consequence of working at Re = 20.

The total volume of fluid displaced during a beating cycle for the different phase lags is compared on figures 13 (a) and (b). For a small cilia spacing (a/L = 1.67), the efficiency of the antipleptic metachrony is obvious. It agrees well with the results of [START_REF] Khaderi | Microfluidic propulsion by the metachronal beating of magnetic artificial cilia: a numerical analysis[END_REF] who observed a larger net flow produced by antipleptic metachrony for this value of cilia spacing. Symplectic waves appear to be less or at best equally efficient than antipleptic motion, except for ∆Φ = -7π/8 for a/L = 1.67 where there is a peak in the total displaced volume of flow. On the contrary, negative peaks are found for a/L = 1.67 and a/L = 2 for ∆Φ = -π/3. There are two neighbouring maxima at ∆Φ = π/4 and ∆Φ = π/2 for a/L = 1.67 and a/L = 2 respectively, indicating that specific phase lags are more able to generate a strong flow. At this point, it is worth remembering that in the present model, the recovery and stroke phases occur in the same plane which is not the case in real configurations. Accordingly to Downton & Stark (2009), who studied both 2D and 3D beating patterns, it is expected that in presence of 3D beating patterns, both the directional pushing efficiency ǫ P N and total displaced volume of flow would increase. It is also expected that a small flow might occur in the y-direction.

Energetic cost of the metachronal waves

The previous sections have been dedicated to the study of the directional transport and the volume displacements of fluids. An energetic perspective is now introduced, to characterize completely the potential benefits of metachrony.

The average power spent by the cilia during a beating cycle is given by:

P cil = s,i V s i • (F i m + F i P CL ) N cilia (3.2)
using the forces illustrated in figure 2. The power spent is averaged over several beating cycles. To have a dimensionless power P * , the power P ∞ spent by an isolated cilium (a/L = 10) is computed such that P * =P cil /P ∞ . Before going any further, it is important to remember here that the only parameter that differs between each value of phase lag ∆Φ for a given cilia spacing a/L, is the size of the domain N x over the x-direction, hence the number of cilia acting in one wavelength. Moreover, between different cilia spacings a/L, the space a between two cilia is modified in both the x and y-directions, and the sizes N x and N y of the domain must be changed accordingly. As a consequence, if the cilia spacing is increased, the density of cilia in both the x and y-directions is decreased. In all the simulations presented in §3.2, all the other parameters (ratio of viscosity r ν , beating frequency f , length of the cilia L, etc.) are fixed.

Figures 14 (a) and (b) show the average dimensionless power as a function of ∆Φ. For small values of phase lags (|∆Φ| π/2), the average power spent decreases to a minimal value for ∆Φ ≈ π/4. For this particular value, the system spends less power than the synchronized case: the antipleptic MCW allows the system composed of all cilia to encounter a smaller resistance from the surrounding fluids.

Nevertheless, it is obvious, from a fluid mechanics point of view, that when cilia beat together in a synchronized way, the viscous resistance felt by each cilium is reduced.

In the case of metachrony, the behaviour is dual. Indeed, cilia in the stroke phase of the antipleptic motion encounter a stronger viscous resistance from the flow due to their added respective remoteness. Therefore, compared to the synchronized motion, they will transfer more energy to the flow; and the energy transferred will be entirely used to propel the fluids in the desired direction. The opposite is true during the recovery stroke of the antipleptic motion: the clustered behaviour of the cilia allows them to experience a lower viscous resistance from the flow, and therefore to limit the amount of power transferred during this phase. But, globally, the system with antipleptic MCW and |∆Φ| < π/2 encounter less viscous resistance than the synchronized motion of cilia, as it can be seen on figures 14 (a) and (b). It results in a better efficiency of the cilia with antipleptic metachronal motion to transfer their momentum to the flow, while in the meantime requiring less power. For the symplectic motion, a similar phenomenon happens: cilia in their stroke phase are clustered and therefore unable to fully exert their pushing action, while cilia in their recovery phase are away from each other and generate a stronger reversal flow compared to the antipleptic motion. It results in a lower capacity of the symplectic MCW to transport mucus, while still allowing the system to spend less power than synchronously beating cilia for |∆Φ| < π/2. An energetic ratio will be introduced in §3.2.5 to quantify, for such value of ∆Φ, the capacity of the cilia with antipleptic or symplectic metachrony to transfer their momentum to the flow. This is not true for large phase lags (|∆Φ| π/2): the system with antipleptic or symplectic metachrony spends more power than the synchronous case. One can suppose that it is a direct consequence of the reduced number of cilia in a wave length. Then, to explain the better efficiency of the antipleptic MCW over the symplectic ones for large phase lags, an investigation of the flow topology is necessary. Figure 15 shows an antipleptic MCW (on the left) and a symplectic MCW (on the right) at the same time for a phase lag ∆Φ = ±2π/3. For such value of ∆Φ, the antipleptic wave is more efficient for transporting the mucus (see figures 13 (a) and (b)), although the average power spent by both system is relatively similar (see figures 14 (a) and (b)). One can easily see the main difference between both cases: in the symplectic case, the cilium in stroke phase encounters the reversal flow generated by the other cilium in recovery phase. When these two secondary flows meet, vortices are generated and the global transport of fluid is less efficient, as a fraction of the energy transferred to the flow is used to cancel this reversal flow. On the contrary, in the antipleptic case, the cilium in stroke phase does not immediately feel the influence of the reversal flow created by the other cilium in recovery phase, which is behind him. Henceforth, this cilium is able to fully exert its pushing effect on the mucus phase. Moreover, for the antipleptic case, there is a suction effect due to the combined motion of the cilia in the stroke and recovery phases, maximizing the propulsion of mucus. One can then expect a weak blowing effect between the cilia in the stroke and recovery phases of symplectic motion.

When the cilia spacing is large (see the case a/L = 3.33 of figures 14 (a) and (b) for example), the power spent by the system is the same for all phase lags: the cilia are too far from each other to be impacted by the phase difference of their immediate neighbours. One can observe that for ∆Φ = 0 (i.e. for synchronized beating), the average power spent by a cilium P * is equal to 1 in figure 14 (a) for the cilia spacings a/L = 1.67 and a/L = 3.33, meaning that synchronous systems spend the same amount of power as an isolated cilium. It shows the beneficial cost of antipleptic metachrony for small phase lags.

Displacement ratio

Inspired by the work of [START_REF] Kim | Pumping fluids with periodically beating grafted elastic filaments[END_REF], a displacement ratio, which can be seen as the transport efficiency of the waves, is introduced to quantify the capacity of a given system to transport particles, with respect to a given amount of power. In that context, η 1 is defined by the mean displacement over the x-direction during one beating cycle, divided by the mean power that a cilium had to spend during this beating cycle. Since the main purpose of mucociliary clearance is to transport mucus, and since experimental data [START_REF] Winters | Roles of hydration, sodium, and chloride in regulation of canine mucociliary transport system[END_REF] report that the total thickness in the vertical direction of the mucus layer is in the range [1.4L; 10L], values for the displacement were taken on an arbitrary plane z/L = 3.2 near the extremity of the domain. To obtain a value for the displacement, the instantaneous average fluid velocity over the x-direction is computed, and the resulting value is then multiplied by the period of a full beating cycle, giving the mean displacement < d x > over one period on the (x,y,3.2L) plane.

By dividing this mean displacement with appropriate quantities, a dimensionless expression of the displacement ratio is obtained:

η 1 = < d x > N cil λ P * (3.3)
In the synchronized case, i.e. ∆Φ = 0, λ is infinite and thus the size of the domain over the x direction was used and divided by the number of cilia.

In figures 16 (a) and (b), one can see, from an energetic point of view, the superiority of the antipleptic wave to transport mucus for small cilia spacings, which confirms the previous findings. More than that, for the smallest cilia spacing (a/L = 1.67), a clear peak of efficiency can be seen for antipleptic MCW with ∆Φ = π/4. If the cilia spacing is increased (a/L = 2), similar results are found. Nevertheless, for a/L = 2.53, a different behaviour occurs: the displacement ratio is found to be the worst for ∆Φ = ±π/2. Then, for cilia far away from each other, the displacement ratio η 1 remains constant for all phase lags (the cilia do not feel the influence of the others anymore).

Similar results regarding the efficiency of the antipleptic MCW are found in [START_REF] Osterman | Finding the ciliary beating pattern with optimal efficiency[END_REF], where the optimal beating pattern is investigated. In their study, they observed that antipleptic MCW were often the most efficient. They also found, as it is the case in the present study, that increasing the density of cilia results in an increase of the efficiency up to a critical point where clustering becomes counter-productive.

Mixing

To investigate how the mixing can be enhanced by metachrony, the average stretching rate over the transport and mixing areas during a beating period is computed. Figure 17 shows the results obtained for all cilia spacings. Clearly, the antipleptic MCW is the most efficient for stretching fluids, hence to mix fluids. Clear peaks are visible for antipleptic waves with ∆Φ = π/4 and a/L = 1.67; with ∆Φ ≈ π/3 and a/L = 2; and with ∆Φ ≈ π/3 and a/L = 2.53. On the contrary, symplectic MCW are almost always less efficient for mixing than antipleptic MCW, except for ∆Φ = -7π/8 and a/L = 1.67 where there is an enhancement in the mixing. This peak is also present for the opposite phase lag ∆Φ = 7π/8 where it reaches approximately the same value. This is perfectly coherent with all previous results, and partially with the ones of [START_REF] Ding | Mixing and transport by ciliary carpets: a numerical study[END_REF] who observed, for an unique layer of fluid with an uniform viscosity and a cilia spacing a/L = 1.67, the existence of two peaks in the shear rate: a weak one for symplectic waves for ∆Φ = -π/2, and a stronger one for antipleptic waves for ∆Φ = π/2. The main difference here is that in the present study, the symplectic MCW with ∆Φ ≈ -π/2 seem to have the worst capacity for mixing, while the antipleptic MCW reach their full mixing capacity for ∆Φ = π/4. The combined study of figures 14 and 17 gives a good insight of the mixing efficiency of the system.

Quantification of inertial effects

In this part, the influence of the Reynolds number is considered. The present objective is to compare the qualitative behaviour of the MCW between Re < 1 and Re > 1. To reach Reynolds numbers of the order of 10 -2 , the size of the cilia has been reduced to half the size of real cilia. The geometry is then divided by a factor 2 in the three spatial directions. The main consequence of such choice is that the computed quantities (fluxes, displaced volume of fluid, etc.) displayed here are approximately 8 times smaller compared to the results previously shown (as in figure 13). However, the qualitative behaviour of the MCW remains the same, as well as the drag exerted by the cilia. Thus, a comparison of the MCW behaviour between Re < 1 and Re > 1 is possible. On Figure 18, one can see the total displaced volume of fluid for Re = 0.02, Re = 1 and Re = 20, and a cilia spacing a/L = 1.67. A transition between Re < 1 and Re > 1 for the synchronous cases is expected, and can be seen on figure 18. As mentioned in §3.2.3, for Re 1, the synchronized cases are less efficient than the symplectic cases. Except for ∆Φ = 0, no other notable quantitative differences are present between the cases Re < 1 and Re = 20. The general behaviour of both the antipleptic and symplectic MCW remains similar for all the other phase lags ∆Φ. Especially, the presence of the peak around ∆Φ = π/4 is present for all the Reynolds numbers tested (Re = 0.047, 2, 5, 10 not shown). The similarity of the metachronal cases (i.e. ∆Φ = 0) for this range of Reynolds numbers is certainly due to the fact that there are always cilia in the stroke phase. Hence, a transport is observed during all the beat cycle, and the inertial effects thus remain minor. This is no longer the case for fully synchronously beating cilia: despite being weak, inertial effects cancel the reversal of the flow.

Conclusions

By using a lattice Boltzmann solver coupled to an Immersed Boundary method, and considering a purely hydrodynamical feedback from the fluid, symplectic and antipleptic metachronal waves can emerge in a 3D two-phase flow configuration with a viscosity ratio of 20. It is known that metachronal waves may emerge due to hydrodynamic interactions [START_REF] Elgeti | Emergence of metachronal waves in cilia arrays[END_REF]. However, for the first time, both types of metachrony are observed to emerge in a two-layer fluid with different viscosities, using a simple feedback law. This feedback depends on a coupling parameter α that can be used to tune the strength and direction of the waves. It is observed that cilia experience weaker torques when being in a clustered configuration, i.e. during the recovery phase of antipleptic motion, as well as during the stroke phase of symplectic motion. The resulting beating pattern is a slow stroke phase and a fast recovery phase for both cases, suggesting that even if a simple hydrodynamic interaction is enough to let emerge metachronal waves, it is not sufficient to reproduce all features of the beating patterns observed in real ciliated surfaces, and other biological issues may play a role in the transport mechanism. The study of [START_REF] Hussong | A continuum model for flow induced by metachronal coordination between beating cilia[END_REF] has shown that metachronal motion can switch from symplectic to antipleptic with increasing inertial effects. Here, the influence of cilia density is highlighted: by assuming a least effort behaviour for the cilia, the present model shows that the metachrony can switch from antipleptic to symplectic by lowering the cilia density. A more detailed study of the quality of synchronization along the x and y-directions is the next step toward a better understanding of the emergence of MCW.

A thorough comparative study of the antipleptic and symplectic MCW has been performed and the results show that the antipleptic MCW are the most efficient ones for transporting mucus. This is in accordance with most recent studies addressing this point [START_REF] Khaderi | Microfluidic propulsion by the metachronal beating of magnetic artificial cilia: a numerical analysis[END_REF][START_REF] Osterman | Finding the ciliary beating pattern with optimal efficiency[END_REF][START_REF] Ding | Mixing and transport by ciliary carpets: a numerical study[END_REF]; and now confirmed in a two-layer environment. In the range of cilia spacing studied, and especially for small phase lags (|∆Φ| < π/2) the antipleptic MCW have a better ability to (i) transport the flow in the same direction compared to symplectic MCW with the same wavelength, (ii) generate higher flow rate, (iii) advect particles at a given power input, and (iv) generate a higher stretching and hence, are more able to mix fluids. On the contrary, symplectic MCW do not appear to have a great impact on the flow and are often less efficient than antipleptic MCW. This is not in agreement with [START_REF] Elgeti | Emergence of metachronal waves in cilia arrays[END_REF] who reported that symplectic waves are almost nearly as efficient as antipleptic waves, using an optimized-efficiency model. This difference may be due to several factors, including the two-layer flow character of the present work, the ratios of viscosity considered, or even the cilia beat shape used. However, these results are in accordance with [START_REF] Khaderi | Microfluidic propulsion by the metachronal beating of magnetic artificial cilia: a numerical analysis[END_REF] who explained this better transport efficiency as being due to vortices obstructing the flow for symplectic motions.

Among the quantities introduced earlier, some are more decisive for the transport of particles. The most important ones are the net volume of mucus displaced, the transport efficiency η 1 , and the mixing capacity of the system. It has been shown all along this work that antipleptic MCW with ∆Φ = π/4 generate the larger flow while in the meantime showing the higher value of the transport efficiency (the power spent P * for ∆Φ = π/4 being minimal). A preliminary study of the mixing also strongly indicates that they might also be the more suitable to mix fluids. The best transport efficiency of the antipleptic MCW is certainly due to the clusterized aspect of the cilia in recovery phase which minimizes their impact on the flow, while cilia in stroke phase are able to fully exert their pushing action. While being more able to advect particles than any other form of coordinated motions, antipleptic MCW with ∆Φ < π/2 also require less power. Note that the results presented all along concern a viscosity ratio r ν of 20. [START_REF] Chatelin | A parametric study of mucociliary transport by numerical simulations of 3D non-homogeneous mucus[END_REF] has shown that mucus transport was maximized for viscosity ratios ranging from 10 to 20. With the present model, and within this interval of viscosity ratios, the same trends are always observed (results not shown), and the best transport capacities are always obtained for antipleptic MCW with ∆Φ = π/4. The transport properties of the waves also remain similar and do not vary significantly. However, it is worth noticing that the best transport capacity was obtained for r ν = 20. Hence, the optimum is probably outside the range of viscosity ratios studied. It also strongly indicates that increasing the mucus viscosity results in an increase of the transport properties. As in previous studies [START_REF] Norton | Model of ciliary clearance and the role of mucus rheology[END_REF][START_REF] Chatelin | A parametric study of mucociliary transport by numerical simulations of 3D non-homogeneous mucus[END_REF], an optimum viscosity for the transport of particles can certainly be found. A more detailed study of how different mucus viscosity affect the transport properties of the MCW is the next step toward a better understanding of the mucociliary clearance.

While in the present study only one prescribed beating pattern was investigated, the results can certainly be generalised to other beating patterns. An interesting path of investigation would be to study how this beat shape is influenced by different mucus viscosities, mucus rheologies, or thicknesses of the PCL. Note that a Reynolds number of 20 was used for computational reasons. While this does not affect the behaviour of both the antipleptic and symplectic MCW, weak inertial effects are seen to change the dynamics of the flow in the particular case of synchronously beating cilia. Finally, the mucus was considered as being a Newtonian fluid in the present work, while it is in reality highly non-Newtonian. It is expected that the clearance velocity of both the antipleptic and symplectic MCW would be increased by considering the mucus as being a viscoelastic fluid. One could also expect the symplectic waves to become more efficient than the antipleptic MCW for highly viscous mucus. Indeed, the clusterized aspect of the cilia during the stroke phase of symplectic motion could help them overcome more easily the viscous resistance of the mucus (Knight-Jones 1954). Nevertheless, note that studies [START_REF] Norton | Model of ciliary clearance and the role of mucus rheology[END_REF][START_REF] Chatelin | A parametric study of mucociliary transport by numerical simulations of 3D non-homogeneous mucus[END_REF] all tend to show, however, that the mucus can not be too "solid" nor too "liquid".

The numerical results obtained here constitute a useful basis of investigation to progress in the understanding of respiratory diseases linked to cilia beating disorders, such as COPD or severe asthma [START_REF] Chanez | Severe asthma is an epithelial disease[END_REF]. It is also of interest for industrial purposes, such as the design of cilia-based actuators for mixing [START_REF] Chen | Magnetically actuated artificial cilia for optimum mixing performance in microfluidics[END_REF] or flow-regulator in microscopic biosensors.

Future works include the implementation of a realistic rheological model [START_REF] Lafforgue | Rheological characterization of macromolecular colloidal gels as simulant of bronchial mucus[END_REF] for the mucus in order to take into account its highly non-Newtonian behaviour, and the study of the mass transfer occurring at the epithelial surface. A long term objective of this work is to build a numerical environment to predict the transport and mixing of drugs inside both the mucus and PCL layers. With this tool, the effect of drugs could be tested virtually before proceeding to clinical trials. For instance, drugs acting on microscopic parameters such as beating frequency, viscosity of the mucus, density of cilia could be tested and their effects on macroscopic quantities such as mucus flow rate and mixing could be explored, in order to progress in the understanding and treatment of respiratory diseases.

The authors would like to thank Dr. Zhe Li and master intern Jean Mercat for the development of the numerical code. S. Chateau and S. Poncet acknowledge also the Natural Sciences and Engineering Research Council of Canada for its financial support through a Discovery Grant (RGPIN-2015-06512). This work was granted access to the HPC resources of Compute Canada and Aix-Marseille University (project Equip@Meso ANR-10-EQPX-29-01).

Figure 1 .

 1 Figure 1. (a) Schematic view of the computational domain. The present case corresponds to an antipleptic MCW. The domain is filled with PCL (in blue) and mucus (in red). (b) Beating pattern of a cilium with the parametric equation used. Steps 1 to 6 correspond to the recovery phase, and steps 7 to 9 to the stroke phase.

Figure 2 .

 2 Figure 2. Schematic view of a cilium with the corresponding forces exerted on the fluids.

Figure 3 .

 3 Figure 3. Domain filled with passive tracers. The present case corresponds to an antipleptic MCW with ∆Φ = π/4. An array of 8 cilia with a cilia spacing a/L = 1.67 is considered on a computational domain of size (Nx = 201, Ny = 26, Nz = 50). The mucus phase (in red) is located above the PCL phase (in blue), and the ratio of viscosity is set such as rν = 10. The colorbar represents the dimensionless velocity magnitude of the fluids on the periodic boundary over the x-direction.

Figure 4 .

 4 Figure4. Symplectic MCW emerging from an initially random state of the cilia. 24 cilia arranged in a 8 × 3 rectangle are considered on a computational domain of size (Nx = 361, Ny = 136, Nz = 50) with a cilia spacing a/L = 3. The mucus phase is in red and the PCL phase in blue. The color bar indicates the phase of a particular cilium within one beating period, which is represented by a circle at its base. Top: 3D view of the system. Bottom: 2D view of the same system in a (x,y) plane to highlight the 3D modulation in the z-direction.

Figure 5 .

 5 Figure5. Antipleptic MCW emerging from an initially random state of the cilia. 1024 cilia arranged in a 32 × 32 square are considered on a computational domain of size (Nx = 161, Ny = 161, Nz = 32) with a cilia spacing a/L = 0.23. The mucus phase is in red and the PCL phase in blue. The color bar indicates the phase of a particular cilium within one beating period, which is represented by a circle at its base. Top: 3D view of the system. Bottom: 2D view of the same system in a (x,y) plane to highlight the 3D modulation in the z-direction.

Figure 6 .

 6 Figure 6. Different kinds of collective coordination for the beating cilia. 32 cilia arranged in a 8 × 4 square are considered on a computational domain of size (Nx = 241, Ny = 121, Nz = 50) in each case, with a cilia spacing a/L=2. (a) Synchronized motion; (b) Symplectic metachronal motion; (c) Antipleptic metachronal motion; (d) No synchronization (random state of cilia). The figures show contours of the magnitude of the dimensionless vorticity ω using a logarithmic scale. The black lines show the frontier between the PCL at the bottom and the mucus layer above.

Figure 7 .

 7 Figure 7. Normalized average displacement in the x-direction as a function of z/L for (a) ∆Φ = ±π/9 and ∆Φ = ±π/4 with a/L = 0.8; (b) ∆Φ = ±π/4 and a/L = 1.67. As in Ding et al. (2014), the mixing zone extends from 0 to 1L, and the transport zone from 1L to 1.4L.In both cases (a) and (b), the mucus-PCL interface is located at z/L = 0.9 and is indicated by a dashed horizontal line.

Figure 8 .

 8 Figure 8. Stretching rate and direction for: (a) Synchronized motion (∆Φ = 0), (b) Symplectic MCW (∆Φ = -π/4), (c) Antipleptic MCW (∆Φ = π/4). In each case, the cilia spacing is a/L = 1.67, and the size of the computational domain is (Nx = 201,Ny = 26,Nz = 50).

Figure 9 .

 9 Figure 9. Directional pushing efficiency ǫP N over one beating cycle for three different phase lags ∆Φ. Results obtained for an antipleptic MCW and a/L = 2.

Figure 10 .Figure 11 .

 1011 Figure 10. Snapshots of the fluid velocity taken at 8 different instants during the beating cycle of the cilium located at the left of the images. The cilia are beating in an antipleptic motion with a phase lag ∆Φ = π/4, and a spacing of a/L = 2. (a) t = 0: the recovery phase of the first cilium begins; (b-e) recovery phase; (f-h) stroke phase. Once (h) is completed, a new cycle begins: (h) → (a). The figures show contours of the norm of the normalized fluid velocity V f such that: V f = V dim f V ref f with V ref = λ/T the reference speed of the present system.

Figure 12 .

 12 Figure 12. Normalized volumetric flow rate Qv as a function of time over a beating cycle for different phase lags ∆Φ and different cilia spacings a/L. The recovery phase occurs for t ∈ [0, π] and the stroke phase for t ∈ [π, 2π].

Figure 13 .Figure 14 .

 1314 Figure 13. Total dimensionless displaced flow volume generated by an array of cilia over a beating cycle for different phase lags and cilia spacings. +: a/L = 1.67;: a/L = 2; * : a/L = 2.53; : a/L = 3.33.

Figure 15 .

 15 Figure 15. Comparison of the flow generated by an antipleptic and a symplectic wave for a the cilia spacing a/L = 1.67. (a): Antipleptic MCW with ∆Φ = 2π/3. (b): Symplectic MCW with ∆Φ = -2π/3. The plane is colored with the magnitude of the dimensionless fluid velocity.

Figure 16 .

 16 Figure 16. Displacement ratio η1 as a function of the phase lag ∆Φ for different cilia spacings a/L.

Figure 17 .

 17 Figure 17. Average stretching rate in the transport and mixing areas as a function of the phase lag ∆Φ for different cilia spacings. (a) +: a/L = 1.67; : a/L = 2; * : a/L = 2.53; : a/L = 3.33 (b) 3D view of the corresponding plot.

Figure 18 .

 18 Figure 18. Total displaced volume of fluid for Re = 0.02, Re = 1 and Re = 20 as a function of the phase lag ∆Φ. The results of these simulations are obtained with L = 7 lu and a/L = 1.67.

† Email address for correspondence: sylvain.chateau@usherbrooke.ca